Potential Associations Among Alteration of Salivary miRNAs, Saliva Microbiome Structure, and Cognitive Impairments in Autistic Children

Recent evidence has demonstrated that salivary molecules, as well as bacterial populations, can be perturbed by several pathological conditions, including neuro-psychiatric diseases. This relationship between brain functionality and saliva composition could be exploited to unveil new pathological mechanisms of elusive diseases, such as Autistic Spectrum Disorder (ASD). We performed a combined approach of miRNA expression profiling by NanoString technology, followed by validation experiments in qPCR, and 16S rRNA microbiome analysis on saliva from 53 ASD and 27 neurologically unaffected control (NUC) children. MiR-29a-3p and miR-141-3p were upregulated, while miR-16-5p, let-7b-5p, and miR-451a were downregulated in ASD compared to NUCs. Microbiome analysis on the same subjects revealed that Rothia, Filifactor, Actinobacillus, Weeksellaceae, Ralstonia, Pasteurellaceae, and Aggregatibacter increased their abundance in ASD patients, while Tannerella, Moryella and TM7-3 decreased. Variations of both miRNAs and microbes were statistically associated to different neuropsychological scores related to anomalies in social interaction and communication. Among miRNA/bacteria associations, the most relevant was the negative correlation between salivary miR-141-3p expression and Tannerella abundance. MiRNA and microbiome dysregulations found in the saliva of ASD children are potentially associated with cognitive impairments of the subjects. Furthermore, a potential cross-talking between circulating miRNAs and resident bacteria could occur in saliva of ASD.


Introduction
Autism Spectrum Disorder (ASD) is a heterogeneous group of complex neurodevelopmental disorders characterized by impaired social communication and the presence of stereotyped and repetitive behaviors [1]. The etiopathogenesis of ASD is still unclear but is believed to be the complex result of a combination of genetic, epigenetic and environmental factors [2][3][4]. ASD is genetically highly heterogeneous. Both inherited and de novo ASD-associated variants have been characterized in hundreds of genes [5][6][7][8][9][10][11][12]. Immune dysregulation and gastrointestinal abnormalities are of particular interest in the light of several papers reporting ASD-associated disturbances in the peripheral, enteric and neuro-immune systems [13,14]. This association of ASD with a great prevalence of immune and gastrointestinal dysregulations led to investigations on the gut microbiome of ASD patients, which is emerging as a key regulator of intestinal physiology, neuroimmunity, and host behavior. Intriguingly, reports on gnotobiotic animals and probiotic studies have shown that microbiome dysregulation can directly induce behavioral and neuropathological endophenotypes of human ASD [15][16][17][18][19]. The gut microbiota represents a barrier against the proliferation of pathogenic organisms [20], playing a key role in the functioning of the host immune system [21], modulating gene expression [22], and reducing inflammation [23]. Alterations of brain structure and function development are related to modifications of the gut microbial composition because the interactions between this and the Central Nervous System (CNS) are already established during fetal life and maintained in adults [24]. Despite several reports about gut microbiota dysbiosis in ASD patients, there is little consensus across independent studies on specific bacterial species that are similarly altered [25]. Moreover, whether microbiome alteration is caused by ASD symptoms or whether it contributes to the ASD gut phenotypes is still unclear. While mounting evidence suggests a key role for the gut microbiota in ASD, the etiopathogenetic contribution of microorganisms living in the oral cavity has not been satisfactorily explored. The oral cavity represents the major gateway of bacteria into the human body: it is a unique and complex habitat [26], which harbors approximately 700 predominant taxa [27]. In addition, biologically relevant interactions between the saliva microbiome and other microbiomes in the human body may occur [28]. Several studies have shown that saliva microbiome is dysregulated in patients affected by systemic diseases, such as liver cirrhosis, diabetes, rheumatoid arthritis, and cancer [29][30][31][32][33]. The oral microbiome is modulated in response to several intrinsic and extrinsic factors that may in turn affect brain functionality [34]. The oral cavity is considered an important source of alterations that manifest at distant body sites, including the neural system [27]. Indeed, some studies reported the effects of the oral microbiome on neural functions. Accordingly, oral dysbiosis has been associated with Parkinson's disease, Alzheimer's disease, multiple sclerosis, and migraine [35][36][37][38][39]. Saliva also represents an important reservoir of molecules, including proteins and microRNAs (miRNAs), associated with the pathological development of neurological diseases and brain dysfunctions [40][41][42][43][44][45][46]. Recently, fourteen miRNAs were found differentially expressed in the saliva of ASD patients and showed significant correlations with Vineland neurodevelopmental scores [47]. Based on these premises, we report the miRnome and microbiome analysis from saliva of ASD children with the aim to verify the existence of a molecular correlation between miRNA dysregulation and salivary dysbiosis and their contribution to ASD pathogenesis.

