Investigation of Gene Regulatory Networks Associated with Autism Spectrum Disorder Based on MiRNA Expression in China

Autism spectrum disorder (ASD) comprise a group of neurodevelopmental disorders characterized by deficits in social and communication capacities and repetitive behaviors. Increasing neuroscientific evidence indicates that the neuropathology of ASD is widespread and involves epigenetic regulation in the brain. Differentially expressed miRNAs in the peripheral blood from autism patients were identified by high-throughput miRNA microarray analyses. Five of these miRNAs were confirmed through quantitative reverse transcription-polymerase chain reaction (qRT-PCR) analysis. A search for candidate target genes of the five confirmed miRNAs was performed through a Kyoto encyclopedia of genes and genomes (KEGG) biological pathways and Gene Ontology enrichment analysis of gene function to identify gene regulatory networks. To the best of our knowledge, this study provides the first global miRNA expression profile of ASD in China. The differentially expressed miR-34b may potentially explain the higher percentage of male ASD patients, and the aberrantly expressed miR-103a-3p may contribute to the abnormal ubiquitin-mediated proteolysis observed in ASD.


Introduction
Autism spectrum disorder (ASD) comprise a group of chronic neurodevelopmental disorders characterized by social and language impairments and restricted and repetitive interests and behaviors [1]. ASD includes autistic disorder, Asperger's syndrome, and pervasive developmental disorder not otherwise specified (PDD-NOS) [2]. Patients with ASD vary greatly in terms of clinical manifestation and may also show associated medical comorbidities. Symptoms begin to appear at the age of three years, and affected individuals often require constant care from family members or professionals [3][4][5]. The Autism and Developmental Disabilities Monitoring (ADDM) Network calculated the prevalence of 8-year-old ASD as 14.7% based on the data from 11 ADDM Network sites in the USA obtained in 2010 [6]. The prevalence rates have increased significantly in the past decade, which may be due to the comprehensive effect of prenatal risk factors, environmental pollution, and heritable factors. It is one of the most heritable of the common neurodevelopmental diseases. Until now, it has been difficult to understand how diverse genetic susceptibilities translate to a clinical phenotype with so many genomic loci contributing to heterogeneous functions [7]. Thus far, studies have focused more on genetic factors, but only a few studies have investigated the role of miRNA in autism spectrum disorder [8][9][10].
MicroRNAs (miRNAs) are small endogenous non-coding regulatory RNAs (typically 21-23 nucleotides) that function as posttranscriptional regulators of gene expression [11,12]. These are known to play a critical role in neurodevelopment, metabolism, neuroplasticity, apoptosis, and other fundamental neurobiological processes. By complementarily base-pairing with the 3'untranslated region (3'UTR) of specific target mRNAs, miRNA can regulate gene expression [13,14]. Recent studies have demonstrated that miRNAs are present in human body fluids, such as serum, plasma, and cerebrospinal fluid. [15][16][17][18]. Several recent studies have reported differential miRNA expression in autism spectrum disorder. However, no specific study of regulatory miRNA expression and the gene regulatory networks of aberrant miRNAs has focused on autism spectrum disorder in China. Therefore, we performed the first investigation of the miRNA expression profile in the peripheral blood of Chinese autism patients and the gene regulatory networks of aberrant miRNAs.

Materials and Methods
Sample collection and RNA purification 5 autism patients (1 female and 4 males; mean age: 4.9±1.917) and 5 controls were recruited in microarray analysis study and another 15 patients (2 females and 13 males; mean age: 4.3± 1.623) and 15 controls were recruited in qRT-PCR analysis study. The diagnosis of ASD was established in the patients with the use of the Diagnostic and Statistical Manual, Fourth Edition, Text Revision (DSM-IV-TR; American Psychiatric Association, 2000) criteria. Patients show typical ASD clinical symptoms, including deficits in communication and restricted patterns of interests or activities, while social deficits are observed in the school-age patients. Peripheral blood samples were obtained from participants (5ml blood for each) at Xiangya Hospital, and written informed consents were obtained from the parents of all of the participants because the participants were underage and did not have the capability to provide signed written informed consent. This study was approved by the Institutional Review Board and the Ethics Committee at Xiangya Hospital, Central South University (Changsha, China).
The total RNA was isolated using the Trizol method (Life Technologies, USA) according to the manufacturer's instructions. The quality of RNA was examined using an Agilent Bioanalyzer (Santa Clara, CA, USA). The RNA quantity and purity were assessed using a K5500 microspectrophotometer (Kaiao, China). Values of A260/A281.5 and A260/A2301 indicate acceptable RNA purity, and a value of RIN (RNA Integrity Number)7 obtained through the Agilent 2200 RNA assay indicates acceptable RNA integrity (Agilent, USA). Genomic DNA contamination was evaluated by gel electrophoresis. The total RNA was stored at -80°C until use.

