Identification of circRNA–miRNA– mRNA regulatory network associated to the autism spectrum disorder in children through integrated bioinformatics analysis

Background Autism spectrum disorder (ASD) is a complex neurological disability with multifactorial etiology. ASD is described by behavior, speech, language, and communication defects. CircRNA is a type of ceRNA that plays an important role in modulating microRNAs (miRNA) in several disorders. However, the potential role of the circRNA/ miRNA/mRNA regulatory network in the pathogenesis of ASD is not fully understood. Therefore, this study aimed to create a circRNA/miRNA/mRNA network associated with ASD to cast light on the pathogenesis of ASD. Methods CircRNA expression profile data were recruited from Gene Expression Omnibus datasets, and the differentially expressed circRNAs (DEcircRNAs) were identified. Then, miRNAs modulated by these circRNAs were predicted and overlapped with differentially expressed miRNAs. Next, the potentially involved genes were identified by overlapping predicted targets, and differentially expressed genes. The enrichment analysis was performed, and a PPI network was projected. Subsequently, ten key genes were selected from the network. Furthermore, a circRNA/miRNA/ mRNA regulatory network was constructed, and probable molecules and drugs with potential anti-ASD effects were predicted. Results 11 DEcircRNAs and 8 miRNAs regulated by 4 circRNAs were identified as being significantly involved. Subse-quently, gene enrichment analysis of 71 overlapped mRNA regulated by these miRNAs showed that they are mostly associated with hippocampal synaptogenesis, neurogenesis, and axon guidance. Additionally, two high-score compounds, GSK3β inhibitor (SB216763) and dexamethasone, and three drugs (haloperidol, nystatin, paroxetine) were confirmed as potential therapeutic options for ASD. Conclusion The results of this study may help gain deeper insight into the pathogenesis of the circRNA/miRNA/ mRNA regulatory network in ASD, providing potential therapeutic management options


Introduction
Autism spectrum disorder (ASD), commonly known as autism, is a common, highly heritable, heterogeneous neurodevelopmental disorder with underlying cognitive properties [1].The incidence of this complex neurobehavioral disorder has increased vastly in recent years.This disease is not limited to race, ethnicity, or specific geographical regions.Autism is 4 times more prevalent in boys than girls [2,3].Based on the Diagnostic and Statistical Manual of Mental Disorders (DSM), autistic patients are classified into four groups: Asperger's syndrome (AS) in which the affected individuals show the mildest symptoms, have normal intelligence, and can live a relatively normal social life; Pervasive developmental disorder in which the affected individuals display more severe symptoms; Real autistic patients lie almost in the middle of the spectrum; Dissociative Identity Disorder (DID) in childhood in which the affected individuals have the most severe symptoms.These children have normal development until the age of 2, however they later undergo a sudden significant regression in many aspects, including social communication, speech, and learning skills [4,5].ASD is characterized by altered early brain development and neural reorganization.However, due to the lack of suitable biomarkers, diagnosis is based on behavioral observations.Non-coding RNAs (ncRNAs) and related genes have been suggested as beneficial biomarkers in the diagnosis of many neurodevelopmental disorders [6][7][8].The ncRNA transcripts are highly conserved and involved in the regulation of gene expression [9].Based on their length, ncRNAs are classified into two classes, including small noncoding RNAs (sncRNA) and long noncoding RNAs (lncRNAs).microRNAs (miRNAs) and circular RNAs (circRNAs) belong to the first and second groups, respectively [10,11].miRNAs are small noncoding sequences (mostly 19-24 nucleotides in length) that negatively regulate gene expression after transcription through binding to the target mRNAs [12].Many studies have proven that Some miRNAs are involved in ASD etiology, and altered level of these molecules in the cerebral neurons leads to changes in the expression of the target genes [13].CircRNAs are a type of non-coding RNA with a circular structure and a length of more than 200 nucleotides.Numerous studies have shown that circRNAs are vital gene regulators expressed widely in different tissues and specific developmental stages [14].Some of these RNA molecules are highly conserved among different species.Binding to miRNA and acting as a sponge, circRNAs play a significant role in the regulation of miRNAs' target genes [15,16].CircRNAs are highly expressed in the mammalian brain compared to other tissues.These molecules are also involved in various processes, including neurogenesis, synaptogenesis, and neuronal development [17].
Thus, these ncRNAs are concluded to play functional roles in neuron development and related diseases.In the present study, some circRNAs, miRNAs, involved genes, and their potential regulatory mechanisms in ASD have been studied using bioinformatics analysis.Figure 1 presents the flowchart of work.First, circRNAs, miRNAs, and differentially expressed genes were obtained from the Gene Expression Omnibus (GEO) database.After identifying the role of each group, enrichment analysis, and the Regulatory Network Enrichment Analysis tool were used to find the key genes in these pathways.Finally, the interactions of drugs with the hub genes were investigated to suggest possible therapeutic choices for ASD.