Microbial Structure of the Saliva Microbiome in Children with ASD and NUCs
To evaluate changes in the microbial composition of the salivary microbiome, a total of 53 ASD and 27 NUC samples were sequenced using the Illumina MiSeq platform.
A total of 10,919,413 valid reads with an average of 127,464 reads per sample were generated (74,469-1,641,234 range). All sequences were assigned to 341 Operational taxonomic Units (OTUs) with at least 97% similarity level showing 10 phyla, 65 genera, and 86 species. The alpha diversity, Chao1 index (richness), Shannon index (diversity), and Shannon Evenness E index (evenness) revealed no significant differences between ASD and NUC groups (p > 0.05) although the Chao1 index showed a slight decrease in the ASD group ( Figure 3A). The overall dissimilarities of microbial community structure between the two groups were calculated by using the weighted UniFrac distance. The beta diversity shown by PCoA (principal coordinate analysis) plot on weighted (accounting for the abundance of OTUs) revealed no clustering differentiation of the bacterial communities in the two groups ( Figure 3B). The composition of the salivary microbiome between ASD and NUC was explored in terms of the relative abundances at different taxonomic levels through White's non-parametric t-test; p-value < 0.05 by STAMP software. Predominant phyla in ASD and NUC groups were Firmicutes (41.1% vs. 42.7%), Bacteroidetes (18.9% vs. 22.3%), Proteobacteria (16.9% vs. 10.8%), Actinobacteria (14.1% vs. 13.3%), and Fusobacteria (7.6% vs. 9.3%) constituting 98.6% of salivary microbes (see Supplementary Figure S1). Overall, Proteobacteria were more abundant in ASD patients compare to control participants and a slight increase was also observed in Actinobacteria, while the other phyla were more abundant in the NUC group.
A total of 65 genera were detected (see Supplementary Figure S2) and among these, 10 genera showed a statistically significant difference in relative abundance between ASD and NC children by applying the two-sided White's non-parametric t-test ( Figure 4A), and also confirmed by applying the Mann-Whitney-Kruskal-Wallis test and t-test/ANOVA ( Supplementary Table S1). and NUC groups were Firmicutes (41.1% vs. 42.7%), Bacteroidetes (18.9% vs. 22.3%), Proteobacteria (16.9% vs. 10.8%), Actinobacteria (14.1% vs. 13.3%), and Fusobacteria (7.6% vs. 9.3%) constituting 98.6% of salivary microbes (see Supplementary Figure S1). Overall, Proteobacteria were more abundant in ASD patients compare to control participants and a slight increase was also observed in Actinobacteria, while the other phyla were more abundant in the NUC group.  , Shannon E index (c, evenness) were calculated for saliva samples. The bars depict mean ± SD of relative abundance rates. NS, p > 0.05. ASD (salivary samples n = 53), NUC (salivary samples collected from healthy controls, n = 27). (B) β-diversity. PCoA plot generated using weighted UniFrac distances shows none differences between the two groups (ASD in red and NUC in blue).
In particular, Rothia, Filifactor, Actinobacillus, Weeksellaceae, Ralstonia, Pasteurellaceae, and Aggregatibacter increased their abundance rates in ASD, while Tannerella, Moryella, and TM7-3 decreased ( Figure 4A). At the species level, we found statistically significant differences for 11 taxa out of 86 species ( A total of 65 genera were detected (see Supplementary Figure S2) and among these, 10 genera showed a statistically significant difference in relative abundance between ASD and NC children by applying the two-sided White's non-parametric t-test ( Figure 4A), and also confirmed by applying the Mann-Whitney-Kruskal-Wallis test and t-test/ANOVA (Supplementary Table S1). In particular, Rothia, Filifactor, Actinobacillus, Weeksellaceae, Ralstonia, Pasteurellaceae, and Aggregatibacter increased their abundance rates in ASD, while Tannerella, Moryella, and TM7-3 decreased ( Figure 4A). At the species level, we found statistically significant differences for 11 taxa out of 86 species ( Figure 4B). Cliff's delta effect size values ranged between −0.4556 and 0.2613 and, according to the Cliff's statistic, could be considered small and mostly moderate.
For microbiome data, we also applied the negative binomial regression. The Negative Binomial (NB) model was positively assessed for seven (Tannerella, Weeksellaceae, Moryella, Filifactor, Ralstonia, Pasteurellaceae, and Actinobacillus parahaemolyticus) of the eleven microbiome taxa and revealed many statistically significant regressions shown in Figure 6. For each of the predictor variable, Supplementary Table S2 reports a regression coefficient (B) with a p-value (p < 0.05) resulting from the Wald Chi-Square test. Effect size measures, reported as Standardized mean difference (SMD) values, showed small and moderate effects (see Supplementary Table S2). Negative binomial regression (NBR) results were partially confirmed applying the multiple and simple linear regression models. These additional analyses confirmed the regression between Total Intelligence Quotient (TIQ), ADOS-A, ADOS-B and Tannerella; miR-141-3p and Weeksellaceae; miR-141-3p, ceruloplasmin and Moryella; PIQ and Filifactor; miR-141-3p, ADI-C, ADI-D, ammonium and Actinobacillus parahaemolyticus (see Supplementary Table S3). −0.34); while, Ralstonia showed a positive correlation with ADI-A (r = 0.34) ( Figure 5B). Moreover, computing the Spearman correlations between salivary DE miRNAs and bacteria, we detected a negative association between Tannerella and miR-141-3p (r = −0.30) ( Figure 5C). For microbiome data, we also applied the negative binomial regression. The Negative Binomial (NB) model was positively assessed for seven (Tannerella, Weeksellaceae, Moryella, Filifactor, Ralstonia, Pasteurellaceae, and Actinobacillus parahaemolyticus) of the eleven microbiome taxa and revealed many statistically significant regressions shown in Figure 6. For each of the predictor variable, Supplementary Table S2 reports a regression coefficient (B) with a p-value (p < 0.05) resulting from the Wald Chi-Square test. Effect size measures, reported as Standardized mean difference (SMD) values, showed small and moderate effects (see Supplementary Table S2). Negative binomial regression (NBR) results were partially confirmed applying the multiple and simple linear regression models. These additional analyses confirmed the regression between Total Intelligence Quotient (TIQ), ADOS-A, ADOS-B and Tannerella; miR-141-3p and Weeksellaceae; miR-141-3p, ceruloplasmin and Moryella; PIQ and Filifactor; miR-141-3p, ADI-C, ADI-D, ammonium and Actinobacillus parahaemolyticus (see Supplementary Table S3).

