Expression and Regulatory Network Analysis of miR-140-3p, a New Potential Serum Biomarker for Autism Spectrum Disorder

Given its prevalence and social impact, Autism Spectrum Disorder (ASD) is drawing much interest. Molecular basis of ASD is heterogeneous and only partially known. Many factors, including disorders comorbid with ASD, like TS (Tourette Syndrome), complicate ASD behavior-based diagnosis and make it vulnerable to bias. To further investigate ASD etiology and to identify potential biomarkers to support its precise diagnosis, we used TaqMan Low Density Array technology to profile serum miRNAs from ASD, TS, and TS+ASD patients, and unaffected controls (NCs). Through validation assays in 30 ASD, 24 TS, and 25 TS+ASD patients and 25 NCs, we demonstrated that miR-140-3p is upregulated in ASD vs.: NC, TS, and TS+ASD (Tukey's test, p-values = 0.03, = 0.01, < 0.0001, respectively). ΔCt values for miR-140-3p and YGTSS (Yale Global Tic Severity Scale) scores are positively correlated (Spearman r = 0.33; Benjamini-Hochberg p = 0.008) and show a linear relationship (p = 0.002). Network functional analysis showed that nodes controlled by miR-140-3p, especially CD38 and NRIP1 which are its validated targets, are involved in processes convergingly dysregulated in ASD, such as synaptic plasticity, immune response, and chromatin binding. Biomarker analysis proved that serum miR-140-3p can discriminate among: (1) ASD and NC (Area under the ROC curve, AUC: 0.70; sensitivity: 63.33%; specificity: 68%); (2) ASD and TS (AUC: 0.72; sensitivity: 66.66%; specificity: 70.83%); (3) ASD and TS+ASD (AUC: 0.78; sensitivity: 73.33%; specificity: 76%). Characterization of miR-140-3p network would contribute to further clarify ASD etiology. Serum miR-140-3p could represent a potential non-invasive biomarker for ASD, easy to test through liquid biopsy.


INTRODUCTION
Autism Spectrum Disorder (ASD) is the name for a heterogeneous group of complex neurodevelopmental conditions, which are clinically defined by: (1) defects in social interaction and communication; (2) fixed interests and repetitive behaviors. Typically, ASD symptoms become fully manifest during school age and have a lifelong impact on everyday functions (American Psychiatric Association, 2013).
The broadening of the autism concept and the resulting changes in ASD categorization have increased ASD awareness and improved its diagnostic surveillance in health and educational care. This has led to an alarming rise in the number of milder cases of ASD, without co-occurring intellectual disability, in developed countries around the world.
Recently, it has been reported that ASD affects one in 68 US children and that approximately four males suffer from ASD for every female (Christensen et al., 2016). Comorbid neuropsychiatric and neurodevelopmental disorders contribute to ASD impairment, being common (70.8%) and frequently multiple (57%) in ASD children (Simonoff et al., 2008). Such conditions include social anxiety disorder, attentiondeficit/hyperactivity disorder (ADHD), oppositional defiant disorder (ODD), chronic tic disorder, and obsessive-compulsive disorder (OCD).
Tourette Syndrome (TS) is a neurodevelopmental disorder characterized by considerable motor as well as behavioral impairment: it affects approximately 1% of the population with a male:female ratio of 3:1. It is clinically defined by childhood onset of multiple motor tics and at least one phonic tic, which collectively must persist for at least 12 months (American Psychiatric Association, 2013). 88% of TS patients also show comorbidity and psychopathology. Comorbidity with ADHD and OCD is very common. Co-existent psychopathologies include depression, anxiety, learning difficulties, personality disorder, ODD, and conduct disorder (Robertson, 2012). TS can cause severe difficulties in social functioning and a reduced quality of life in patients suffering from it. Characteristic, but not essential for diagnosis, symptoms include complex tics, such as echolalia and echopraxia (copying others' vocalizations and actions, respectively), palilalia and palipraxia (repeating own words/phrases and actions, respectively), coprolalia (inappropriate involuntary swearing) as well as repeating of complex words (Robertson, 2015).
It has been observed that 4.8% of ASD children also suffer from TS (Simonoff et al., 2008) and that 6-11% of TS cases show comorbidity with ASD (Robertson, 2012). TS and ASD share clinical symptomatology and many behavioral features. Genetic studies also support the existence of common susceptibility genes in both disorders (Clarke et al., 2012). The exact etiology of both disorders is still elusive.
Strong evidence suggests that ASD may arise from genetic, epigenetic and environmental factors (Abdolmaleky et al., 2015;Nardone and Elliott, 2016;Sun et al., 2016). ASD is genetically highly heterogeneous. Both inherited and de novo ASD-associated variants have been characterized in hundreds of genes. Both inherited and de novo rare genetic variants can be detected in 10-30% of ASD cases. Single common inherited variants can be found in approximately 1.1-1.2% of ASD cases; when considered cumulatively, these can explain 15-50% ASD cases. However, no common risk loci have been identified to date for ASD (Vorstman et al., 2017).
Resultant micro-and macro-structural and functional abnormalities, which emerge during brain development in ASD, create the dysfunction of neural networks involved in socio-emotional processing (Vissers et al., 2012;Maximo et al., 2014;Kern et al., 2015;Ecker, 2017).
Since common risk loci for ASD have not been proposed yet (Vorstman et al., 2017), a molecular test for non-syndromic ASD is not available and diagnosis relies on clinical assessment and confirmation. Clinical diagnosis of ASD depends on behavioral observations, according to the Diagnostic and Statistical Manual of Mental Disorders (5th ed.; DSM-5;American Psychiatric Association, 2013). Accepted gold standard tools for diagnostic assessment of ASD are the Autism Diagnostic Observation Schedule (ADOS) and the Autism Diagnostic Interview-Revised (ADI-R) (Falkmer et al., 2013). Considering the clinical variation and etiological heterogeneity of ASD, a precise diagnosis can be very difficult. It is further complicated by ASD association with comorbid disorders (Constantino and Charman, 2016). Therefore, there is an urgent need for potential ASD biomarkers that could support clinical discrimination of patients.
MicroRNAs (miRNAs) are conserved post-transcriptional regulators of gene expression, collectively increasing the precision and robustness of gene-regulatory networks and affecting all cellular pathways, from development to metabolism (Berezikov, 2011). Other than in neoplastic and degenerative diseases, their dysregulation plays a role in several neurodevelopmental and neuropsychiatric disorders (Geaghan and Cairns, 2015).
Since their identification and characterization in serum and plasma of humans and other animals (Chen et al., 2008;Chim et al., 2008), extracellular miRNAs have attracted researchers for their potential as new non-invasive tools for diagnosis, prognosis, and treatment evaluation of many human diseases. Extracellular miRNAs can be detected in all mammalian body fluids. Stability and general consistency of levels among individuals, along with the existence of specific expression signatures in association with both physiological and pathological conditions, make circulating miRNAs appropriate biomarkers; these also suggest the prospective alternative use of liquid biopsies as sources of biomarkers in the clinic (Weiland et al., 2012;Larrea et al., 2016). Numerous studies have identified either serum or plasma circulating miRNAs as promising biomarkers for neurodevelopmental and neuropsychiatric disorders, as ADHD (Wu et al., 2015), TS (Rizzo et al., 2015), depression/anxiety disorder (Wang X. et al., 2015), posttraumatic stress disorder (Balakathiresan et al., 2014), bipolar disorder (Rong et al., 2011), schizophrenia , and epilepsy (Wang et al., 2015a,b;An et al., 2016), and other brain pathological conditions, like traumatic brain injury (Di Pietro et al., 2017) and vascular dementia (Ragusa et al., 2016).
Although ASD research is progressively and actively growing, only two papers have characterized circulating miRNAs in body fluids from ASD patients with a high-throughput approach (Mundalil Vasu et al., 2014;Hicks et al., 2016); none of these studies has also focused on patients affected by other neurodevelopmental disorders comorbid with autism.
We hypothesized that the serum profile of circulating miRNAs may contain some specific fingerprints for ASD which could also support in the discrimination among it and comorbid neurodevelopmental disorders. Aiming to gain more knowledge about ASD biomolecular basis and identify new potential biomarkers for this disorder, we exploited a high-throughput approach to analyze the expression of circulating miRNAs in serum from ASD, TS, and TS+ASD patients.
Following profiling of serum miRNAs through our previously published protocol (Rizzo et al., 2015;Ragusa et al., 2016), we validated serum miR-140-3p as significantly upregulated in ASD patients compared to unaffected controls (NCs), TS patients, and TS+ASD patients. Also, we demonstrated that its serum expression levels are correlated with scores from the tic scale YGTSS (Yale Global Tic Severity Scale). Then, we observed that miR-140-3p network node genes are involved in biological processes convergingly dysregulated in ASD (i.e., synaptic plasticity, immune response, and chromatin binding). Finally, through biomarker analysis, we proved that serum miR-140-3p can discriminate ASD from NC, TS, and TS+ASD.

Patient Selection
From a database of more than 2,000 patients (from the Section of Child and Adolescent Psychiatry, Department of Clinical and Experimental Medicine, University of Catania), 79 Caucasian patients, aged 3-13 years from various socio-economic contexts, were randomly recruited and studied from January to November 2016. Thirty patients affected by ASD [mean age 6.5 (standard deviation, SD 3.5); M:F 22:8], 24 patients affected by TS [mean age 8.7 (SD 5.2); M:F 21:3], and 25 patients affected by TS+ASD [(mean age 9.3 (SD 6.7); M:F 25:0] were included in the study. They were compared to 25 neurologically intact unaffected negative control (NC) individuals [mean age 9.5 (SD 3.9); M:F 16:9], recruited from local schools, without any history of either ASD or TS and who suffered from neither chronic diseases nor psychiatric disorders ( Table 1).
The study was approved by the local Ethics Committee. All parents gave written informed consent.
Diagnoses of ASD, TS and other clinical conditions were made according to both DSM-IV-TR (Diagnostic and Statistical Manual of Mental Disorders, IV edition-Text Revision) and DSM-5 criteria by a child neurologist (RR). All the participants were evaluated at the University Hospital Policlinico -Vittorio Emanuele of Catania. The three clinical groups (ASD, TS, and TS+ASD) and the NCs were assessed using the following scales/schedules: ADOS and ADI-R to evaluate ASD symptoms; YGTSS to evaluate presence and severity of tics. Moreover, the three clinical groups (ASD, TS, and TS+ASD) and the NCs were also assessed by a psychologist through WISC-III (Wechsler Intelligence Scale for Children, III edition) as an evaluation of both IQ (Intelligence Quotient) and cognitive functioning. Neuropsychological features of participants are summarized in Table 1.
In a previous study (Rizzo et al., 2015), we reported that only miR-429 was significantly differentially expressed (DE) in the serum of TS patients compared to NCs: TS patients were included in this experimental series aiming to compare them with the other classes of neuropsychiatric patients.

Sample Processing
Peripheral blood samples from all participants were taken in the morning using a butterfly device into serum separator collection tubes, provided with Clot activator and gel for serum separation as additives (BD Biosciences). Collection tubes were treated according to current procedures for clinical samples. In order to separate serum from blood cells, tubes were rotated end-over-end at 20 • C for 30 ′ and then centrifuged at 3,500 rpm at 4 • C for 15 ′ in a Beckman J2-21. Supernatants were aliquoted into 1.5 ml RNasefree tubes and stored at −80 • C. Prior to RNA extraction, stored supernatants were centrifuged again at 3500 rpm at 4 • C for 15 ′ to remove circulating cells or debris. Serum samples were aliquoted into 1.5 ml RNase-free tubes and they were either immediately used for RNA extraction or stored at −80 • C until analysis (Rizzo et al., 2015).

RNA Extraction
RNA was extracted from 400 µl serum samples using Qiagen miRNeasy Mini Kit (Qiagen, GmbH, Hilden, Germany), according to Qiagen Supplementary Protocol for purification of total RNA, including small RNAs, from serum or plasma. RNA was eluted in a 40 µl total volume of RNase-free water with two consecutive steps of elution (30 µl followed by another 10 µl of RNase-free water) performed in the same collection tube.

MiRNA Profiling
We  (Ragusa et al., 2016). We individually carried out the analysis on microfluidic cards A and B. We used a customized normalization approach for the relative quantification analysis. Supplementary File 1 reports Ct values for TLDA panels A and B and shows our procedure step by step. For each comparison, a Ct value matrix (miRNAs in rows, samples in columns) was created. In a similar way to the GMN (global median normalization) method (Park et al., 2003), for each sample of the comparison, the median and mean Ct values within the array, reflecting the loaded mass of template cDNA, were calculated. However, all Ct values representing a specific miRNA were kept out of these calculations if even just one of them corresponded to a flagged well. Then, using the Pearson correlation, miRNAs whose expression profile was closer (more positively correlated) to these values were identified as the best endogenous controls within the arrays. We normalized MiRNAs to the top three stable miRNAs within the arrays. MiR-146a and miR-223 * were the most frequently stable miRNAs for cards A and B, respectively, and the most abundant among those we could select.
Therefore, for each comparison, three Ct value matrices (miRNAs in rows, samples in columns) were produced according to the 2 − Ct method (Schmittgen and Livak, 2008). DE circulating miRNAs were obtained performing SAM (Significance of Microarrays Analysis) statistical analysis on these matrices with MeV (Multi experiment viewer v4.8.1) statistical analysis software 1 . For each pairwise comparison, we used a two-class unpaired test, based on at least 100 permutations per miRNA, with a FDR (False Discovery Rate) cut-off of 0.15, in order to detect dysregulated miRNAs. This analysis identified many DE miRNAs for each comparison. However, we have used very strict criteria to select miRNAs for further validation (i.e., number of SAM tests in which they were identified as DE, number of comparisons in which they resulted as DE, their abundance and the quality of their amplification curves during the profiling runs) in order to investigate only the most promising ones.

MiRNA Profiling Data Validation
RNA from sera of 30 ASD, 24 TS, and 25 TS+ASD patients and 25 NCs was used to perform miRNA-specific reverse transcription reactions producing miRNA-specific cDNAs for real-time PCRs. These RT-PCR analyses were performed using TaqMan MicroRNA Assays (Applied Biosystems) specific for the most interesting miRNAs identified as DE, miR-30d, miR-140-3p, miR-148a * , and miR-222, and for the selected endogenous control, miR-146a. At first, the ASD group was composed of 32 patients. We checked if some of those samples should be considered as outliers, within this original ASD group, for: (1) the serum expression of miR-140-3p; (2) the severity of ASD symptoms. We have looked at their Ct values for miR-140-3p and at scores that they obtained for the four items of the ADOS scale (A: Communication; B: Social interaction; C: Imagination; D: Repetitive and restricted behaviors). For these expression values and ADOS scores, we defined the corresponding mean ± 2 * (SD) ranges and we considered patients with a value and/or score outside of those ranges as outliers. Two ASD patients were excluded from the original ASD group since: (1) both were outliers for miR-140-3p expression; (2) one was an outlier for the Imagination item (1/4 ADOS items), whereas the other one was an outlier for the Imagination and Repetitive and restricted behaviors items (2/4 ADOS items). All the following analyses were performed with GraphPad Prism for Windows v6.01 (GraphPad Software, La Jolla California USA 2 ). D'Agostino-Pearson omnibus K2 test and Shapiro-Wilk test were performed to check if data from every small group were normally distributed. Ordinary one-way ANOVA was used to test the differential expression of the selected miRNAs between the four groups. Statistical significance was established at a p ≤ 0.05. Tukey's multiple comparisons test was performed to identify which groups differed in the selected miRNAs' expression. Statistical significance was established at a multiplicity adjusted p ≤ 0.05. Expression FC (Fold Change) values of miRNAs were calculated by applying the 2 − Ct method (Schmittgen and Livak, 2008).

Correlation between miR-140-3p Expression and Neuropsychiatric Parameters
Correlation between Ct values for miR-140-3p, obtained from the normalization to miR-146a, and neuropsychiatric parameters was analyzed in both a general (all patients and controls) and a class-specific (just one class of patients and controls) way, since some of these parameters were related to a certain class of neuropsychiatric disorders. IQ, ADOS items regarding communication, social interaction, imagination, and repetitive and restricted behaviors (ADOS items A-D), and YGTSS (Yale Global Tic Severity Scale) were the neuropsychiatric parameters chosen for this analysis. Either Pearson or Spearman correlation was computed on GraphPad Prism software when analyzing normally and not normally distributed data, respectively. Two-sided p-values from this correlation analysis were corrected for multiple comparisons by using three different approaches: Bonferroni correction, Holm-Bonferroni correction, and Benjamini-Hochberg (BH) FDR procedure. Statistical significance was established at a p-value ≤ Bonferroni corrected α = 0.05/16 = 0.003, at a Holm-Bonferroni corrected p ≤ 0.05, and at a Benjamini-Hochberg FDR adjusted p ≤ 0.01. Linear regression analysis was also carried out on GraphPad Prism software only for significant correlations. Statistical significance was established at a p ≤ 0.05.

Computational Analyses
Reconstruction of the miR-140-3p-mediated Regulatory Network MiR-140-3p targets whose validation was based on strong evidence were retrieved by DIANA-TarBase v7.0 3 (Vlachos et al., 2015) and miRTarBase 4 (Chou et al., 2016) databases. The biological network, composed of MIR140, these targets, and their first neighbors, was built retrieving interactome data through BisoGenet v3.0.0 Plug-in (Martin et al., 2010) in Cytoscape v3.4.0 (Shannon et al., 2003). Network centralities analysis, permitting the identification of the nodes that, more than others, were good candidates as regulators of the underlying biological processes in which the network is involved, was carried out through CentiScaPe v2.1 Plug-in (Scardoni et al., 2014).
Network Functional Analysis clusterProfiler v3.2.11 R package (Yu et al., 2012) was used to perform functional enrichment analyses on miR-140-3pmediated regulatory network node genes in R v3.3.2 5 (R Core Team, 2016). We searched for the gene annotation terms from the GO (Gene Ontology), DO (Disease Ontology), KEGG (Kyoto Encyclopedia of Genes and Genomes), and Reactome databases that were over-represented in the list of network node genes compared to the entire genome. Statistical significance for the hypergeometric test was established at a BH adjusted p ≤ 0.05. gofilter() and simplify() functions in clusterProfiler were employed in order to select level-specific GO terms and to remove the most redundant ones, respectively.

Network Expression Analysis
In order to investigate if deregulation of network node genes was implicated in ASD, we searched for raw high-throughput gene expression datasets produced from the analysis of samples of ASD patients on two public repositories, GEO (Gene Expression Omnibus) DataSets (Edgar et al., 2002) and ArrayExpress (Kolesnikov et al., 2015). Datasets retrieved by GEO DataSets were analyzed performing limma tests with the GEO2R tool 6 ( Barrett et al., 2013), whereas datasets retrieved by ArrayExpress were analyzed performing Tusher SAM tests with MeV software. We reported only network node genes whose log 2 FC expression was significantly higher than 1 and lower than -1 as upregulated and downregulated, respectively, within ASD datasets. MeV software was also used to produce the curated ASD expression heatmap. Supplementary File 2 reports all the datasets selected for network functional analysis along with their references and IDs (Gregg et al., 2008;Kuwano et al., 2011;Voineagu et al., 2011;Ginsberg et al., 2012;Kong et al., 2012; see Supplementary File 2).

ROC Curve Analysis and Biomarker Performance Evaluation
Ct values for miR-140-3p, obtained from the normalization to miR-146a, served as input data to perform a classical univariate ROC (Receiver Operator Characteristic) curve analysis for each of the comparisons where we found this miRNA to be dysregulated on the server Metaboanalyst 3.0 7 (Xia and Wishart, 2016). An appropriate Ct cut-off point maximizing both sensitivity and specificity (that is, the threshold that maximizes the distance to the diagonal line) was found for each curve by calculating the maximum Youden index J (max [(sensitivity + specificity)-1]). GraphPad Prism software was used to create Figure 4. The true positive rate (y-axis) was plotted in function of the false positive rate (x-axis), for different Ct cut-off points.
Since these ROC curves were based on a miRNA already identified as differentially expressed between the compared groups (miR-140-3p), through them we could only assess its idealized discriminative power. It is possible that this miRNA only accurately predicts outcomes in the initial data set and that minor fluctuations in the training data could markedly lower its predictive performance.
Therefore, after these preliminary ROC curve analyses, we built corresponding logistic regression models for the expression of miR-140-3p and we tested them through CV (cross-validation) and permutation testing, once again, by using the server Metaboanalyst 3.0. CV gives an indication of how accurate a given model might be in predicting new samples, validating its general applicability (Xia et al., 2013). 100-time repeated random sub-sampling CV was used to test the performance of the built logistic regression models. At each CV, 2/3 of samples are used for model training and the remaining 1/3 of samples are used for testing. Permutation testing indicates if a given model is significantly different from a random guessing model for the sample population, validating the proposed model structure (Xia et al., 2013). Permutation testing on the performance measure AUC (Area under the ROC curve) was used to calculate the significance of the built logistic regression models. The permutation tests use this procedure: random label re-assignment to each sample; 3-time repeated random subsampling CV; comparison of the performance measures between the models obtained by using the original and the permuted sample labels. This procedure was repeated 100 times. If the performance measure of the original data lies outside the normal distribution of the one of the permuted data, then the tested model is significant. Statistical significance was established at a p ≤ 0.05.

High-Throughput Expression Analysis of Circulating miRNAs in ASD, TS, and TS+ASD Patients
By using TLDA technology, we investigated the expression levels of 754 miRNAs in sera from four ASD patients, five TS patients, four TS+ASD patients, and three unaffected NCs. We identified miR-146a and miR-223 * as the best endogenous controls for panels A and B, respectively. Supplementary File 1 reports Ct values for TLDA panels A and B.

Dysregulated Expression Levels of miR-140-3p in Serum from ASD Patients
We selected miR-30d, miR-140-3p, miR-148a * , and miR-222 for further validation through single TaqMan assays. MiR-146a was used as endogenous control in all the analyses carried out.

Correlation between miR-140-3p Expression and Neuropsychiatric Parameters
In order to test if any link existed between serum expression of miR-140-3p and commonly used ASD and TS neuropsychiatric parameters, we computed the correlation between Ct values for miR-140-3p and various neuropsychiatric scores of study participants ( Table 2). Supplementary File 4 reports Ct value for miR-140-3p and neuropsychiatric scores of each participant.
When we included all patients and controls in our analysis, we found a positive correlation (Spearman r = 0.33; twosided p = 0.0005, significant according to Bonferroni correction; Holm-Bonferroni corrected p = 0.008; Benjamini-Hochberg FDR adjusted p = 0.008) and a linear relationship (y = 4.537x-3.269, y: YGTSS score, x: Ct value for miR-140-3p, two-sided p = 0.002) between miR-140-3p expression levels and YGTSS scores (Figure 2). We could assume that lower serum levels of miR-140-3p, which have been observed in TS+ASD patients, are potentially associated with both occurrence and worsening    of motor and phonic tics and that higher serum levels of miR-140-3p, which we found in ASD patients, could be linked to the absence of tics. This analysis confirmed that serum expression of miR-140-3p correlated with a crucial neuropsychiatric scale for the clinical diagnosis of TS. We infer that miR-140-3p could prove to be useful to strengthen the behavior-based diagnosis of either ASD or TS+ASD, which can be particularly challenging in some clinical cases.

Reconstruction of miR-140-3p-mediated Regulatory Network: Functional and Expression Analyses of Network Node Genes
By searching on online databases of miRNA-mRNA interactions for validated targets of miR-140-3p, we retrieved CD38 (CD38 molecule) and NRIP1 (nuclear receptor interacting protein 1) as its only targets whose validation was based on strong evidence. Through network analysis, we reconstructed the regulatory network composed of MIR140 (microRNA 140, the gene encoding miR-140-3p), CD38, NRIP1, and their first neighbors. This network had 111 nodes and 821 edges. NRIP1, POLR2A (RNA polymerase II subunit A), EP300 (E1A binding protein p300), E2F1 (E2F transcription factor 1), ESR1 (estrogen receptor 1), PHF8 (PHD finger protein 8), and TAF1 (TATA-box binding protein associated factor 1) were the nodes with the highest degree within it (Supplementary Figure 1).
In order to investigate the potential etiological role of this miRNA-mediated network in ASD, we performed functional enrichment analysis of network node genes using GO, DO, KEGG, and Reactome gene annotation databases (Supplementary Figures 2-4). Supplementary File 5 reports all the data from network functional analysis.
Genes from miR-140-3p-mediated regulatory network played a role in various mechanisms within the nervous system (i.e., neurogenesis, regulation of synaptic plasticity, long term synaptic depression, cellular response to nerve growth factor, neuron differentiation, dendrite development, and neuronal death). In addition to their role in nervous system development, they were also involved in growth regulation, endocrine system development, heart development, respiratory system development, and tongue development (Table 3, Supplementary Figure 5). MIR140 gene was annotated with the over-represented DO term physical disorder, that refers to diseases determined by a genetic abnormality, error  Term database ID, term description, corresponding BH adjusted p-value generated by the hypergeometric test, Gene Ratio, and Background Ratio values are reported for all the GO Biological Process terms regarding nervous system and development. Gene Ratio: the ratio of number of genes of interest that are annotated with a certain term from the database used to perform the analysis to number of genes of interest that are annotated with terms from the same database. Background Ratio: the ratio of number of genes in the genome that are annotated with a certain term from the database used to perform the analysis to number of genome genes that are annotated with terms from the same database. GO, Gene Ontology; BP, Biological Process; BH, Benjamini-Hochberg.
with embryonic development, infection or compromised intrauterine environment (DOID:0080015; BH adjusted p = 0.027; Gene Ratio = 0.070; Background Ratio = 0.017). Among the most interesting terms whose enrichment was determined by CD38, we found those regarding response to estradiol, retinoic acid, drugs, hypoxia, ketone, and oxidative stress, activation and proliferation of immune cells, regulation of protein localization, cellular calcium ion homeostasis, and blood circulation. We found bacterial infectious disease as the only DO term among them. Finally, CD38 was directly involved in the regulation of synaptic plasticity after and in long-term synaptic depression (Table 4, Supplementary Figure 6). Among the most interesting terms whose enrichment was determined by NRIP1, there were those related to response to estradiol and steroid hormones, reproductive system development, and development of primary sexual characteristics. NRIP1 was annotated with many molecular functions (i.e., histone deacetylase binding, nuclear hormone receptor binding, core promoter sequence-specific DNA binding, retinoic acid receptor binding, and retinoid X receptor binding). Finally, NRIP1 regulated the transcription of genes involved in circadian rhythms by interacting with RORA (RAR related orphan receptor A) (Table 5, Supplementary Figure 7).
To verify if dysregulation of network node genes was implicated in ASD, we used publicly available raw highthroughput gene expression datasets produced from the analysis  of ASD samples. Ten genes were found to be downregulated in whole blood of ASD patients compared to NCs (log 2 FC <-1): CCDC85B (coiled-coil domain containing 85B), CD3E (CD3e molecule), CIB1 (calcium and integrin binding 1), CTBP1 (C-terminal binding protein 1), LCK (LCK proto-oncogene, Src family tyrosine kinase), MAP3K7 (mitogen-activated protein kinase kinase kinase 7), NR1H2 (nuclear receptor subfamily 1 group H member 2), RARA (retinoic acid receptor alpha), STAT3 (signal transducer and activator of transcription 3), and ZAP70 (zeta chain of T cell receptor associated protein kinase 70); two were found to be overexpressed (log 2 FC > 1): PHF8 and TAF1 (Figure 3). Supplementary File 2 reports all the data from network expression analysis.

Serum Levels of miR-140-3p in the Discrimination of ASD Patients
We used Ct values for miR-140-3p to perform a classical univariate ROC curve analysis for each of the comparisons where we found this miRNA to be dysregulated. The univariate ROC plots revealed an AUC of 0.71 for the comparison ASD vs. NC (p = 0.006), 0.73 for ASD vs. TS (p = 0.002), and 0.78 for ASD vs. TS+ASD (p = 0.00007) (Figure 4). We used Ct value cut-offs corresponding to the sensitivity/specificity pair with the highest Youden index J for every computed ROC curve to perform a blind diagnosis on all the 104 analyzed samples (Figure 5). Then, we built a logistic regression model for miR-140-3p expression in each comparison and we tested those predictive models through CV and permutation testing. 100-time repeated random sub-sampling CV was used to test the performance of the logistic regression models. MiR-140-3p continued to perform at a good level for the comparison ASD vs. NC, with an average AUC of 0.70, a sensitivity of 63.33%, and a specificity of 68% (Figures 6 1C,1D). MiR-140-3p continued to perform at a good level also for the comparison ASD vs. TS, with an average AUC of 0.72, a sensitivity of 66.66%, and a specificity of 70.83% (Figures 6 2C,2D). MiR-140-3p continued to perform at a very high level for the comparison ASD vs. TS+ASD, with an average AUC of 0.78, a sensitivity of 73.33%, and a specificity of 76% (Figures 6 3C,3D). CV results demonstrated the general applicability of these predictive models. 100-time repeated permutation tests on the performance measure AUC were carried out to validate the structure of these models. Permutation testing results were significant and quite stable in different runs for all the models tested (Figures 6 1E, 2E, and 3E).
These data proved that serum miR-140-3p could be used in the discrimination of ASD patients. In particular, it could potentially support the differential behavior-based diagnostic process of two classes of neurodevelopmental disorders, ASD and TS+ASD.

DISCUSSION
MiRNAs are conserved post-transcriptional regulators of gene expression, coordinating development, homeostasis, and response to external factors in organisms (Berezikov, 2011). MiRNAs are main players also in the brain, where they control many neurodevelopmental processes, including patterning, cell specification, local translational control of neuronal plasticity, neurogenesis, and neuronal apoptosis (Kosik, 2006). Some miRNAs show differential spatio-temporal and sex-biased expression patterns in the developing human brain and regulate targets that are highly enriched for genes related to transcriptional regulation, neurodevelopmental processes, and common neurodevelopmental disorders, such as ASD, schizophrenia, and bipolar disorder (Ziats and Rennert, 2014). Since miRNAs are critical for various developmental mechanisms and pathologic transcriptional processes, many studies have investigated their role in neurodevelopmental and psychiatric disorders (Omran et al., 2012;Geaghan and Cairns, 2015;Scott et al., 2015), providing some clues to the etiology of these complex disorders.
Many factors, including disorders comorbid with ASD, like TS, complicate ASD behavior-based diagnosis and make it vulnerable to bias. Stability, general consistency of expression among individuals, and condition-specific expression profile make circulating miRNAs appropriate non-invasive diagnostic biomarkers (Weiland et al., 2012;Larrea et al., 2016). Expression profiling of miRNAs circulating in either serum or plasma has already been used as a biomarker discovery approach for several psychiatric disorders (Rong et al., 2011;Balakathiresan et al., 2014;Rizzo et al., 2015;Wang et al., 2015a,b;Wang X. et al., 2015;Wei et al., 2015;Wu et al., 2015;An et al., 2016;Ragusa et al., 2016;Di Pietro et al., 2017).
We hypothesized that serum profile of circulating miRNAs may contain specific fingerprints for ASD: these could provide some hints to the molecular basis of ASD and also be used as supportive means to the clinical diagnostic process, especially in the discrimination among neurodevelopmental disorders.
FIGURE 5 | The potential use of serum miR-140-3p as a biomarker: criteria for ASD diagnosis. The graphs show the distribution of Ct values of all the 104 analyzed samples, for which we already had a clinical diagnosis. We used data from classical univariate ROC curve analyses to perform a blind diagnosis of all study participants. In (A), the Ct ≤ 2.427 criterion divides ASD patients from NCs and determines the correct discrimination of 19/32 ASD patients and 20/25 NCs. In (B), the Ct ≤ 2.447 criterion divides ASD patients from TS patients and determines the correct discrimination of 19/32 ASD patients and 19/24 TS patients. In (C), the Ct ≤ 2.824 criterion separates ASD patients from TS+ASD patients and determines the correct discrimination of 22/32 ASD patients and 19/25 TS+ASD patients.
Using TLDA technology, we profiled serum expression of 754 human miRNAs in a discovery set of samples, including four ASD, five TS, and four TS+ASD patients and three NCs. The undeniable limit of this profiling approach is that, despite analyzing a total number of 16 samples in this exploratory step, we could have missed some significantly dysregulated circulating miRNAs because of the small number of samples in each compared group. However, we also applied very strict selection criteria for the identification of those 9 miRNAs (miR-30d, miR-30e-3p, miR-140-3p, miR-148a * , miR-222, miR-454, miR-483-5p, miR-1274B, and miR-1290) as differentially expressed, in order to spot more robust findings, likely to be confirmed in the next validation step. Nevertheless, through single TaqMan assays we observed the dysregulation of just one (miR-140-3p) of the 4 miRNAs selected for validation (miR-30d, miR-140-3p, miR-148a * , and miR-222) even if, according to results from profiling, those showed the most interesting and marked expression trends. This could be explained by the fact that TLDA approach is still a highthroughput approach: the preliminary step of preamplification suggested for TLDA analysis may have inserted an amplification bias, leading to not very accurate results. Knowing how crucial validation is, we have also performed this step by TaqMan assay, a probe-based system that is designed to specifically detect the expression of miRNAs of interests.
There is no consensus about optimal normalization strategies and accurate reference genes for both intracellular and extracellular miRNAs (Schwarzenbach et al., 2015). Suggested endogenous controls for miRNAs differ depending on: (1) the species considered; (2) either the tissue or body fluid analyzed; (3) either the physiological or pathological condition investigated; (4) different sample preparation methods, especially for circulating miRNAs (Marabita et al., 2016). On one hand, researchers can select reference miRNA genes according to reports from similar studies in the literature, and then, validate them for the set of samples under analysis (Schwarzenbach et al., 2015). On the other hand, they can screen the specific samples for more suitable endogenous controls through the profiling of a large number of genes. In this case, it is suggested to use either the mean expression values of whole miRNAs in each of the samples or genes with expression levels similar to these values as the screened endogenous controls (Schwarzenbach et al., 2015;Marabita et al., 2016). At the moment, accurately described normalization approaches and validated reference genes for serum miRNAs in ASD patients lack: this, together with different ethnicity of participants and other potential diverse analytic variables in studies similar to ours, led us to prefer our customized normalization approach over reference gene selection from literature. This approach, which is described in Supplementary File 1 and is inspired to the recommended miRNA array normalization strategy reported above (Schwarzenbach et al., 2015;Marabita et al., 2016), has proved its value in other works on neurological disorders published by our group (Rizzo et al., 2015;Ragusa et al., 2016). Thanks to it, we have selected miR-146a and miR-223 * as the most appropriate and accurate reference genes for our system. Some studies have identified miR-146a as either dysregulated in human ASD tissues (Mor et al., 2015;Nguyen et al., 2016) or associated to inflammation and immune response observed in neurodegenerative and neurological disorders (Iyer et al., 2012;Kiko et al., 2014;Müller et al., 2014;Wang et al., 2015b;An et al., 2016;Romano et al., 2017): however, this neither discourages its use as reference miRNA in our dataset nor affects the applicability of our results from differential expression analysis. Probability scores more than 0.5 belong to the ASD group, those less than 0.5 belong to the other group. Incorrectly classified subjects are identified by their ID number. (E) Results from the permutation tests on the model performance measure AUC. Average ROC curve and corresponding p-value are reported. AUC, Area under the ROC curve; 95% CI: 95% Confidence Interval; CV, cross-validation.
Through our expression analysis, we have identified miR-140-3p as dysregulated in serum from ASD patients. It is upregulated in ASD patients compared to NCs, TS patients, and TS+ASD patients: its levels are the highest in the ASD group (Figure 1). We observed that miR-140-3p levels are the lowest in the TS+ASD group. It is interesting that the two groups ASD and TS+ASD showed such a different expression trend for this miRNA. It is likely that the lower miR-140-3p expression reflects the presence of phonic and motor tics due to the comorbidity of TS with ASD. It may also depend on other either physiological conditions or comorbidities: this needs to be further investigated. On the contrary, according to our previous study on TS (Rizzo et al., 2015), we expected not to find any difference in miR-140-3p expression between TS patients and NCs.
We found a positive correlation and a linear relationship between Ct values for miR-140-3p and YGTSS scores of all study participants (Figure 2). Lower serum levels of miR-140-3p, which we observed in TS+ASD patients, could be potentially associated with both occurrence and worsening of motor and phonic tics, whereas higher serum levels of miR-140-3p, which we found in ASD patients, could be linked to the absence of tics. This finding indicates that expression analysis of serum miR-140-3p could strengthen the clinical diagnostic process of either ASD or TS+ASD. Moreover, this result gives strength to the hypothesis that the presence of phonic and motor tics, determining the comorbidity of TS with ASD, may be responsible for the different levels of serum miR-140-3p between ASD and TS+ASD patients.
Through network functional analysis, we observed that the regulatory network mediated by miR-140-3p is partly involved in managing structural and functional integrity of the nervous system and in the development of several human systems and organs (Table 3, Supplementary Figure 5). In particular, our computational data showed that CD38 and NRIP1, validated targets of miR-140-3p, take part in a set of biological processes convergingly dysregulated in ASD, like synaptic plasticity, immune response, and chromatin binding (Voineagu and Eapen, 2013;Gokoolparsadh et al., 2016;Ansel et al., 2017; Tables 4, 5, Supplementary Figures 6, 7).
CD38 was initially identified as an activation marker of immune cells, but it is now considered a virtually ubiquitous multifunctional molecule, involved in signaling and cell homeostasis. It is highly expressed in the brain, particularly in the hypothalamus (Quarona et al., 2013). Thanks to its ADP-ribosyl cyclase activity, CD38 regulates the mobilization of calcium ion from intracellular stores and therefore, it is mainly involved in proliferation, contraction, and secretion. It has also been found in a soluble form, maintaining this enzymatic activity, in body fluids and in exosomes (Quarona et al., 2013). CD38 has been linked to HIV infection, cancer, type 2 diabetes mellitus, and asthma (Quarona et al., 2013). Its enzymatic activity is responsible for the secretion of oxytocin and makes it one of the principal regulators of the social brain (Jin et al., 2007). Other than in social behavior, CD38 plays a role also in hippocampus-dependent learning and memory (Kim et al., 2016) and in postnatal glial development (Hattori et al., 2017). It is not surprising that much evidence has linked CD38 to ASD. Two genetic variants of CD38, the intronic SNP rs3796863 and the common Japanese SNP rs1800561 (causing the R140W mutation), have been associated with ASD (Munesue et al., 2010). In a study on two young sisters with ASD, a deletion of 4p15.32, resulting in a BST1 (bone marrow stromal cell antigen 1)-CD38 fusion transcript and in disruption of CD38 expression, was identified only in the girl affected by more severe ASD and asthma (Ceroni et al., 2014). CD38 expression is markedly reduced in LBC derived from ASD patients compared to their unaffected parents (Lerer et al., 2010); all-trans retinoic acid can upmodulate CD38 expression in these cells (Riebold et al., 2011). Finally, ASD patients have a higher absolute number of B cells per volume of blood and number of B cells expressing the cellular activation marker CD38 (Ashwood et al., 2011). NRIP1 (also known as RIP140, receptor-interacting protein 140) is a widely expressed, multifaceted transcription coregulator. Its primary physiological action is to trigger hormonecontrolled gene suppression (Nautiyal et al., 2013). NRIP1 mainly controls female fertility in the ovary, promoting ovulation, and energy homeostasis in metabolic tissues. It exerts a co-activator function in the regulation of circadian rhythms, inflammatory cascade, and mammary gland development (Nautiyal et al., 2013). NRIP1 is also expressed in the cortical and hippocampal areas of the brain (Nautiyal et al., 2013). It plays a crucial role in brain development and functioning and in cognitive and emotional processes (Duclot et al., 2012;Flaisher-Grinberg et al., 2014;Feng et al., 2015). Downregulation of its nuclear form controls mice brain aging (Ghosh and Thakur, 2009) whereas its cytosolic form acts as a neuroprotector in mice brain, preventing endoplasmic reticulum stress-induced neuronal apoptosis  and maintaining brain cholesterol homeostasis (Feng et al., 2015). In mice hippocampus, increased levels of NRIP1 expression are associated with depression-like symptoms (Chunhua et al., 2016). NRIP1 protein levels are considerably increased in the hippocampus from Down Syndrome patients (Gardiner, 2006). NRIP1 has been suggested as a potential candidate gene in autism, on the basis of in silico analysis of chromosomal regions involved in an unbalanced rearrangement del(21)(q11.2q21.2), identified in an ASD patient from a cohort of 126 ASD patients through high-resolution comparative genome hybridization (Iurov et al., 2010). Its potential etiological role in ASD has not been further investigated, but evidence suggests that NRIP1 and ASD could be indirectly linked. NRIP1 represents a molecular bridge between circadian rhythms and metabolism. It is part of a feedback mechanism regulating the circadian clock: it is under circadian regulation and it can alter basal levels of other clock genes, also by acting as a coactivator for the nuclear receptor RORα, known to be a stimulator of clock genes' transcription (Poliandri et al., 2011). RORA has been identified as dysregulated in many tissues from ASD patients (Cook et al., 2015). In general, mutations affecting the function of circadian-relevant genes are more frequent in ASD patients than in unaffected controls . All these findings give strength to other results that have previously linked ASD and circadian clock genes, leading to the interpretation of ASD as a neurodevelopmental disorder arising from atypical biological and behavioral rhythms : NRIP1 contribution to this association is worthy of further investigation.
When studying circulating miRNAs in pathologies, the biggest challenge is to elucidate the relationship between the diseased tissue and the corresponding expression levels of these molecules observed in liquid biopsies. MiRNAs in circulation could either be passively released non-specific by-products of cellular activity and cell death or actively secreted cell-cell signaling messengers (Ragusa et al., 2014Larrea et al., 2016). In ASD, this challenge is further complicated by the fact that a specific and unique diseased tissue has not been identified yet. In a gene expression dataset obtained from the analysis of whole blood samples of ASD patients, we found a marked dysregulation of twelve node genes from miR-140-3p-mediated network (Figure 3, Supplementary File 2). Even though we have detected the differential expression of some network node genes in all six ASD datasets reported, no one of those genes showed a marked and consistent expression trend in two or more datasets. Although it is not unexpected that independent microarray studies, using different technologies and platforms, give inconsistent results, we could not demonstrate a striking involvement of dysregulation of miR-140-3p-mediated network in ASD. Given the expression of miR-140-3p, CD38, and NRIP1 in the brain, it is also conceivable that brain tissues are responsible for serum dysregulation of miR-140-3p. Focusing on the contribution of serum extracellular vesicles to the expression of circulating miR-140-3p in ASD patients might help clarify potential tissue-serum links (Witwer, 2015). Nevertheless, our computational analysis of the potential functional role of intracellular miR-140-3p and of its possible involvement in ASD etiology suggests to the scientific community new processes, molecules and mechanisms to further investigate in the context of ASD.
Through ROC curve analyses and performance evaluation of predictive models, we proved that serum levels of miR-140-3p could be used in the discrimination of ASD patients from NCs, TS patients, and TS+ASD patients (Figures 4-6). We obtained the higher performance of serum miR-140-3p as a biomarker for the discrimination among ASD and TS+ASD patients. It was crucial to evaluate the performance of the biomarker through CV and permutation testing, since these predictive models were based on a miRNA already identified as differentially expressed between compared groups. We suggest that serum miR-140-3p could serve as a potential non-invasive biomarker to complement and support the behavior-based diagnosis of ASD, especially the differential one between ASD and TS+ASD. The main limit of our biomarker analysis is that we had only one miRNA to test for its predictive accuracy. It is well known that a single biomarker can hardly be as performing as a combination of them.
Further studies on larger cohorts and on participants of lower age would be necessary in order to get compelling evidences on miR-140-3p discriminatory power and prove its value in supporting early diagnosis of ASD. In this context, it would be interesting to identify which factors can be responsible for ASD patient misclassification (Figures 6 1-3D) by miR-140-3p, in order to optimize its predictive performance.
Our study, just like many others on circulating miRNAs to be used as biomarkers, has some limitations. The first one regards diagnostic specificity. As reported above, circulating miR-140-3p has already been associated with multiple pathological conditions and this denote that it could be simply indicative of a general disease state (i.e., inflammation and response to stress). The second one is reproducibility. There is little overlap between circulating miRNAs reported as biomarkers from independent investigators and this challenges their clinical utility. Just one of two other independent studies on circulating miRNA in ASD (discussed below) is consistent with ours. That is why results should be validated in larger cohorts and experimental conditions should be carefully standardized (Witwer, 2015).
Our study is the third high-throughput one profiling circulating miRNAs in a body fluid from ASD patients in order to discover some potential biomarkers.
The first study (Mundalil Vasu et al., 2014) was carried out on serum from a Japanese cohort of 55 ASD patients and 55 unaffected controls. The authors identified and validated 13 circulating miRNAs as dysregulated in serum from ASD patients and showed the accurate predictive power of 5 of them in discriminating ASD patients (Mundalil Vasu et al., 2014). None of circulating miRNAs from this study matches those from our profile. Our ASD and NC sample size is smaller than the one from this work, but we have used an array technology that allowed us to profile the expression of many more miRNAs than the 139 that the authors analyzed.
We tested the expression of these 5 predictive miRNAs (miR-19b-3p, miR-130a-3p, miR-181b-5p, miR-320a, miR-572) (Mundalil Vasu et al., 2014), together with miR-429 from our previous investigation on TS (Rizzo et al., 2015), in sera from 15 ASD patients, 15 TS patients, and 15 NCs. We did not observe any difference in miRNA expression among groups with the exception of miR-429, which we confirmed as DE between TS patients and NCs. Differences in sample size (55 samples per group vs. 15 samples per group) did not alter the result on miR-429 expression (Rizzo et al., 2015). However, functional enrichment analyses from both studies (Mundalil Vasu et al., 2014; this paper) demonstrated over-representation of the same neurological pathways, as TGF-β signaling, Hedgehog signaling, Wnt signaling, and regulation of synaptic plasticity. This observation suggests that discrepancies can be explained with differences in pre-analytic variables, such as genetic structure of studied populations, cohort composition, sample processing, validation technique, and data normalization (Witwer, 2015). Ethnicity of participants, cohort size, and miRNA panel and intercalating dye-based system used may have determined the inconsistencies observed.
The second paper (Hicks et al., 2016) describes a pilot study on saliva from a US cohort of 24 ASD patients and 21 unaffected controls, whose results partly match with ours. By RNA-sequencing, the authors identified 14 circulating miRNAs as dysregulated in saliva from ASD patients and showed the discriminative accuracy of this molecular signature (Hicks et al., 2016). MiR-140-3p is part of this ASD molecular fingerprint and is upregulated in ASD patients compared to unaffected controls, as in our study (Hicks et al., 2016;this paper). Moreover, in agreement with our results, their functional enrichment analysis detected significant over-representation of target genes related to neuronal development and transcriptional activation (Hicks et al., 2016). Our ASD and NC sample size is slightly bigger than the one from this work. Also, sequencing data from it have not been validated through miRNA-specific qPCR assays. Discrepancies between this work and our study can also be explained, other than with all the factors listed above, with the fact that they investigated two different human body fluids. It is interesting that miR-140-3p shows the same expression trend in both saliva and serum. This observation may suggest that shared mechanisms could determine the increased levels of miR-140-3p in both body fluids.
These two studies, differently from ours, have identified more than one biomarker for ASD. Furthermore, they have focused only on ASD patients and NCs, whereas our study is the first high-throughput one profiling circulating miRNAs also in patients suffering from another neurodevelopmental disorder comorbid with ASD, TS+ASD patients.
Through the identification of a serum biomarker, our study provides insight into concealed molecular mechanisms determining ASD and a potential complementary and supportive mean for a simpler, faster, and unbiased ASD diagnosis. The network that miR-140-3p regulates is involved in a set of biological processes convergingly dysregulated in ASD. Molecular characterization of miR-140-3p network would contribute to further clarify the heterogeneous molecular basis of ASD. Moreover, serum miR-140-3p could potentially be used as a non-invasive biomarker for ASD, easy to test through liquid biopsies.

ETHICS STATEMENT
This study was carried out in accordance with the recommendations of the local Ethics Committee (Comitato Etico dell' Università degli Studi di Catania, CEU) with written informed consent from all subjects.

AUTHOR CONTRIBUTIONS
MP and RR conceived the project with contributions by CDP, MR, DB, and MC. MP, CDP, MR, DB, and MC designed the experiments. MC and CB performed them. RR, MG, CND, and RB carried out patients' recruitment and clinical data analysis. RR performed clinical diagnosis. MC performed computational analysis and carried out statistical data analysis. MP, RR, and MC wrote the paper. All authors contributed to the critical revision of the data, read and approved the final manuscript.

FUNDING
This project was financially supported by Ministero dell'Università e della Ricerca Scientifica e Tecnologica (MIUR) and by Università degli Studi di Catania.

ACKNOWLEDGMENTS
The participation of Dr Duilia Brex to the final phase of the project is acknowledged. We thank the Scientific Bureau of the University of Catania for kindly editing the paper. We thank patients, unaffected controls and their families for accepting to take part in this project.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnmol.

2017.00250/full#supplementary-material
Supplementary Figure 1 | MiR-140-3p-mediated regulatory network. Nodes representing MIR140 gene and direct target genes of miR-140-3p are depicted as red diamonds; those symbolizing protein-coding genes are drawn as yellow circles; those standing for miRNA-coding genes are illustrated as orange triangles. Size of both nodes and node labels, except for CD38, MIR140, and NRIP1, is directly proportional to node degree within the network Supplementary Figure 2 | Over-represented GO terms in miR-140-3p-mediated regulatory network compared to the entire genome. (A) Dot plot for GO Biological Process terms. (B) Dot plot for GO Molecular Function terms. (C) Dot plot for GO Cellular Component terms. Each circle in the plots symbolizes an over-represented term: its x-axis coordinate reflects the Gene Ratio value; its size is directly proportional to the Count value; its color represents the Benjamini-Hochberg adjusted p-value generated by the hypergeometric test. Gene Ratio: the ratio of number of genes of interest that are annotated with a certain term from the database used to perform the analysis to number of genes of interest that are annotated with terms from the same database. Count: number of node genes within the network that are annotated with a certain term. GO: Gene Ontology Supplementary Figure 5 | Over-represented GO Biological Process terms regarding nervous system and development in miR-140-3p-mediated regulatory network. Dot plot for all the most interesting GO Biological Process terms regarding nervous system and development. Each symbol in the plot symbolizes an over-represented term: its x-axis coordinate reflects the Gene Ratio value; its y-axis coordinate reflects the Background Ratio value; its color represents the Benjamini-Hochberg adjusted p-value generated by the hypergeometric test. Gene Ratio: the ratio of number of genes of interest that are annotated with a certain term from the database used to perform the analysis to number of genes of interest that are annotated with terms from the same database. Background Ratio: the ratio of number of genes in the genome that are annotated with a certain term from the database used to perform the analysis to number of genome genes that are annotated with terms from the same database. GO: Gene Ontology.
Supplementary Figure 6 | Over-represented GO, DO and KEGG terms associated with CD38 in miR-140-3p-mediated regulatory network. Dot plot for all the most interesting GO, DO and KEGG terms whose enrichment is determined by CD38. For an in-depth description of the plot see Supplementary Figure 5 legend. For more information about what symbol shapes stand for see figure legend. GO terms whose enrichment is determined by both CD38 and NRIP1 are represented by bigger symbols. GO: Gene Ontology; DO: Disease Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes. Figure 7 | Over-represented GO and Reactome terms associated with NRIP1 in miR-140-3p-mediated regulatory network. Dot plot for all the most interesting GO and Reactome terms whose enrichment is determined by NRIP1. For an in-depth description of the plot see Supplementary Figure 5 legend. For more information about what symbol shape stands for see figure legend. GO terms whose enrichment is determined by both CD38 and NRIP1 are represented by bigger symbols. GO: Gene Ontology.

Supplementary
Supplementary File 1 | MiRNA profiling raw data: Ct values for TLDA panel A and B. Customized normalization approach.
Supplementary File 2 | Data from network expression analysis.
Supplementary File 4 | Ct values for miR-140-3p and neuropsychiatric scores of all participants.
Supplementary File 5 | Data from network functional analysis.