Differentially expression identification
Limma, a Bioconductor package in R software, was implemented to identify the differentially expressed circRNAs (DEcircRNAs), miRNAs (DEmiRNAs), and Genes (DEGs) between ASD and healthy control.The log-fold change (FC) in expression and p value (P) were calculated.The p value < 0.05 and |logFC|> 2.0 were considered statistically significant for DEcircRNAs, and p value < 0.05 and |logFC|> 0.5 were considered for DEmiRNAs and DEGs, applying the Benjamini-Hochberg (BH) technique [18].The overlapping miRNA and genes were examined using the Venn diagram web tool.Visual hierarchical cluster analysis was also conducted to present the volcano plots of DEcircRNAs, DEmiRNAs, and DEGs.

Prediction of miRNA response elements
The Circular RNA Interactome (CircInteractome) and cancer-specific circRNAs database (CSCD, https:// gb.whu.edu.cn/ CSCD) were used to predict the target miRNA response element (MRE) sites and RNA binding protein (RBP) sites for each DEcircRNAs.DIANA-miRPath v3.0 (https:// www.micro rna.gr/ miRPa thv3) is an online database to assess miRNAs regulatory roles and estimate the related regulation pathways.Overlapped miRNAs of the DEmiRNAs and miRNAs in MRE positions were projected as potential target miRNAs of the DEcircRNAs.DIANA-miRPath also predicted the function of the overlapped miRNA.The overlapped miRNAs were then selected for further target predictions.

Prediction of miRNA targets
miRNA-gene interactions were predicted using Integrated Target Prediction of miRNA MiRWalk 3.0 software (http:// mirwa lk.umm.uni-heide lberg.de/), miRmap [19], and TargetScan [20], which identify miRNAs binding sites within the complete gene sequence.All overlapped genes in three databases were selected and intersected with the DEGs to obtain candidate target genes.

The circRNA-miRNA-mRNA regulatory network
DEcircRNAs, interactions DEmiRNAs, and interactions DEGs were selected to establish a circRNA-miRNA-mRNA regulatory network.Cytoscape 3.9.1 version was used to present this network.

Enrichment analysis
The Metascape database (http:// metas cape.org/) and R packages (clusterProfiler, org.Hs.eg.db,AnnotationDbi) were used to perform Gene Ontology (GO) analysis including biological process (BP), molecular function (MF), cellular component (CC), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the overlapped genes in the network.Metascape is a freely accessed web-based online bioinformatics resource with useful tools for the functional interpretation of large lists of genes including GO functional, KEGG pathways, and Hallmark Gene Sets.

The protein-protein interaction network
Protein-protein interaction (PPI) analysis was performed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) for genes in the regulatory network and visualized by Cytoscape 3.9.1.PPI is related to the interactions between proteins in a biological process.Then, the network was used for topological analysis.The CytoHubba plugin with the degree method in Cytoscape software was used to detect the top ten hub genes in the network.CytoHubba offers 11 approaches for the topological analysis of networks from different viewpoints.In addition, critical modules consisting of hub genes and several connective modules in the PPI networks were revealed by the molecular complex detection (MCODE) plugin (degree cutoff = 2, node score = 0.2, Max depth = 100).

Connectivity map (CMap) analysis
Connectivity map (CMap) analysis is a gene expression profiling database in which gene expression is connected to diseases.CMap facilitates the process of comparing small molecule compounds or drugs that are significntly related to diseases using gene expression profiles.CMap analysis according to dysregulated genes in the network was conducted to reveal candidate drugs for ASD treatment.The connectivity scores (− 100 to 100) were adjusted to indicate closeness among the genes and compounds; a positive score represented a positive connection, and a negative score indicated a negative connection with the uploaded genes.

Drug-ASD gene interaction analysis
Drug-gene interaction analysis of hub genes was carried out using DGIdb (https:// dgidb.org/) [21].DGIdb is a web resource that provides information about genes correlated with approved or potential drugs.Key significant genes uploaded to DGIdb for equivalent with FDA-approved drugs to identify the possible targets for ASD treatment.

Identification of circRNA-miRNA interactions
CircRNA can play pivotal roles in many disorders by functioning as a "sponge" to miRNAs.Therefore, some affected miRNAs were predicted based on this ceRNA theory.The structures of the 11 circRNAs are shown in Fig. 3, all of which have MRE, RBP, and ORF.Two online databases, CSCD and CircInteractome, were used to predict potential miRNAs targeted.A total of 658 MREs were identified for DEcircRNAs, and 8 significant interacting miRNAs were obtained by overlapping 239 MREs with 20 DEmiRNAs.Then, 4 circRNA-miRNA interactions, including hsa_circ_0053004-miR-1276/miR-940, hsa_circ_0125982/miR-330, hsa_circ_0088217-miR-330, and hsa_circ_0129247-miR-564/miR-330/miR-142/ miR144/miR-577/miR-921 were identified.Moreover, DIANA-miRPath was used to investigate the signaling pathways in which 14 miRNAs may be involved.According to Fig. 4, these miRNAs are involved in the GABAergic synapses, thyroid hormone signaling pathway, FoxO signaling pathway, and TGF-beta signaling pathway.

Identification of overlapped genes
The target prediction was done for 8 miRNAs associated with 4 DEcircRNAs to investigate how miRNAs are involved in autism disorder.In total, 5930 target mRNAs were identified using TargetScan and miRDB databases.Then, 71 overlapping interacting genes were selected by the intersection of 216 DEGs from gene expression data and target genes.

Functional and pathway enrichment analyses
To better explore the potential mechanisms of genes in the ceRNA regulatory network, KEGG pathway enrichment analysis and GO analysis were accomplished, which comprises biological process (BP), molecular function (MF), and cellular component (CC).The top 10 highly enriched GO and KEGG pathways by Metascape database are shown in Fig. 6A, and the top 20 highly enriched GO and KEGG terms analyzed by R software are presented in Fig. 6B.GO and KEGG analysis showed that the interacting genes were involved in transcriptional regulation by MECP2, hippocampal synaptogenesis and neurogenesis, cell junction assembly, neuron projection development, commissural neuron axon guidance, and axon guidance.GO analysis showed that, in the case of BP, these overlapped genes were mainly involved in peptidyl serine phosphorylation, regulation of metal ion transport, neuron projection development, sympathetic ganglion development, and brain development.For MF, overlapped genes were mainly enriched in transmembrane receptor protein tyrosine kinase activator activity, commissural neuron axon guidance, regulation of mitotic cell cycle, calmodulin-dependent protein kinase activity, and cell junction assembly.For CC, neuron projection, endocytic vesicle, intracellular membrane-bounded organelle, ciliary landscape, and axon were predicted.

Construction of the protein-protein interaction and hub genes identification
The NetworkAnalyst tool linked with the STRING database was used to construct a PPI network based on the 86 DEGs (Fig. 7A).The PPI network contained 984 nodes and 2607 edges.Then, the cytoHubba plugin in Cytoscape 3.9.1 was used to screen the hub genes.Ten genes (CAMK2A, CAMK2D, CAMK2G, MECP2, CAMK 4, NRG1, DLG1, BDNF, MAPK1, and PRKCA) were identified as hub genes in this cluster, including ten nodes and 26 edges (Fig. 7B).

Candidate compounds from connectivity map analysis
The regulatory networks of circRNAs/miRNAs/mRNAs were constructed, and dysregulated genes in the network were uploaded to the CMap website.Then, compounds with potential influence on ASD were predicted.GSK3β inhibitor (SB216763), with the highest positive scores, was considered a potential drug for ASD treatment.

Pharmacogenomics analysis for hub genes to find potential drug
The DGIdb website was used to perform drug-gene interaction analysis for the ten hub genes identified in the PPI network.Finally, it was found that Haloperidol, nystatin, and paroxetine are potential drugs and might be effective in treating ASD since they were closely associated with the four hub genes (CAMK2G, BDNF, NRG, MAPK1).

Discussion
Autism fundamentally affects the patient's quality of life with growing prevalence and incidence worldwide [22].Although many studies have been conducted on autism in recent years, the pathogenesis of autism is still unknown.The recent developments in bioinformatics facilitated more efficient research on pathogenesis mechanisms and treatment strategies [23].
In this study, we have evaluated the regulatory network between circRNA, miRNA, and their associating genes using the GEO data sample of children with autism.Cir-cRNAs, as non-coding RNAs, acting like ceRNAs, can affect various cellular biological processes [24].Also, various computational algorithms, such as inductive matrix completion, Locality-Constrained Linear Coding algorithm, and KATZ algorithm have proven that circR-NAs and miRNAs are strongly correlated with the pathogenesis of several diseases [25][26][27][28].Therefore, circRNAs and miRNAs can be used as diagnostic biomarkers and help improve novel treatments.In the current study, eleven circRNAs were selected, including 3 upregulated and 8 down-regulated circRNAs in ASD patients.To our knowledge, DEcircRNAs have not been studied yet, but their expression is significantly changed in autistic patients.Thus, they may be excellent diagnostic biomarkers or therapeutic targets in ASD.One of the most prominent features of circRNAs is that they act as a sponge to absorb associating miRNAs, using their ability to interact with miRNA binding sites.Therefore, they can indirectly regulate gene expression.Among these 11 circRNAs, 4 of them can regulate the expression of 8 miRNAs based on the ceRNA theory.Surprisingly, overexpressed miRNAs are associated with downregulated circRNAs, and the underexpressed miRNAs have binding sites on upregulated circRNAs.Limited studies have investigated a correlation between these miRNAs and neural development and neuronal-related disorders.Mor et al. showed that the expression of miR-142-5p and miR-144 is remarkably increased in the brain tissue of autistic patients.The dysregulation of these miRNAs can affect the expression of other genes that control the synaptic processes.Also, miR-496 can regulate neurogenesis and neuronal migration through fine-tuning of N-cadherin levels [29], and miR-921 was significantly increased in a microarray analysis accomplished on Chinese children with ASD [30].Go and KEEG analysis have been carried out on the overlap genes between the DEGs and miRNAs target genes.One of the suggested pathways is the axon guidance pathway, which was previously demonstrated in autism [31].Studies have shown poor expression of axon guidance proteins, such as PLXN4 and ROBO2, in the brain of autistic patients [32].Cell junction assembly pathway in synaptogenesis was also projected by GO analysis.The synaptogenesis step is vital in brain development in children.In ASD patients, this process is interrupted by the mutation or altered expression of associating genes, such as SHANK3 and CDH10 [33].
After studying the pathways involved in the regulatory networks of circRNA-miRNA-mRNA, ten hub genes were selected with high scores, five of which are calcium/calmodulin-dependent protein kinase (CAM-KII) (encoded by CAMK2A, CAMK2B, CAMK2G, and CAMK2D genes).There is growing evidence regarding the prominent regulatory role of CAMKII in the synaptic processes and synaptic plasticity, which can affect different functions, such as learning and memorizing.Numerous mutations in CAMK2A, CAMK2B, and CAMK2G genes have been identified that correlate with neural development disorders and neurological disorders, such as ASD [34].Another significant protein is MECP2, which is vital for neural development and regulates the expression of BDNF protein.MECP2 plays a critical role in synaptic homeostatic plasticity, and the absence of the encoding gene caused Rett syndrome in animal models, which is classified as a severe form of autism spectrum disorder [35].In addition, in human samples, the expression of MECP2 protein was significantly decreased.MECP2 down-regulation is due to various mutations or methylation in the promotor region.Another key protein, neuregulin 1 (NRG1), is vital for central nervous system development and activation of neurotransmitter receptors, such as glutamate.NRG1 protein is involved in neuronal migration, synaptogenesis, glycogenesis, dendritic growth, and communications between neurons and glial cells.Neuronal activities specifically regulate the expression of NRG1 protein, leading to significantly reduced expression in autistic children [36].Two bioactive compounds, dexamethasone, and GSK3β inhibitors were identified for treatment strategies according to the CMap analysis.The GSK3β inhibitor disturbs the function of the GSK3 protein in different signaling pathways, including the neuronal signaling pathways, and can improve ASD symptoms [37].Haloperidol is a classic antipsychotic that is used to improve symptoms of schizophrenia.Haloperidol is an FDA-approved drug for the treatment of schizophrenia, severe behavioral disorders in children, Tourette syndrome, and hyperactivity in children.Haloperidol interacts with two of the hub genes and can reduce autistic symptoms in children (names of genes) [38,39].Nystatin, a compound frequently used to preserve intracellular dialyzable components, specifically inhibits the potassium channel Kv1.3 [40].Paroxetine is a selective serotonin reuptake inhibitor (SSRI) used to treat behavioral disorders including panic disorder, major depressive disorder, posttraumatic stress, social anxiety disorder, and generalized anxiety disorder [41,42].

Conclusion
In this study, two regulatory networks of circRNAs and their functions was investigated by bioinformatics analysis.In addition, some compounds and drugs were predicted which may have potential use in the management of ASD and improving symptoms.Such studies can broaden our knowledge of ASD pathogenesis and guide future research.Finally, the results of the present study are based on bioinformatics analyses and require further investigation.

Fig. 1
Fig. 1 Flowchart of current study about creating a circRNA/miRNA/mRNA network and forecasting potential therapeutic management for ASD

Fig. 2
Fig. 2 Results of differential analysis shown by volcano plot for gene, miRNA, and circRNA.A Volcano plot showing fold differences in circRNA expression and p value.B Fold differences in genes.And C for miRNAs.Red points represent dysregulated expression

Fig. 3 Fig. 4 Fig. 5
Fig. 3 structures of the 11 DEcircRNAs.The various colors in the outer and inner ring represent the diverse exons and the positions of MRE, RBP, and ORF

Fig. 6
Fig. 6 GO enrichment analysis and KEGG enrichment analysis of differential genes in ASD and control groups.A The top 20 GO and pathway enrichment analyses of biological processes, cellular components, and molecular functions by Metascape database.B KEGG and GO enrichment analysis of differential genes in the ASD conducted by clusterProfiler package R software

Fig. 7 A
Fig. 7 A The PPI network of overlapped interacting genes.Interacting DEGs are shown by the nodes, and interactions between genes by the edges.B Top ten hub genes constructed thereby Cytohubba plugin.The network shows the degree of significance of hubs through a color scale ranging from red (the most important) to yellow (less important)