Salivary microRNA profiling dysregulation in autism spectrum disorder: A pilot study

Introduction Autism spectrum disorders (ASD) are the most prevalent neurobiological disorders in children. The etiology comprises genetic, epigenetic, and environmental factors such as dysfunction of the immune system. Epigenetic mechanisms are mainly represented by DNA methylation, histone modifications, and microRNAs (miRNA). The major explored epigenetic mechanism is mediated by miRNAs which target genes known to be involved in ASD pathogenesis. Salivary poly-omic RNA measurements have been associated with ASD and are helpful to differentiate ASD endophenotypes. This study aims to comprehensively examine miRNA expression in children with ASD and to reveal potential biomarkers and possible disease mechanisms so that they can be used to improve faction between individuals by promoting more personalized therapeutic approaches. Materials and methods Saliva samples were collected from 10 subjects: 5 samples of children with ASD and 5 from healthy controls. miRNAs were analyzed using an Illumina Next-Generation-Sequencing (NGS) system. Results Preliminary data highlighted the presence of 365 differentially expressed miRNAs. Pathway analysis, molecular function, biological processes, and target genes of 41 dysregulated miRNAs were assessed, of which 20 were upregulated, and 21 were downregulated in children with ASD compared to healthy controls. Conclusion The results of this study represent preliminary but promising data, as the identified miRNA pathways could represent useful biomarkers for the early non-invasive diagnosis of ASD.