Functional Enrichment Analyses
Pathway enrichment analyses were computed to investigate the potential biological impact associated with the differential expression of the miRNAs reported in this study. We computationally identified a list of statistically over-represented gene ontologies and biological pathways (Figure 7) potentially related to ASD, as Pre-NOTCH transcription and Translation (False Discovery Rate (FDR) corrected p = 0.0015) and Pre-NOTCH Expression and Processing signaling pathway (FDR corrected p = 0.00471), signaling by NOTCH (FDR corrected p = 0.0221), signaling by NGF (FDR corrected p = 0.0291) through which cell fate decisions in neuronal development are regulated. Moreover, we found Parkinson's disease related pathways (FDR corrected p = 0.0002), involved in synaptic and mitochondrial dysfunction and neuroinflammation; Hippo signaling pathway (FDR corrected p = 3.991× 10 -6 ), FoxO (Forkhead box O) signaling pathway (FDR corrected p = 0.0001), PI3K-Akt signaling pathway (FDR corrected p = 0.0008), and mTOR signaling pathway (FDR corrected p = 0.02), involved in many biological events such as apoptosis, cellular stress and cell-cycle control. In addition, bacterial invasion of epithelial cell (FDR corrected p = 1.19 × 10 -6 ) related pathways were also found. These data would represent an interesting link between the differential expression of miRNAs and the bacterial abundance alterations found in saliva.
Int. J. Mol. Sci. 2020, 21, x FOR PEER REVIEW 10 of 24 regression coefficient is indicated by a color gradient from green (negative prediction) to red (positive prediction), as shown in the colored bar. Statistically significant p-values are indicated by asterisks.