miRNA expression profiling by microarray
RiboArray miDETECT Human Array (A10101-1-12-19, 1×12K) microarrays were used to screen the miRNA expression profile (RIBOBIO, China). The microarrays contained 2,578 assay probes corresponding to the entire set of annotated human and nonhuman primate miRNA sequences from miRBase20.0. The internal control of the microarray assay was from RiboArray internal controls database referenced from multiple public probe database including Exiqon miRCURY LNA, GeneChip miRNA Array and Agilent miRNA. In this study, five autism patients aged 2.5 to 7 years of age and five age-matched controls were enrolled, and all of these were used as unique samples. For hybridization, 2.5 μg of total RNA labeled with Cy5 was used for each sample. To yield adequate microarray results, the labeling efficiency should be between 1.0 and 3.6. After overnight hybridization and washing, the fluorescence images were scanned using a GenePix 4000B laser scanner (Molecular Device, USA) and digitized using the R software. After normalization using the quantile normalization method [19], the p-values were calculated with the Rank Product Method [20]. miRNAs with p<0.05 were selected for cluster analysis. The raw data is available at Gene Expression Omnibus (accession number: GSE67979).

miRNA qRT-PCR analysis
The total RNA was reverse transcribed to cDNA in a volume of 20 μL according to the manufacturer's instructions. The five most differentially expressed miRNAs in the 15 autism patients and 15 health controls were screened for further validation by Bulge-Loop miRNA qRT-PCR assays (CFX96, Bio-Rad, USA). Has-miR-16-5p was used as an internal reference in the qRT-PCR experiments [15]. The miRNA qRT-PCR Primer Set (RiboBio Co. Guangzhou, China) was used to detect and quantify the expression of the five microRNAs (let-7a-5p, let-7d-5p, miR-34b-3p, miR-103a-3p, and miR-1228-3p) according to the manufacturer's instructions [21]. Every qRT-PCR assay was performed in triplicate in a volume of 20 μL containing 1 μL of cDNA. The miRNA expression level was calculated using the 2 -ΔΔCt method [11].