Introduction
Autism spectrum disorder (ASD) is a neurological and developmental disorder characterized by altered social interaction, restricted/repetitive behavior, and medical and mental health conditions that result in lifelong functional and social impairments. It is estimated that ASD is associated with intellectual disability in 70% of cases, with seizures in 30% of cases (Minshew and Williams, 2007). Furthermore, ASD leads to alterations in the circadian cycle: it has been reported that people affected by ASD present a reduced total sleep, with an increased sleep latency. Studies have shown a higher frequency or absence of circadian variation of melatonin and cortisol in ASD patients (Tordjman et al., 2015).
Autism spectrum disorder presents a multifactorial etiology, comprising both genetic and environmental factors (Hallmayer et al., 2011;Singh et al., 2013;Sehovic et al., 2020;Wu et al., 2020). A metanalysis published in 2016, including a cohort of 6,413 twin pairs, estimated that the heritability of ASD to be 74-93% (Tick et al., 2016). Similar values were recently reported by Wu et al. (2020), while moderate heritability values (37-38%) were obtained by Hallmayer et al. (2011) in a Californian study (2011). Both genetic and epigenetic factors may be involved, such as histone modification, DNA methylation, and non-coding sequences of RNA (miRNA) (Zhang et al., 2020). Parental age, in utero exposure to several substances including pollution and inflammation, and low birth weight, may also increase the risk (Hallmayer et al., 2011;Sehovic et al., 2020).
The diagnosis of ASD, according to DSM-5, is mainly based on the examination of behavioral characteristics. A person must present deficits in three areas of social interaction: socialemotional reciprocity, non-verbal communicative behaviors, and the ability to maintain and understand relationships (Hallmayer et al., 2011;Wu et al., 2020;Qin et al., 2022). Unfortunately, all these characteristics can be interpreted differently, especially in young patients, leading to delays in diagnosis. Objective hallmarks are therefore needed (Hicks et al., 2020).
The number of children with ASD has grown dramatically over the last 2 decades with a current estimated prevalence of approximately 1 in 60 worldwide, an increase of 154% and even higher figures in future (Sehovic et al., 2020;Salloum-Asfar et al., 2021;Qin et al., 2022). Although the etiology of ASD is still mainly unknown, major genetic aberrations have been recognized in only 25-30% of patients with ASD and are mainly related to pathological mechanisms that start during gestation (Wiśniowiecka-Kowalnik and Nowakowska, 2019). Even if ASD could be diagnosed as early as 18 months of age, on average the diagnosis is done at ∼60 months of age (van't Hof et al., 2021). Furthermore, as early diagnosis of ASD is often unreliable, effective therapy is usually delayed. This might have negative consequences, considering that the effectiveness of therapy often depends on earliness of treatment beginning (Hicks et al., 2020;Salloum-Asfar et al., 2021). Therefore, recent research on early diagnosis of ASD has focused on large screening programs that can be efficient and easy to implement on a population scale (Hicks et al., 2018a).
Macroscopically, among other brain abnormalities observed in patients with ASD, alterations have been described in Broca and Wernicke's areas, which are associated with language, and in the frontal lobe, superior temporal cortex, parietal cortex, and amygdala, which are strongly implicated in social behaviors (Ha et al., 2015). However, their use in the clinical practice and to monitor treatment is difficult and expensive as magnetic resonance is needed (Squarcina et al., 2021).
Numerous biomarkers have been studied, measured and analyzed in order to identify potential ADS indicators for both early diagnosis and pathogenesis (Ballini et al., 2020;Hicks et al., 2020;Isacco et al., 2021;Salloum-Asfar et al., 2021). ASD presents typical features of synaptic dysfunction and synaptopathy (Wu et al., 2020). Indeed, aberrant overexpression of the postsynaptic density 95 protein (PSD95), a member of the membrane-associated guanylate kinase that promotes synaptic stability in excitatory synapses, has been observed in patients with ASD. This alteration leads to increased dendritic spines and pathologic synaptic hyperconnectivity (Coley and Gao, 2018). Another important feature of synapses that have been found Frontiers in Neuroscience 02 frontiersin.org altered in ASD is SHANK (SH and multiple ankyrin repeat domains protein). SHANK proteins are scaffolding proteins that organize intermediate scaffold proteins. Their location is principally in the excitatory synapse, where they allow the synapse to develop and function properly. Alterations in the SHANK gene have been reported as a possible cause and predictive feature of ASD (Wu et al., 2020). Among others, miRNA have been proposed as plausible new biomarkers with high predictivity. More than 91 miRNAs associated with ASD have been identified. For instance, alterations of miRNA related to cellular respiration have been found in the serum of patients with ASD. Furthermore, miRNA-500a-5p and miRNA-197-5p expressions in the serum have been suggested as useful tools to diagnose ASD in children (Kichukova et al., 2021). Through analysis of miRNAs, genetic correlations were observed between clinically unaffected parents of children with ASD and their siblings (Ozkul et al., 2020). A great number of miRNA are altered in ASD, and several are related to alteration in central nervous system (Salloum-Asfar et al., 2021;Qin et al., 2022). Preclinical and clinical models of ASD have demonstrated up-and down-regulations of several miRNAs. For instance, upregulation of miRNA-29b, miRNA329, miRNA-199b, miRNA-382, miRNA-296, miRNA-221 and miRNA-92 and downregulation of miRNA-146a, miRNA-146b, miRNA-130, miRNA-122a, miRNA-342 and miRNA-409 were described. Among these targets, for example, miRNA-199 is known to be able to regulate BDNF (Brain Derived Neurotrophic Factor), while miRNA-146a and miRNA-146b have been shown to be involved in the aetiopathogenesis of ASD (Urdinguio et al., 2010). Ivanov et al. (2022), developed a highly sensitive "silicon-on-insulator"-based nanosensor (SOI-NS) for the detection of ASD-associated miRNAs with a lower limit of 10 −17 M.
Nowadays, among all body sampling tissues, saliva is gaining a pivotal role as it is extremely easy and non-invasive to collect and monitor. It should be noted that when choosing a body fluid for biomarker measurement, it is important to avoid stress as much as possible, as stressful procedures can alter the production of hormones such as cortisol that influence the secretion of other hormones through feedsidewards interaction. Urine or blood collection is more stressful than saliva collection and this makes saliva an appropriate fluid for biological sampling in either patients or healthy subjects. Furthermore, its informative value is considered particularly high, as the salivary proteome is estimated to include more than 3,000 expressed proteins and peptides, comprising biomarkers of neurological diseases (Hicks et al., 2018b(Hicks et al., , 2020. After the first reports by Hicks and colleagues, several large-scale screening studies have focused on salivary poly-omic RNA measurements to identify children with ASD (Hicks et al., 2018a;Sehovic et al., 2020). The rise of the ASD epidemic and numerous hypothetical pathways to be elucidated both diagnostically and pathogenetically have given impetus to related research, resulting in an enrichment of recent evidence.
The complexity of ASD makes clinical diagnosis difficult, therefore, by identifying the biomarkers associated with ASD severity and combining them with the diagnosis, it is possible to better divide the factions within the spectrum and devise more targeted therapeutic strategies. To date, there are no reliable biomarkers to diagnose ASD or define its severity. This study aims to comprehensively examine miRNA expression in children with ASD and to reveal potential biomarkers and possible disease mechanisms so that they can be used to improve faction between individuals by promoting more personalized therapeutic approaches.

Samples collection
This prospective comparative and study included 10 children aged between 3 and 7 years. The ASD group comprised 4 males and 1 female (mean age 5.6 years) and the non-ASD group, healthy controls, comprised 3 males and 2 females (mean age 5.4 years) ( Table 1). Date on birth weight and parental age were retrieved for all participants and reported on Table 2. For approval, written informed consent was obtained from the parents/caregivers of each participant for approval. Exclusion criteria for the ASD participants included global developmental delay (cognitive disability with IQ < 70), g-tube dependence, confirmed chromosomal deletion or duplication, epilepsy, psychiatric diagnosis and premature birth of more than 6 weeks. Exclusion criteria for the non-ASD groups included a family history of ASD in a first degree relative or a chronic medical condition requiring daily medical care by a specialist. Common exclusion criteria among groups included acute illness, sensory disorders such as hearing or visual impairments, upper respiratory tract and oral cavity infections and wards of the state. Saliva samples were collected from all children in a non-fasting state using Oracollect RNA swabs (DNA Genotek; Ottawa, ON, Canada) through a simple, non-invasive process that consists of applying a highly absorbent swab for a few seconds under the base of the tongue. Each sample collected at least 500 µl of saliva. Samples were immediately sent to the laboratory for analysis. The clinical investigation protocol for this study was submitted to the ethical committee of the Policlinico Fondazione Cà Granda Milano Review Board, and approval has been successfully obtained (727_2021). miRNAs extraction miRNAs were extracted from saliva samples using NucleoSpin R miRNA Plasma kit (MACHEREY-NAGEL,  GmbH & Co. KG, Germany) and quantified by Agilent 2100 Bioanalyzer RNA assay (Agilent Technologies, Santa Clara, CA, USA). The NucleoSpin miRNA Plasma kit isolates total RNA including miRNA from saliva without the need to resort to the cumbersome phenol/chloroform extraction of a timeconsuming proteinase digest. The sample material is denatured in Lysis Buffer MLP. The protein is then precipitated by Protein Precipitation Buffer MPP and pelleted by centrifugation. After the removal of protein, the binding conditions for nucleic acids are adjusted by adding isopropanol. Total nucleic acids are bound to the NucleoSpin miRNA Column. The remaining nucleic acids are washed and eluted with minimal amounts of RNase-free water.

Library preparation and sequencing
QIAseq miRNA library kit (QIAGEN, Hilden, Germany) has been used for library preparation following the manufacturer's instructions.
RNA samples were quantified and quality tested by Agilent 2100 Bioanalyzer RNA assay. In addition, final libraries were checked with both Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and Agilent Bioanalyzer DNA assay. The libraries were sequenced on single-end 150 bp mode on NovaSeq 6000 (Illumina, San Diego, CA, USA).

Bioinformatic analysis
The bioinformatics analysis of small RNA-Seq includes: • Base calling and demultiplexing. Processing of raw data for format con-version and demultiplexing by Bcl2Fastq1 2.20 version of the Illumina pipe-line.
• The sRNAbench tool, a sRNAtoolbox2,3 tool, gives expression profiling of small RNAs, prediction of novel microRNAs, analysis of isomiRas, genome mapping, and read length statistics.

Pathways analysis and gene ontology
To better understand the functional role of the 41 selected miRNAs, a pathway prediction analysis was performed. For this purpose, the bioinformatics tool DIANA-mirPath (version 3) (Vlachos et al., 2015) and miRDB (Chen and Wang, 2020) were used. With this computational approach all target genes and related molecular pathways altered by the selected miRNAs were identified by using this computational approach. Furthermore, GO enrich-ment analysis was performed using the tool GO PANTHER version 14.0. 1 GO PANTHER analysis was conducted for the lists of genes obtained from DIANA-mirPath analyses and using the NPinter tool to evaluate the interactors of the top 10 dysregulated miRNAs.

Statistical analyses
Fold change values of miRNA expression levels were calculated through differential analysis. Student's t-tests were performed to select the differentially expressed miRNAs with a statistical significance, as reported in Tables 3, 4. For the Kaplan-Meier survival analysis, the non-parametric log-rank test was used to compare the survival distributions of the patients with ASD according to the down-regulation or overexpression of selected miRNAs. Data with a P-adj ≤ 0.05 were considered to indicative of statistically significant differences.

Identification of the miRNAs targeted genes and correlation analysis
For the pathway prediction analysis, all the 41 miRNAs were entered into the search bar of DIANA mirPath. The analysis revealed that modulated pathways and target genes existed for these miRNAs according to TarBase version 7.0 the reference database of DIANA mirPath. By performing a differential analysis of all the molecular pathways altered by miRNAs upand downregulated in ASD, it was possible to establish that all reported miRNAs could modulate several biological and pathological pathways, as reported in Tables 4, 5. In detail, up and downregulated miRNAs in ASD can alter 53 and 50 different pathways, respectively. Among the altered pathways, 11 were directly related to ASD. Among these 11 pathways, the most affected pathways were the Wnt signaling pathway and Axon Guidance (20 miRNAs up-regulated and 16 downregulated), PI3K-Akt signaling pathway (20 up-and15 downregulated), GABAergic synapse and TGF-beta signaling pathway (19 up-and 14 downregulated), Jak-STAT signaling pathway (19 miRNAs downregulated), circadian entrainment (13 miRNAs downregulated) and circadian rhythm (16 upregulated), Prion diseases (9 upregulated and 9 downregulated), Hedgehog signaling pathway (only 18 upregulated), Neurotrophin signaling pathway (only 20 miRNAs upregulated). Moreover, Hippo pathway, which is regulated by all 21 upregulated miRNAs and 16 downregulated miRNAs, is involved in neurodegenerative diseases. Therefore, it is evident that a high number of selected miRNAs were involved in the development of ASD and neurodegeneration (miRTarBase) (Tables 4, 5). Subsequently, all genes related to each miRNA were merged. In total, 12 genes were identified in common among all the 41 selected miRNAs were identified using this approach (mirPATHv.3). As shown in Table 6, the target genes of the 41 dysregulated miRNAs are PSAT1, STAT1, PIM1, MYO6, USP48, TWSG1, IL6, RBBP9, MMP14, ARL15, TRPC4, BRAF ( Table 6). Using the same criteria, Table 7 lists only the top 10 dysregulated miRNAs; target genes are ZNF106, GFPT2, JUNB, VPS33A. The subsequent analysis evaluated the interactors of the miRNAs. Twenty-three interactors were revealed by the analysis of the top 10 dysregulated miRNAs ( Table 8) (NPinter). Involved in ASD were: AASDHPPT (hsa-miRNA-1246; hsa-miRNA-199b-5p; hsa-miRNA-199a-3p; hsa-miRNA-199b-3p); AUTS2 (hsa-miRNA-199b-5p; hsa-miRNA-199a-3p; hsa-miRNA-199b-3p). IRDF1 is a predicted gene in ASD (hsa-miRNA-199b-5p; hsa-miRNA-199a-3p; hsa-miRNA-199b-3p) ( Table 8). The last step of the analysis consisted of the GO enrichment analysis by PANTHER, listing the 15 genes. As shown in Figure 2, the selected genes were grouped according to molecular function, biological process, and protein class (Figure 2). As regards the molecular function, it was observed that the majority of genes were involved in catalytic activity (GO:0003824) functions (46% of genes) (Figure 2A). When considering the biological processes, 68.8-43.8% of the genes were involved in cellular processes (GO:0009987) and metabolic processes (GO:0008152) ( Figure 2B). The most represented protein classes were protein modifying enzyme (PC00260), transporter (PC00227), and metabolite interconversion enzyme ( Figure 2C). Finally, pathway analysis revealed the association to molecular mechanisms for only a few genes. Indeed, only 4 genes out of the 15 recognized were assigned to a molecular pathway ( Figure 2D). The most represented pathways were the Inflammation mediated by chemokine and cytokine signaling pathway (P00031), Interleukin signaling pathway (P00036), PDGF signaling pathway (P00047), Alzheimer's diseasepresenilin pathway (P00004), and EGF receptor signaling pathway (P00018).

Discussion
Autism spectrum disorders is a highly heritable neurodevelopmental disease with a sharply increasing incidence (Santini et al., 2014). The diagnosis of ASD is Frontiers in Neuroscience 05 frontiersin.org   Prion diseases (hsa05020) 5.28E + 04 11 9 Fatty acid biosynthesis (hsa00061) 1.28E-01 6 6 Wnt pathway, axon guidance and GABAergic synapse pathways adversely affect neurodevelopment and lead to the pathogenesis of ASD (Giudice et al., 2017;Somekh, 2021). The data obtained show that many miRNAs are involved in the TGF-beta signaling pathway Dana et al., 2020), involved in ASD development, justifying the presence of dysregulated miRNAs in these other pathways.
The present findings show that the adherents junctions, FoxO, Hippo and the Hedgehog signaling pathway are altered.
Interestingly, miRNA-199b-5p increased significantly in saliva of children with ASD, which is involved in the regulation of the AUTS2 gene («miRDB» database). The miRNA-199a-3p, which was observed to be strongly increased in the present work, is considered to be involved in the regulation, among all the other targets, of the CB1 receptor (Cannabinoid receptor type 1) («miRDB» database). Indeed, CB1 receptor, responsible for the social reward, has been observed reduced in the brain of people with ASD postmortem (Shohami et al., 2011). The miRNA-4516 was also found upregulated in the serum of ASD patients (Huang et al., 2015), which is considered to be involved in the sodium voltage-gated channel («miRDB» database). Altered gene expression of this channel is considered to be a solid predictor of ASD (Yamakawa, 2016). As previously Bacterial invasion of epithelial cells (hsa05100) 0.04960 32 12 Prion diseases (hsa05020) 1.54E + 02 14 9 Dorso-ventral axis formation (hsa04320) 0.045132 15 9 Mucin type O-Glycan biosynthesis (hsa00512) 9.65E + 04 14 8 Fatty acid biosynthesis (hsa00061) 2.11E-01 6 4 Glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate (hsa00532) 0.03317 7 3 Frontiers in Neuroscience 09 frontiersin.org observed in the blood of children with ASD (Huang et al., 2015), an increase of miRNA-1246 in saliva was also found. miRNA-1246 is predicted to be implicated in the regulation of SYN II (Synapsin II) («miRDB» database), which was identified as an ASD predisposition gene (Corradi et al., 2014). Also, miRNA-1246 together with the miRNA-4516, is involved in circadian cycles, and metabolic processes («miRDB» database). Furthermore, miRNA-1246 is predicted to regulate the expression of the voltage-gated sodium channel («miRDB» database). In this sense, the discovery of an altered miRNA-1246 in the saliva, which follows the same trend as observed in the peripheral blood, suggests that it could be a suitable peripheral, non-invasive biomarker for the early diagnosis of ASD. Furthermore, a decrease in miRNA-454-5p, which is among the other functions is involved in CREB (cAMP Response Element-Binding) regulation («miRDB» database) was observed. miRNA-660-5p was also found decreased in saliva of children with ASD. Interestingly, this miRNA is implicated in regulating insulin receptor substrate 1 («miRDB» database). miR-4284 was found to be decreased in patients with ASD. This miRNA is in fact implicated, as is miR-1246, in the regulation of SYN II gene («miRDB» database). Miao et al. (2021), identified miR-4284 as an anticancer in colon cancer. It inhibits colon cancer tumorigenesis by reducing PLIN5 and inhibiting EMT. Thus, miR-4284 could be a potential therapeutic target in metastatic colon cancer. Another miRNA that resulted downregulated is miRNA-203a-5p. This miRNA is predicted to be implicated in the regulation of ERK1 («miRDB» database) Cossu et al., 2019). Lastly, a decrease of the miR-1973 was observed. This miRNA is associated, among other functions, with the regulation of the potassium gated voltage channel.
Furthermore, analyses performed on the genes regulated by the first 5 up-and downregulated miRNAs showed the main pathways related to these miRNAs. Among other pathways, it was found that the EGF (epidermal growth factor) receptor signaling was affected by two miRNAs. Consistent with the present findings, the EGFR signaling was observed to be altered in children with ASD (Vallés and Barrantes, 2021). Another interesting pathway observed is that of the nicotinic receptor's pathway, which resulted regulated by three different miRNAs. Indeed this receptor was observed to be altered in people with ASD, and more in particular in the frontal and parietal cortex, as well as in the cerebellum, which is deeply implied in social cognition, which means that cerebellar alteration could be involved in eye avoidance and similar behaviors. PDGF (Platelet-Derived Growth Factor), was the single gene most regulated by the miRNA studied in the present work, with a total of four different miRNAs. Although there are only a few papers are available that demonstrate a connection between PDGF and ASD, it should be noted that PDGF was indeed observed increased in the plasma of children with ASD (Zakareia et al., 2012). Another noteworthy element that was found regulated by a miRNA from the studied pool, was the JAK/STAT signaling. The tumor suppressor p53 was also observed regulated by two different miRNAs of the first five analyzed. In fact, in a preclinical model of ASD, it was suggested that dysregulation of p53 was correlated with altered DNA repair capacity and genomic instability that might lead to ASD (Wong et al., 2016). Some of the miRNAs in our study are found to be associated with inflammatory and proliferative diseases ( Table 9). In particular, hsa-miRNA-1246 has been reported as a biomarker secreted in different biological samples in many neoplastic diseases. In fact, studies in colorectal cancer cell lines have shown oncogenic role of miR-1246. In this cancer, the methyltransferase METTL3 oncogene has been shown to increase methylation of pri-miR-1246 to enhance maturation of pri-miR-1246. Notably, miR-1246 has been predicted to suppress expression of the Sprouty Related EVH1 Domain Containing 2 (SPRED2) tumor suppressor, by increasing MAPK pathway (Peng et al., 2019). In this study, we would like to report that this miRNA is upregulated with highest fold change.
Furthermore, our analysis showed that genes in significant pathways that are predicted to be targeted by deregulated miRNAs control apoptosis,  Frontiers in Neuroscience 11 frontiersin.org cell death of immune cells, cell cycle progression, epilepsy and neuronal cell proliferation. Our results reinforce the findings of other groups regarding important pathways, such as the regulation of cell growth, the role of organs besides the brain in ASD, and overlap with other non-neurodevelopmental diseases, such as cancer.

Conclusion
In our pilot study, we screened differentially expressed miRNAs in the saliva sample of ASD children and healthy control. Our sequencing investigation identified the significant number of miRNAs that were up-or downregulated in the ASD group. The targets of miRNAs under study have been observed through specific tools to predict associated pathways and its regulation. The dysregulated miRNAs and their target genes were found to be involved in developmental pathways already known as the cause of typical ASD and neurobiological dysfunctions. This pilot study can be improved on biostatistical grounds by increasing the size of ASD cohort. Also, the lack of proteomic and genetic validation by RT-PCR of the specific miRNAs, and target genes, remains an unfinished task to claim these miRNAs as biomarkers for ASD.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement
The study was conducted in accordance with the Declaration of Helsinki and approved by Ethics Committee "DIAGNOSI PRECOCE DEI DISTURBI DELLO SVILUPPO NEUROLOGICO IN SALIVA, INDIVIDUAZIONE DI MARCATORI PREDITTIVI ATTRAVERSO BIOSENSORI" Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico, Milano, U.O.C. Chirurgia Maxillo-Facciale e Odontostomatologia (727_2021). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author contributions
ZK, MM, AS, and GTo participated in the conceptualization and writing-original draft preparation. MM and GTo participated in the methodology and software. AA participated in the data curation. MC and TB participated in the formal analysis. AC and FI participated in the resources. MB, MD, and GTa participated in the validation, supervision, and writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Funding
This study was funded by internal funds of the University of Milan, Fondazione IRCCS Cà Granda, Ospedale Maggiore Policlinico, Milan, Italy.