Functional Enrichment Analyses
Pathway enrichment analyses were computed to investigate the potential biological impact associated with the differential expression of the miRNAs reported in this study. We computationally identified a list of statistically over-represented gene ontologies and biological pathways (Figure 7) potentially related to ASD, as Pre-NOTCH transcription and Translation (False Discovery Rate (FDR) corrected p = 0.0015) and Pre-NOTCH Expression and Processing signaling pathway (FDR corrected p = 0.00471), signaling by NOTCH (FDR corrected p = 0.0221), signaling by NGF (FDR corrected p = 0.0291) through which cell fate decisions in neuronal development are regulated. Moreover, we found Parkinson's disease related pathways (FDR corrected p = 0.0002), involved in synaptic and mitochondrial dysfunction and neuroinflammation; Hippo signaling pathway (FDR corrected p = 3.991× 10 -6 ), FoxO (Forkhead box O) signaling pathway (FDR corrected p = 0.0001), PI3K-Akt signaling pathway (FDR corrected p = 0.0008), and mTOR signaling pathway (FDR corrected p = 0.02), involved in many biological events such as apoptosis, cellular stress and cell-cycle control. In addition, bacterial invasion of epithelial cell (FDR corrected p = 1.19 × 10 -6 ) related pathways were also found. These data would represent an interesting link between the differential expression of miRNAs and the bacterial abundance alterations found in saliva. These computational data would suggest a functional involvement of DE salivary miRNAs in molecular signaling pathways related to the development of cognitive functions, often reported to be dysfunctional in ASD. This observation is particularly intriguing because it hints at the possibility that some molecular alterations of ASD neuronal circuits could be mirrored in saliva by means of differential secretion of miRNAs.