Investigation of gene regulatory networks based on miRNA expression
To validate the five most differentially expressed miRNAs validated by qRT-PCR analysis, several specific algorithms have been designed to predict the target genes of miRNAs. We first used four publicly available databases, namely TargetScan, miRanda, CLIP-Seq, and miRDB, for the prediction of potential target genes. Only those candidates that appeared in at least three databases were cinsidered candidate miRNA target genes. A Kyoto encyclopedia of genes and genomes (KEGG) biological pathway (http://www.genome.jp/) analysis and a Gene Ontology enrichment analysis of gene function were used to identify the gene networks that may be modulated by the dysregulated miRNAs. Subsequently, miRNA transcription factors were identified through Chip-Seq next-generation genome sequencing and based on data included in the ENCODE and GEO databases [22,23]. Integrated miRNA regulatory networks were then constructed based on the above analysis.

Statistical analysis
The microarray data were normalized using a quantile normalization method. The p-values were calculated with the Rank Product Method [19,20], and differences with p< 0.01 were considered statistically significant. For analysis of the qRT-PCR results, a Wilcoxon rank-sum test was used, and differences with p< 0.05 was considered statistically significant [11]. Fisher's exact test (p< 0.05) was used to identify significant biological pathways, and the false discovery rates (FDR) were used as significant corrections [24]. A hypergeometric distribution was used for the Gene Ontology enrichment analysis of gene function, and p< 0.01 was considered statistically significant [25].

Result
Microarray assay of dysregulated miRNA expression in autism spectrum disorder The RNA extracted from the peripheral blood was analyzed using a RiboArray miDETECT Human Array. After normalization with the quantile normalization method, the hybridization data revealed 24 upregulated and 20 downregulated miRNAs (Table 1; more details in S1 Table). Interestingly, the dysregulated miR-15a and miR-15b previously discovered in autistic cerebellar cortices were found in this study to be also dysregulated in the peripheral blood [9]. miR-103, miR-92, and miR-19-2 have been reported to be aberrantly expressed in autistic lymphoblasts (Table 2), whereas the related miR-103a-3p, miR-92a-3p and miR-19b-3p were found to be abnormally expressed in autistic peripheral blood [8,10,26]. Thus, the differential expression of these miRNAs in the peripheral blood may partially reflect systemic changes in ASD.
A subset of differentially expressed miRNAs was selected for clustering analysis. The correlation of the expression profiles between the biological replicates and treatment conditions was demonstrated through unsupervised hierarchical clustering analysis (S1 File).

qRT-PCR confirmation of the expression of selected miRNAs
Based on the dysregulated miRNA profiles obtained from the microarray analysis combined with recent research on the function of miRNAs, the five most aberrantly expressed miRNAs were selected for further qRT-PCR validation. As a result, the miR-34b expression level was significantly increased, whereas the levels of miR-let-7a, miR-let-7d, miR-103a, and miR-1228 were decreased (Fig 1, p< 0.05).

Investigation of gene regulatory networks based on miRNA expression
It has been predicted that the regulation of miRNAs is not a simple process. More than one target in different functional categories can be modulated by one miRNA, and more than one miRNA may target the same specific gene [10,18]. Thus, the statistical data for the miRNA candidate target genes predicted by the four network databases (TargetScan, miRanda, CLIP-Seq, and miRDB) are presented in Fig 2. The Kyoto encyclopedia of genes and genomes (KEGG) biological pathway database was employed to extract the most relevant target genes on different biological pathways to predict the underlying connections among the target genes and their associated functions. The three most significant pathways for the five miRNAs were taken into account (Table 3). A Gene Ontology enrichment analysis of gene function was performed to fully explore the specific functions of the target genes. Interestingly, the results generated by the KEGG and GO analyses indicated correlations between the candidate targets of the miRNAs and autism and revealed other neurological functions previously and newly reported to be related to ASD, such as circadian rhythm-mammal [27], ubiquitin-mediated proteolysis (UPS) [28], dopaminergic synapse [29,30], MAPK signaling pathway, development of primary male sexual characteristics and male sex differentiation in male gonad development.
We also included miRNA transcription factors in our study. Therefore, based on the abovedescribed results, the S2 File shows the complex and integrated miRNA regulatory networks, which present the relationship between the miRNAs and mRNAs and the relationship between the miRNAs and miRNA transcription factors. Additionally, we calculated the number of associated genes for both candidate target genes and transcription factor genes in the miRNA regulatory interaction networks and ranked the number of associated genes, which was greater than 10 and is shown in

Discussion
Autism spectrum disorder (ASD) is one of the most heritable neurodevelopmental diseases with a markedly increasing incidence in China. The characteristics of ASD include repetitive behaviors accompanying deficits in social interaction and communication [31]. Because the diagnosis of ASD mainly relies on psychometric evaluations, its diagnosis is often difficult and subjective. Therefore, molecular evidence for early and objective ASD diagnosis particularly hsa-miR-15a down-regulated up-regulated cerebellar cortex Abu-Elneel et al. [9] hsa-miR-15b down-regulated up-regulated cerebellar cortex Abu-Elneel et al. [9] hsa-miR-92a-3p down-regulated down-regulated lymphoblastoid cell lines Talebizadeh et al. [10] doi:10.1371/journal.pone.0129052.t002 non-invasive miRNA/gene biomarkers from the peripheral blood and genetic associations or interactions, is indispensable. As described above, miRNAs are small non-coding regulatory RNAs that play an important role in fundamental neurobiological processes, particularly neurodevelopmental diseases. In Rett syndrome, multiple aberrantly expressed miRNAs, including miR-30a/d, miR-381, and miR-495, negatively regulate target gene translation [32]. In Down syndrome samples, five miRNAs (miR-155, miR-802, miR-125b-2, let-7c, and miR-99a) have been found to be overexpressed in fetal hippocampus and exclusively expressed in neurons [33]. In Fragile X syndrome, miR-367 negatively regulates FXR1P expression in human cell lines, and a loss of FXR1P downregulates brain-specific miRNAsthat are critical in neurodevelopment [34,35].
A series of recent studies have shown that miRNAs in serum and plasma are considered biomarkers for several diseases [36]. Additionally, several studies have identified the expression profiles of miRNAs in ASD, and some of these studies presented data from lymphoblastoid cell cultures [5,8,10]. One of these studies included the analysis of samples from the cerebella of ASD patients [9]. We explored the miRNA expression profiles and the gene regulatory networks of aberrantly expressed miRNAs in autism spectrum disorder in China. In this study, five miRNAs (miR-103a, miR-34b, miR-let-7a, miR-let-7d, and miR-1228) were confirmed to be aberrantly expressed in the peripheral blood of ASD patients. miR-34 modulates aging and neurodegeneration in Drosophila (dme-mir-34 and has-mir-34b are orthologous) [37], and miR-34b is elevated in spinocerebellar ataxia type 3/Machado-Joseph disease (SCA3/MJD) as demonstrated in our previous study [11]. These data illustrate the roles of miR-34b in both neurodevelopmental and neurodegenerative diseases. Mutiple results obtained in the current study indicate a relationship between miR-34b and ASD. The candidate target genes of miR-34b are mainly associated with central nervous system neuron development, long-term memory and tight junction pathways whose breakdown may result in diseases such as systemic inflammatory response syndrome (SIRS), inflammatory bowel disease, and autism [38]. Additionally, the strong association between the candidate target genes of miR-34b and ASD implies that miR-34b may play roles in the physiopathology of ASD. A linkage between autism and the EIF4E region on chromosome 4q has been found in genomewide linkage studies [39], and functional variants modulating EIF4E expression have been previously indicated as risk factors for ASD [40]. Dysregulation of eukaryotic translation initiation factor 4E (eIF4E)-dependent translational control may contribute to ASD-like phenotypes [41]. MAP2 encodes a protein that is enriched in dendrites and is implicated in determining and stabilizing dendritic shape during neuron development. Depleted MAP2 neuronal expression and reduced dendrite numbers have been found in autistic adults [42]. In addition, the MAP2 gene has been identified a good candidate for the generation of a preserved speech variant (PSV) in a patient with autism and Rett-like features [43]. Interestingly, the candidate target genes of miR-34b also affect the development of primary male sexual characteristics, male sex differentiation and male gonad development, which may explain the higher percentage of males in the ASD group. miR-103 and miR-107 have been reported to be highly expressed in brain tissue [44]. Sarachana et al. previously hypothesized that the dysregulation of miR-103/7 may result in abnormal lipid and fatty acid metabolism in ASD [8]. Several findings of our study provide plausible evidence that miR-103a-3p plays a role in the physiopathological mechanism of ASD. The candidate target genes of miR-103a-3p reveal close relationships with the dysfunctional physiological pathways related to ASD, such as central nervous system development, neuron projection morphogenesis and skeletal muscle tissue and organ development, and with various pathways, including ubiquitin-mediated proteolysis (UPS), circadian rhythm-mammal [8,26] and dopaminergic synapse [28,29]. Directly, HECT-type E3 ubiquitin ligase (UBE3A) dysfunction appears to induce a large spectrum of pediatric cognitive dysfunctions, such as autism and Rett and Angelman syndromes, depending on the mechanism and severity of UBE3A loss [27]. Indirectly, UPS is involved in neural functions and co-morbid diseases with ASD, such as synaptic dysfunctions and muscle diseases. Ubiquitination is a highly specific method for manipulating protein expression at the neural synapse [45,46] and overloading or other perturbations of ubiquitin-dependent proteolysis is likely involved in inclusion myopathies [27]. Based on the results of various previous studies [8], we hypothesized that the aberrantly expressed miR-103a-3p may contribute to abnormal ubiquitin-mediated proteolysis in ASD. Specifically, the candidate target gene BTRC is involved not only in the mammal circadian rhythm and mammal and ubiquitin-mediated proteolysis but also in the Wnt signaling pathway [47], all of which have been acknowledged to play roles in ASD. Additionally, one candidate target gene of miR-103a-3p, BDNF, encodes brain-derived neurotrophic factor (BDNF). Increasing evidence suggests a direct or indirect involvement of BDNF in autism [48], and this factor plays a key role in neuronal survival, morphology, differentiation and synaptic strength [49]. It is increased in some ASD patients as well as in some animal models of ASD [50,51].
Other candidate miRNAs have also associated with neuronal physiology and pathology. A previous study suggested that embryonic development may be relevant to ASD [8]. In addition, miR-let-7a expression have been found to be high in early human embryonic tissues and abruptly decreases thereafter, which suggests that miR-let-7a plays critical roles in early embryonic development [52]. Therefore, miR-let-7a may be indirectly associated with ASD through the embryonic development process. miR-let-7d suppresses neural stem cell proliferation by reducing TLX expression in neural stem cells, which may be a novel strategy for identifying potential interventions in relevant neurological diseases [53]. Notably, the candidate target genes of both miR-let-7a and miR-let-7d are associated with the MAPK signaling pathway, which directly or indirectly play roles in the physiopathology of ASD. Ghahramani Seno et al. previously suggested that mitogen-activated protein kinases (MAPK) may be involved in the development of autism [5]. In addition, several genes in the MAPK signaling pathway have been reported to be relevant to ASD. NGF encodes nerve growth factor (NGF) and influences inflammatory responses and neuron projection development. The serum NGF concentrations have been reported to be significantly higher in autistic children [54]. Another analysis suggests that NGF is a potential risk gene for nonverbal communication (NVC) impairments [55]. Recent studies have implied an association between mitochondrial dysfunction and ASD [56]. The protein levels of caspase-3, which is encoded by CASP3, are also increased in ASD patients [57]. NRAS plays roles in the regulation of long-term neuronal synaptic plasticity and transmission and is considered a candidate gene for autism. This finding may provide a novel insight into the etiology of autism [58].
Interestingly, the candidate target genes of miR-1228-3p were found to not be involved in brain-specific or brain-related functions but involved in the process of hepatocyte growth factor (HGF) production. A recent study reported that decreased HGF levels are found in autistic children and suggested that the serum HGF concentration may be a useful biomarker of autistic children [59]. Additionally, our findings show an association between the miR-1228-3p target genes and the insulin signaling pathway, which are hypothesized to contribute to the development of autism in genetically susceptible individuals by promoting PI3K/Tor pathway activation in neurons [60].
Taken together, four of the five aberrantly expressed miRNAs are related to the brain, suggesting that the gene expression differences observed in the peripheral blood may indicate similar alterations in the brain that are likely attribute to a system-scale differential expression of miRNAs.

Conclusion
In summary, we investigated the global miRNA expression profiles of ASD in China, confirming that the use of alternative and more easily accessible peripheral blood for the study of miRNA expression in ASD is valid. The differentially expressed miR-34b may potentially explain why ASD is markedly more common in males than in females. The aberrantly expressed miR-103a-3p may contribute to abnormal ubiquitin-mediated proteolysis in ASD. Thus, these five miRNAs are beneficial for better understanding ASD-specific miRNA expression profiles in peripheral blood and provide new insights into future researches for ASD.
Supporting Information S1 Table. More details of microarray analysis data. (XLSX) S1 File. Hierarchical cluster analysis. An unsupervised hierarchical cluster analysis of 150 significantly differentially expressed miRNAs among all of the autistic individuals (A6, A357, A14, A15, A16) and controls (4,567,28,29,30) shows the distinct miRNA expression pattern of the two groups (p < 0.05). (PDF) S2 File. Complex and integrated miRNA regulatory networks. The image shows not only the relationship between miRNAs and mRNAs or transcriptional factors (TFs) but also the interactive relationships between the proteins of the target genes. Orange rhombi: miRNAs, blue circles: mRNAs, aubergine circles: TFs, pink lines: miRNA-mRNA interactions, blue lines: mRNA-TF interactions, green lines: miRNA-TF interactions. (PDF)