Circulating miRNAs and Microbiome Structure are Altered in Saliva of Pediatric ASD Patients
Recent studies have shown that perturbations in the normal structure of microbiome have a close relationship with human health and disease [19,31,32] and there is mounting evidence that alterations of the GI microbiome may influence neurological disorders and autistic behavior [48]. The shift of microbial communities in the oral cavity in these diseases is being increasingly studied to look at possible early diagnostic markers [25,35,36,39]. Several studies have revealed that miRNAs play very important roles in neural developmental processes [49][50][51][52][53], resulting in neuropsychiatric disorders [54,55] and neurodegenerative diseases [56,57], although their exact pathophysiologic role remains unclear [58]. In the present study, the salivary miRNA content and the oral microbiome were evaluated in a group of pediatric patients with ASD syndrome, compared with neurotypical subjects. MiRnome and microbiome dysregulation in saliva of ASD patients was previously reported by other groups, but there is only a partial concordance with the data presented here [47,[59][60][61]. The salivary microbiome structure in human beings is strongly affected by the geographical, climatic and ethnic origin of samples [62][63][64]. This could be due to the different ratio of dietary protein/ carbohydrate intake that, in turn, modulates salivary pH [65]. This biochemical variability, such as the inconstancy of desquamated oral epithelial cells in saliva among people, would lead to a heterogeneous population of miRNAs in saliva among individuals [66].
In the present paper, a higher relative abundance of Proteobacteria and Actinobacteria phylotypes and a lower amount of Bacteroidetes, Firmicutes, and Fusobacteria were found in ASD patients. In addition, an increased abundance of Rothia mucilaginosa, Filifactor, Actinobacillus parahaemolyticus. parahaemolyticus, Weeksellaceae, Ralstonia, Pasteurellaceae, Agregatibacter segnis, and Haemophilus parainfluenzae, as well as a reduced amount of Tannerella, Moryella and TM7-3 were observed in the saliva of ASD patients. These observations could contribute to define the dysbiotic signatures of disease. In particular, a similar trend for Haemophilus and Rothia is supported by a previous observation [60].These taxa are often involved in several diseases, especially in immunocompromised hosts [67,68]. In the present study, many groups of microorganisms differentially abundant in ASD children were previously associated with CNS (Central Nervous System) disorders, such as Parkinson's and Alzheimer's Disease [69,70]. In some cases, these microorganisms are also relevant pathogens of periodontal diseases. It is noteworthy that Filifactor and Tannerella were previously reported as strong indicators of dysbiosis [71,72].
Previous papers reported that the five DE miRNAs we identified in ASD saliva were associated to neurological diseases, including ASD, in cellular and extra-cellular contexts. In Alzheimer's patients, miR-29a was upregulated in cerebrospinal fluid (CSF) [73], whereas, miR-451a and miR-16-5p were downregulated in exosomes from CSF of young-onset (YOAD) subjects [74] and miR-141 decreased in the plasma fraction enriched in exosomes [75]. Moreover, the levels of miR-29a-3p were decreased in whole blood [76] and those of miR-141-3p in serum of Parkinson's patients [77]. Levels of let-7b-5p were decreased in the saliva of children with mild traumatic brain injury [78]. Intriguingly, a downregulation of miR-451a and miR-16-5p was also identified in peripheral blood of ASD children [79]. Furthermore, miR-451a and miR-16-5p were upregulated and downregulated, respectively, in lymphoblastoid cell lines from autistic monozygotic twins [80]. Finally, miR-451 also showed an increased abundance in post-mortem ASD brain [81]. These findings would suggest a critical role for these miRNAs in the development of cognitive functions. Molecular pathways identified in our functional enrichment analysis would support the involvement of these miRNAs in ASD related mechanisms. Dysregulation of IGF-I/PI3K/AKT/mTOR signaling is associated with ASD pathogenesis by affecting the myelination, synaptic plasticity, mechanisms of social interactions, learning and immune functions [10,82,83]. The Hippo signaling pathway is related to neuropsychiatric disorders such as schizophrenia, bipolar disorder, obsessive-compulsive disorder and ASD, playing an important role in neural development and neuronal maintenance [84]. Finally, NOTCH, FoxO, and NGF pathways play multiple roles in neurodevelopment in the CNS regulating neurogenesis, synaptic plasticity, learning, memory, and behavior [85][86][87].

Potential Associations Among Cognitive Impairments, Salivary miRNA Expression and Microbiome Alteration in ASD Children
We found several linear associations between salivary miRNA expression and anomalies of social interaction and communication, expecially for let-7b, miR-451a, and miR-29a-3p ( Figure 5A). Impairments in social-emotional reciprocity and relationships, verbal and nonverbal communicative behavior, cognitive skills, lead to increase the risk of social isolation and rejection of ASD children in daily social contexts. Some studies reported an association between social isolation and consequent psychological stress and miRNA dysregulation. For instance, miR-141-3p expression progressively increases in a time-dependent manner in mouse models of post-stroke social isolation [88], and members of the miR-29 family increased their expression during the healing process of oral palatal mucosal wounds [89]. Moreover, decreasing plasmatic let-7b-5p correlated strongly with cognitive impairment in the presence of severe alcohol use disorder [90] and it was also related to high stress conditions [91]. Levels of circulating miR-451a in the blood of depressed patients were found reduced and negatively associated with the Hamilton Depression Scale that is used to rate the severity of depression: a depression-like phenotype is often observed in autistic patients [92,93].
In blood of ASD patients elevated levels of lactate and its synthesizing enzyme, lactate dehydrogenase (LDH), reflect the mitochondrial energy metabolism dysfunctions [94]. Interestingly, miR-141-3p was found to induce mitochondrial dysfunctions in obese mice by inhibiting PTEN [95]. Moreover, circulating miR-141-3p positively correlated with LDH levels in rectal cancer [96], as well as, miR-141 positively regulated expression of LDH by inhibiting MAP4K4 in breast cancer [97]. These observations corroborate the positive correlation we found between miR-141 and lactate and suggest that it could represent a combined extracellular phenotype of ASD metabolic abnormalities.
Moryella showed a signature that negatively correlated with anomalies in neurodevelopment arising before the age of 36 months (ADI-D); the same genus was also linked to the Verbal Intelligence Quotient. Ralstonia positively correlated with a worst Qualitative anomaly in social interaction (ADI-A). Finally, the low abundance of phylotype TM7-3, often associated with human inflammatory mucosal diseases [98], negatively correlated with the increasing of both Qualitative anomalies in communication (ADI-B) and Repetitive and restricted behavior (ADI-C).
The negative binomial regression analysis revealed that the abundance of seven species was significantly related to cognitive impairments, especially for Tannerella [99] (Figure 6). In particular, VIQ, PIQ, and TIQ together with behavior (ADI-C), speech and communication anomalies (ADOS-A) were significant predictors of Tannerella abundance. In accordance with these results, it has been reported that periodontal disease is related to cognitive decline [100], disability, speech and communication impairment, low self-esteem and quality of life [101]. Interestingly, the regression analyses uncovered a significant relationship among all the Intelligence Quotients (i.e., VIQ, PIQ, and TIQ) and Weeksellaceae and Ralstonia abundances.
These results pave the way to suggestion that miRNA expression and microbiome structure alteration could be a consequence of ASD symptomatology.

Potential Crosstalk between miRNAs and the Microbiome in Saliva
It is known that microbiota can secrete bioactive compounds able to modify the host epigenome; as well as, miRNAs from the host could selectively regulate the functions of microbiota.
Linear correlation analysis between miRNA expression and microbiome data in saliva led to the identification of a negative relationship between miR-141-3p and Tannerella. This inverse correlation has already been reported. It has been demonstrated that miR-141-3p was underexpressed in the gingiva of patients affected by periodontal disease compared to healthy gingival tissue [102], while, Tannerella is one of the major Gram-negative periodontal pathogens and it is already identified as a marker of autism [60,99]. The negative binomial regression sheds a light to other considerations regarding the potential interaction between: (i) miRNAs and neuropsychological parameters (as predictors) and (ii) microbiome (as outcome variables) in terms of expected change in the outcome variable for a one-unit change in the predictor variable. For example, NB regression revealed new potential interactions between the miRNAs previously associated with cognitive impairments (based on linear correlation analyses) and specific bacteria: let-7b-5p and miR-16-5p with Tannerella; miR-451a and Weeksellaceae; miR-141-3p and Tannerella, Weeksellaceae, Moryella, Filifactor ( Figure 6). MiR-141-3p is known to be expressed in intestinal epithelial cells and a potential biomarker for gut dysbiosis [103].
Although this issue is quite unexplored in the literature, recent reports demonstrated that fecal miRNAs could contribute to shaping the composition of the gut microbiome [104,105], suggesting a mechanism by which host cells can regulate the microbial community.
The significant relationships found in the regression analysis may suggest that the potential crosstalk between oral bacteria and salivary miRNAs could be based on a general response of host to dysbiosis involving several miRNA families, as well as bacterial phyla [106][107][108][109][110][111][112][113].
Taken together, all these findings suggest that in the saliva of autistic children quantitative differences in miRNA expression and bacteria abundance are present and potentially associated with anomalies in social interaction and communication.

Ethics Approval and Consent to Participate
All experiments were approved by the local Ethics Committee, Comitato Etico Catania 1, University of Catania (ID: 002430-36) prior to sample collection. Written, informed consent was obtained from all parents and each participant who were able gave their informed assent.
All experimental methods were in accordance with Helsinki Declaration.

Participant Selection
From a database of more than 2000 patients, 76 treatment-naïve patients affected by ASD were recruited and studied from January to October 2018 in the outpatient service of the Child and Adolescent Psychiatry Unit (Department of Clinical and Experimental Medicine, University Hospital of Catania).
The inclusion criteria were clinical diagnosis of ASD, according to the criteria of the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-5, APA 2013), and the absence of other medical, neurological, genetic or metabolic condition such as epilepsy, cytogenetically visible chromosomal abnormalities, copy number variants (CNVs) or single-gene disorders. They were compared to 39 neurologically unaffected controls (NUC), without any history of ASD and who suffered from neither chronic neurological, metabolic or genetic diseases nor psychiatric disorders. All participants to the study were Caucasians from Sicily, randomly recruited from various socio-economic contexts.
The entire cohort was split into 2 sets of samples: (a) discovery set, composed by 23 ASD and 12 NUC; (b) validation set, made by 53 ASD and 27 NUC. The discovery set was used to perform expression profiling by NanoString technology, while the validation set was analyzed by real time PCR single assays and 16S rRNA microbiome screening.

Assessment
All participants were assessed by a child and adolescent neuropsychiatrist expert in ASD (RR) with the following instruments: WISC-III (Wechsler Intelligence Scale for Children, III edition) [114] or WPSSI (Wechsler Preschool and Primary Scale of Intelligence) [115] as an evaluation of both IQ (Intelligence Quotient) and cognitive functioning, ADOS (Autism Diagnostic Observation Schedule) and ADI-R (Autism Diagnostic Interview-Revised) to evaluate ASD symptoms.
Based on these scales and schedules, the assessment procedure was carefully conducted assigning specific scores which provided a measure of autism severity. Neuropsychological features of participants are summarized in Table 1.

Sample Collection
All participants were instructed to refrain from eating or drinking for at least 3 h prior to saliva collection. Saliva samples were collected in a good status of oral hygiene (i.e., brushing teeth once/twice daily), without any tooth decay. Saliva collection was always performed between 8:30 and 10:30 a.m. to avoid any potential microbiome and miRnome oscillation due to the circadian rhythms. From a minimum of 800 µL to a maximum of 4 mL of non-stimulated and naturally outflowed saliva was collected into 50 mL conical centrifuge tubes. Saliva samples were centrifuged at 10,000 rpm for 15 min at 4 • C to separate the pellet and supernatant for microbiological analysis and miRNA expression assays, respectively. The pellets were immediately processed, while supernatants were aliquoted into 2 mL RNase-free tubes and stored at −80 • C until analysis.
To test the hematological parameters of each study participant, all participants were instructed to refrain from eating or drinking for at least 3 h prior to blood collection, which was performed between 8:30 and 10:30 a.m. Peripheral blood samples were collected through a butterfly device into a 5 mL collection tube. Collection tubes were treated according to current and standard procedures for clinical samples.

RNA Extraction
Extraction of total RNA was carried out from 800 µL of saliva samples using Qiagen miRNeasy Mini Kit (Qiagen, GmbH, Hilden, Germany), according to Qiagen Supplementary Protocol for purification of RNA (including small RNAs) from serum or plasma [46,116]. RNA was eluted in 200 µL RNAse-free water and then precipitated by adding 20 µg glycogen, 0.1 volumes 3 M sodium acetate and 2.5 volumes ice-cold 100% ethanol. After incubation at −80 • C overnight, RNA was centrifuged and washed twice in ice-cold 75% ethanol and resuspended in 7 µL RNAse-free water. The yield and quality of the RNA samples were assessed by using NanoDrop Lite Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA).

MiRNA Profiling by NanoString Technology
To profile the expression of circulating miRNAs from saliva, the NanoString nCounter system assays were performed using the NanoString platform and the nCounter Human v3 miRNA Expression Assay Kits (NanoString Technologies, Seattle, WA, USA), according to the manufacturer's instructions. MiRNA profiling was performed on 23 ASD patients and 12 NUCs, starting from 3 µL of isolated RNA (approximately 150 ng). Samples were processed using the automated nCounter Prep Station; following hybridization, they were purified and immobilized on a sample cartridge for quantification and data collection by using the nCounter Digital Analyzer. The nSolver 3.0 software was used for data analysis. The endogenous control was selected from the arrays in a similar way to the GMN (global median normalization) method [55,117]. By this approach, we identified miR-21-5p as the best endogenous control for our experimental model.

Processing and Analyses of Sequencing Data
QIIME pipeline (Quantitative Insights into Microbial Ecology) v.1.9.1 was used to process the generated raw FASTQ files [120]. V3-V4 16S rRNA FASTQ were de-multiplexed using the barcodes. The paired-end sequences were assembled to form a single read using FLASH [121] and quality-filtered ≥80% bases in a read above Q30 (see Supplementary Table S4). Merged reads were length-filtered based on 445 bp (the expected length). The ends of retained (not-merged) forward reads were clipped to a total read length of 270 bp to remove low quality bases. The high-quality reads were clustered against a reference sequence collection with QIIME. To focus only on the prominent taxa, a filtering step of 0.01% at the Operational taxonomic Unit (OTU) level was performed by running a workflow on QIIME (filter_otus_from_otu_table.py). The taxonomy of each 16S rRNA gene sequence was collapsed to OTUs using the open reference-based OTU picking method against Greengenes database at 97% of sequence similarity [122]. Chimeras were identified and removed by Chimera Slayer and the UCHIME algorithm. Any reads that did not match the reference sequence collection were subsequently clustered de novo. To avoid sample size biases in downstream analyses, rarefaction curves were generated with QIIME (alpha_rarefaction.py workflow) and calculated by applying Explicet and a maximum depth of 74.469 sequences/sample [123]. The OTU tables were used for assessing α-diversity indices (Chao-1, Shannon and Shannon evenness) calculated from the taxonomic profiles and compared across the ASD and NUC groups by QIIME algorithms. Independent Student's t-test and Mann-Whitney U test were used to evaluate α-diversity among the taxonomic profiles and compared across the ASD and NUC groups. β-diversity between ASD and NUC groups was analyzed by weighted and unweighted UniFrac distance matrices (beta_diversity.py workflow) and visualized through tridimensional PCoA plot by using EMPeror (http://boocore.github.io/emperor/). The differentially abundant OTUs across two sample categories were identified by QIIME scripts (differential_abundance.py). The core microbiome was determined by QIIME algorithms (compute_core_microbiome.py) and the diversity analysis was performed with the script core_diversity_analyses.py.
DNA sequences were deposited in the Sequence Read Archive under BioProjects PRJNA518756 and PRJNA518760.

Correlation and Negative Binomial Regression Analyses
To investigate whether a linear relationship exists between the differential expression of salivary miRNAs, microbiome taxa and the neuropsychological and hematological parameters for each selected participant, correlation analysis was performed. Additionally, we applied a Negative Binomial (NB) regression to explore the potential effects of the clinical parameters and the miRNA expression (potential predictors) on the relative microbial abundances (outcome variables). In particular, relationships were evaluated by using generalized linear models (GLM) assuming a negative binomial distribution and log link function, since our response variable is represented by over-dispersed count data that didn't have an excessive number of zeros [124].

Computational Enrichment Analysis
To investigate the functional meaning and potential ASD association of the identified DE miRNAs, DIANA-mirPath v.3 web server [125] and miRNet tool [126] were used for pathway enrichment analysis from KEGG (Kyoto Encyclopedia of Genes and Genomes) and Reactome gene annotation databases.

Statistical Approach
For miRNA profiling analysis, SAM (Significance of Microarrays Analysis) statistical tests were computed by using MeV (Multi experiment viewer v4.8.1) statistical analysis software. We computed a two-class unpaired test, based on 100 permutations. Fold change (FC) values were obtained by calculating the ratio between the normalized count mean of each group. MeV was also used to generate a heat-map of identified DE-miRNAs.
Concerning microbiome analyses, OTU frequencies across sample groups were performed by the Kruskal-Wallis test. Statistical analysis of taxonomic profiles was performed using STAMP (Software Testing AMPlification) [127][128][129] by a two-sided White's non-parametric t-test. Extended error bar plots were produced by STAMP (White's non-parametric t-test and p-value < 0.05) showing the bacterial taxa with a significant difference (p-value < 0.05). Mann-Whitney-Kruskal-Wallis, t-test/ANOVA were used for confirmation (Supplementary Table S1).
The correction for multiple tests of high throughput analyses was performed by applying Benjamini-Hochberg FDR (False Discovery Rate).
The Mann-Whitney U test was applied to evaluate the differential expression of tested miRNAs between the two groups in the validation analysis (p-value < 0.05). Cliff's delta statistic [130,131] was used to estimate the non-parametric effect sizes and were calculated by using the Cliff's Delta Calculator by Macbeth et al. [132].
Correlation analyses were performed by applying Spearman test with two-sided p-values corrected for multiple comparisons by using the Bonferroni-Šídák approach. Statistical analyses were computed using GraphPad Prism version 7.00 for Windows (GraphPad Software, La Jolla, CA, USA). To determine the strength of associations between variables, the correlation coefficients themselves were interpreted as index of effect size [133].