Analysis of Faecal Microbiota and Small ncRNAs in Autism: Detection of miRNAs and piRNAs with Possible Implications in Host–Gut Microbiota Cross-Talk

Intestinal microorganisms impact health by maintaining gut homeostasis and shaping the host immunity, while gut dysbiosis associates with many conditions, including autism, a complex neurodevelopmental disorder with multifactorial aetiology. In autism, gut dysbiosis correlates with symptom severity and is characterised by a reduced bacterial variability and a diminished beneficial commensal relationship. Microbiota can influence the expression of host microRNAs that, in turn, regulate the growth of intestinal bacteria by means of bidirectional host-gut microbiota cross-talk. We investigated possible interactions among intestinal microbes and between them and host transcriptional modulators in autism. To this purpose, we analysed, by “omics” technologies, faecal microbiome, mycobiome, and small non-coding-RNAs (particularly miRNAs and piRNAs) of children with autism and neurotypical development. Patients displayed gut dysbiosis related to a reduction of healthy gut micro- and mycobiota as well as up-regulated transcriptional modulators. The targets of dysregulated non-coding-RNAs are involved in intestinal permeability, inflammation, and autism. Furthermore, microbial families, underrepresented in patients, participate in the production of human essential metabolites negatively influencing the health condition. Here, we propose a novel approach to analyse faeces as a whole, and for the first time, we detected miRNAs and piRNAs in faecal samples of patients with autism.


Introduction
Autism spectrum disorders (ASD) refers to a group of complex neurodevelopmental conditions whose core symptoms are a deficit in communication and social interaction, restricted interests, and repetitive behaviours (https://www.who.int/news-room/factsheets/detail/autism-spectrum-disorders, accessed on 26 November 2021). Comorbidities, such as mental retardation, epilepsy, anxiety, sensory, sleep, and gastrointestinal disorders, as well as food selectivity often occur in ASD [1][2][3]. ASD manifests during the first years of age, and in the last decades, its prevalence has continued to increase, reaching the frequency of one in 54 children aged eight years in the USA in the 2016 [4]. ASD aetiology is multifactorial and has not yet completely elucidated although the interaction between genetic susceptibility and environmental factors is emerging as the most consistent cause of ASD development and severity [2,5,6]. To date, hundreds of risk genes have been associated with ASD as reported in Simons Foundation Autism Research Initiative (SFARI) database Interestingly, these miRNAs can enter bacteria and act at the DNA or RNA level to affect gene expression and control the microbial growth [34]. By means of this mechanism, the host can regulate the composition of its own microbiota that, in return, can influence host miRNA expression via MyD88-dependent pathway [11,35,36]. Furthermore, microbiotaderived metabolites and microbiota-derived EVs also participate in the host-microbiota cross-talk regulating gene expression and intestinal homeostasis of the hosts [11,37,38]. The potential cross-talk between faecal microbiota and miRNA expression in pathological conditions has been identified in inflammatory bowel disease and colorectal cancer, also underlying their potential clinical relevance as biomarkers and therapeutic targets [35,39]. For instance, miR-515-5p and miR-1226-5p induce the growth of Fusobacterium nucleatum and E. coli, respectively [40], while commensal microbiota-induced miR-21-5p over-expression is involved in intestinal permeability via ARF4 [40], thus representing a therapeutic target to restore an intestinal barrier.
Recently, a study found an association between salivary miRNAs and salivary microbiota dysregulation in ASD [41], but none of the studies reported host miRNA-gut microbiota interaction in this condition. PiRNAs are small ncRNA that epigenetically and post-transcriptionally silence the expression of transposable elements integrated in eukaryotic genomes [42]. At present, no piRNAs in stool samples are available either in healthy or in pathological conditions.
In this scenario, it emerges that microbes and small non-coding RNA (sncRNA) in faecal samples should be considered and studied as a whole to comprehend how microbial strains interact among each other and with the host. In this pilot study, we defined, through "omics" technologies, the faecal micro-and mycobioma profile as well as the sncRNA profile of a small group of individuals with ASD and neurotypical controls. Our aim was to find markers of ASD among stool microbial and transcriptional modulators for the possible relationship between them and the host. We applied a bioinformatics approach to correlate gut bacteria and fungi composition with host miRNAs and piRNA expression for the first time in stool samples from patients with ASD and attempted to highlight possible mechanisms of microbiota-host bidirectional cross-talk in ASD.

Ethical Committees
The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of Istituto San Vincenzo (protocol code 312, 28 December 2018). The Ethics Committee prepared the informed consent, including the instructions for sample collection; and enrolled the children. The collected informed consents were signed by the parents and/or legal guardians since the individuals were underage. All methods were performed in accordance with relevant guidelines and regulations regarding observational studies.

Sample Collection
Naturally evacuated stool samples were obtained from all individuals and collected by previously instructed parents. Stools were collected in stool nucleic acid collection and transport tubes, then returned refrigerated to the Institute for Biomedical Technologies, CNR. Samples were stored at -80 • C until RNA/DNA extraction.

DNA and SmallRNA Extraction from Stool
Total DNA was extracted from frozen stool samples (200 mg) using commercial kit and the relative protocol for pathogen detection (QiAamp DNA stool mini kit, Qiagen GmbH, Hilden, Germany) with minor modification, and then, pre-lysis mechanical grinding was performed to increase sample homogenization, microbial lysis, and DNA extraction strength. Zirconia beads (100 µm in diameter) were added to ASL buffer (300 mg/mL buffer) before incubation at 95 • C, and three cycles per 1 min in a bead beater were performed before thermal lysis. The DNA quality (280/260 ratio) was checked by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and quantity measured by Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Wilmington, DE, USA).
RNA was isolated from frozen stool samples (200 mg) using the commercial kit RNeasy Power Microbiome (Qiagen GmbH, Hilden, Germany) according to the manufacturer's protocol. RNA quality was assessed by Agilent RNA 6000 Nano on Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA), and RNA concentration measured by Qubit RNA assay (Thermo Fisher Scientific, Wilmington, DE, USA).

SmallRNA Sequencing
Small RNA-sequencing libraries were generated directly from total RNA isolated from the stool samples and performed with the TruSeq Small RNA Library Preparation Kits (Illumina, Inc., San Diego, CA, USA) based on the manufacturer's protocol. Libraries were sequenced on NextSeq 500 (Illumina, Inc., San Diego, CA, USA), 1 × 75 bp and 30 million reads per sample.

Metataxonomic Bioinformatics Analysis
Metataxonomic analysis was performed in R (4.0.3) using Dada2 [44] pipeline using phyloseq package [45] against DADA2-formatted reference databases latest available version (July 2021), that is, Silva v138 [46] for 16S and Silva v132 [47] for 18S. DADA2 plug-in was used to filter, trim, dereplicate, merge, remove chimaeras, and assign taxonomy to all produced sequences to obtain the Operational Taxonomic Units (OTUs). Variance Stabilising Transformation was applied to normalise across samples on OTUs with DESeq2 package [48] as described by McMurdie and Holmes [45]. For each sample, the number of observed OTUs and the percentages of relative abundances of phyla, orders, classes, and families were determined. To evaluate statistically significant differences between ASD and controls (Ctrl) at genus level, the univariate DESeq2 method was used [45,48]. Default Wald test was applied in DESeq2, and significance threshold was set to p-value < 0.05 for the 16S analysis, while for the 18S analysis, all results were considered.

Small, Non-Coding RNA Data Analysis
SmallRNA reads were processed according to a custom bioinformatics pipeline that we developed [49]. Summarising the main pipeline steps, smallRNA reads were checked for quality control using FastQC package (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc, accessed on 22 July 2020), filtered, and then mapped against Arena-Idb [50], a reference database representing a comprehensive and non-redundant dataset of public ncRNA sequences and annotations, using Bowtie aligner [51], with one mismatch in the leftmost 20 bp of the read. In order to obtain reliable read counts and to fix the problem of multireads [52], we used the RSEM tool [53] for accurate expression estimations of identified ncRNAs.
An evident heterogeneity in the expressions of several individual references required an accurate management of the expression normalization step. A reference-free clustering of the sequences was performed with the SEED [54], an algorithm for clustering very large NGS sets. Sequences were joined into clusters that differ by up to three mismatches and three overhanging residues from their virtual centre. The cardinalities of the clusters resulted to be more stable (by showing higher correlations among the samples) and were used to compute the scaling factors of TMM normalization. Such factors were applied to normalise the ncRNA reference expression counts. Expression data were analysed with edgeR (https://bioconductor.org/packages/release/bioc/html/edgeR.html, accessed on 22 July 2020). EdgeR package applied Robinson and Smyth exact statistical methods for multigroup experiments [55,56]. The Benjamini-Hochberg multiplicity correction method was used on the p-values to control the false-discovery rate (FDR).
Analysis of 16S was performed considering all ASD samples against all Ctrl samples. As displayed in Figure 1a, no relevant differences can be evidenced between the two groups from the alpha diversity analysis. Instead, alpha diversity analysis, in particular Shannon and Simpson indices, of 18S samples ( Figure 1b) highlight a major species diversity for ASD samples compared to controls.
Metataxonomic comparative analysis was performed at family/genus level. ASD samples displayed a reduced microbiota variability compared to controls. In particular, 24 families, corresponding to 38 bacteria genera, were significantly detected in all samples; among these, nine families are prevalent in healthy samples and five in ASD ( Figure 2). Considering differences in microbiota composition at genus level, ASD samples displayed a prevalence of Bifidobacterium, Desulfovibrio, Coprococcus, Alistipes, and Sutterella. Instead, in neurotypical samples, the most commonly represented genera are Phascolarctobacterium, Akkermansia, Barnesiella, Enterorhabdus, Lachnospiraceae_NK4A136, Ruminococcus, Prevotellaceae_UCG-001, and Streptococcus.
Metataxonomic comparative analysis was performed at family/genus level. ASD samples displayed a reduced microbiota variability compared to controls. In particular, 24 families, corresponding to 38 bacteria genera, were significantly detected in all samples; among these, nine families are prevalent in healthy samples and five in ASD ( Figure 2
Metataxonomic comparative analysis was performed at family/genus level. ASD samples displayed a reduced microbiota variability compared to controls. In particular, 24 families, corresponding to 38 bacteria genera, were significantly detected in all samples; among these, nine families are prevalent in healthy samples and five in ASD ( Figure 2  Due to high heterogeneity of the sample compositions (see Supplementary Figure S1), we performed the metataxonomic analysis comparing each ASD sample to the whole control set (n = 6) and thus evidenced that twenty bacteria families were prevalent in the most Ctrl samples but were not so frequent in ASD ( Figure 3, Supplementary Figure S2 and Table S1). These included the previously identified Akkermansiace, Barnesiellaceae, Eggertellaceae, and Tannerellaceae characterising Ctrl samples, while Sutterellaceae for the ASD samples and Acidaminococaceae and Ruminococcaceae were detectable in both sample types. Instead, several families displayed a different trend: Bacterioidaceae, Bifidobacteriaceae, Christensellaceae, Erysipelotrichaceae Lachnospiraceae, and Streptococcaceae were shown as characteristic of Ctrls, while Pasturellaceae and Prevotellaceae were detectable only in ASD. Moreover, Coriobacteriaceae and Veillonellaceae were not previously identified as well as Lactobacillaceae and Synergisticaceae, which were respectively under-and over-represented in ASD samples.
Among these families, Prevotella_9 and Ruminococcus genera are mainly present in ASD samples, while Phascolarctobacterium, Akkermansia, Bacteroides, CAG-352, and Dialister genera were principally detectable in Ctrl samples.

Fungi Profiling: Metataxonomic Analysis
Although the 18S analysis was performed on the Fungi kingdom, in order to reduce the huge amount of vegetable sequence detected in these samples, low abundances can be detected in all samples. Moreover, the mycobiome comparative analysis returned even less variability in ASD samples. Only five fungi families can be significantly identified: Saccharomycetaceae and Debariomycetaceae in both sample types, while Aspergillaceae, Malassenziaceae, and Cladosporaceae are detectable only in controls ( Figure 2b). Although Saccharomycetaceae was the most diffused family, it was not further classifiable at genus level. Among ASD, the Debaryomycetace family was predominant and represented by Candida-Lodderomyces_clade and Meyerozyma-Candida_clade. In neurotypical individuals, the predominant genera were Penicillium and Malassezia, belonging to Aspergillaceae and Malasseziaceae, respectively.
The analysis performed comparing the mycobiome profile of each ASD sample to those of the whole Ctrl group (Supplementary Table S2) returned Aspergillaceae and Malasseziaceae as prevalent in control samples, while Debaryomycetace and Dipodascaceae were predominant in ASD samples ( Figure 3). This analysis confirms for the Debaryomycetace family Candida-Lodderomyces_clade and Meyerozyma-Candida_clade genera and annexes Geotrichum from Dipodascaceae characterising ASD mycobiota and Penicillium, belonging to the Aspergillaceae family, detectable only in Ctrl samples, while Malassezia were seriously reduced in ASD samples.

sncRNA Profiling
The analysis of human sncRNAs in stool samples was performed by the bioinformatics pipeline as described in [49]. Quality control trimming and filtering returned sample reads in the average of 26.2 million and 15-51 base length. About 2.4 million reads per sample (9.2%) mapped to the human genome (GRCh38), and an average of 1.2 million reads per sample were assigned to ncRNA classes. Among these reads, 58.2% represented long ncRNAs, and 41.8% sncRNAs. The 0.5% fraction of ncRNA reads were miRNAs, according to previously published literature [61], and 9.6% piRNAs.
We identified a total of 11,596 sncRNAs in stool samples, and no relevant differences were observed between cases and controls; a mean of 1271 piRNAs (53.4%) and 444 miRNAs (18.6%) per sample were detected. The distribution of the different classes of sncRNA did not differ between samples of ASD patients and controls ( Figure 4). The most expressed miRNAs in both patients and controls are hsa-miR-182-5p and hsa-miR-681; hsa-miR-657 and hsa-miR-2110 are mainly present in all ASD, while hsa-miR-1203 is mainly present in neurotypical individuals. The most represented piRNAs in all samples are hsa-piR-27489, hsa-piR-32912, hsa-piR-32921, hsa-piR-23722, and hsa-piR-19705; in addition, piRNAs hsa-piR-28059 is highly expressed in ASD subjects, while hsa-piR-33182 and hsa-piR-33031 are highly expressed in healthy subjects.  As for the metataxonomic analysis, the sncRNA analysis was also performed comparing each ASD sample to the whole controls collection (n = 6) (Supplementary Table S3). A total of 42 miRNAs were up-regulated in three out of six ASD samples although none resulted common to all the three samples. Target gene analysis of the dysregulated miR-NAs identified 18 target genes common to three ASD stool samples ( Table 1). The differential expression analysis conducted on piRNAs returned no down-regulated piRNAs  Table 1. List of miRNAs, piRNAs and miRNA + piRNAs common targets. "Single-cell type specificity" was obtained by single-cell transcriptomics; "Tissue specificity (RNA)" and "GI RNA expression" were obtained by RNA-seq; "GI protein expression" was obtained by immunocytochemistry investigations. All these data are from ProteinAtlas (https://www.proteinatlas.org, accessed on 2 December 2021). Biological process and Molecular function data are from UniProt (https://www.uniprot.org, accessed on 2 December 2021). SFARI score categories are reported in the following website (https://gene.sfari.org/about-gene-scoring/, accessed on 2 December 2021). Abbreviation: GI, gastrointestinal; ASD, autism spectrum disorders.  As for the metataxonomic analysis, the sncRNA analysis was also performed comparing each ASD sample to the whole controls collection (n = 6) (Supplementary Table S3). A total of 42 miRNAs were up-regulated in three out of six ASD samples although none resulted common to all the three samples. Target gene analysis of the dysregulated miRNAs identified 18 target genes common to three ASD stool samples ( Table 1). The differential expression analysis conducted on piRNAs returned no down-regulated piRNAs and a total of 84 up-regulated piRNAs in four out of six patients with ASD. Among these, only hsa-piR-21363 was common to two patients, while three target genes (CFLAR, GOLGA6L2, and SLC2A4) were common to four patients. Moreover, considering both miRNA and piRNA gene targets, five genes that resulted common in more than two samples (N4BP1, SLC2A4, SLC12A6, TTN, and ZNF33A) were identified. Table 1 summarises the 26 target genes of identified transcriptional modulators, considering the tissue expression, protein annotation, related disease, and SFARI score for proteins mutated in ASD.
The miRNA and piRNA target genes common to at least three samples were annotated with KEGG database and MSigDB-Hallmark gene sets. Results are displayed in Figure 5, grouped by pathway (a) or Hallmark gene set (b). This analysis returns pathways or biological processes involving one or more target genes. These are related to cell-cell junction, bacterial invasion, inflammation, and metabolite signalling and are pathways linked to autism.

Case Study: Analysis of Siblings
Within the analysed stool samples, four were obtained from two couples of siblings: couple #1 (male-male), one subject with ASD and the other neurotypical control, and couple #2, a male with ASD and a neurotypical female. We separately analysed the data from these samples to highlight differences of gut microbial community and sncRNAs between patient and control with common genetic background and similar diet.
The microbiota composition of the two couples of siblings was compared; abundance fractions are reported in Figure 6. Couple #1 displays a similar microbiota composition, and the only relevant differences concern the Bacteroidaceae and the Rumicococcaceae family (Figure 6a), which, respectively, increased and reduced in the ASD individual, in contrast to Debaryomycetaceae (Figure 6c), which was noticeable only in the ASD sibling mycobiota to the detriment of Saccharomycetaceae. A different composition was identified in the couple #2, with the ASD sibling displaying increased Akkermansiaceae, Bacetroidaceae, Lachnospiraceae, and Muriobaculaceae, while Rikenellaceae and Veillonellaceae were decreased in the ASD individual (Figure 6b). The mycobiota composition displays a relevant increase in Debaryomycetaceae and Malasseziaceae in ASD to the detriment of Saccharomycetaceae (Figure 6d). By comparing the results obtained within each couple of siblings, we identified upregulated and down-regulated miRNAs and piRNAs by statistical analyses and defined as statistically significant those ncRNA with log2 fold change ≤ −1 or ≥ 1 and p-value < 0.05 (Fisher's test) (Supplementary Table S4). A total of 750 miRNAs were identified in couple #1, and among these, 93 have a significantly differential expression. Among these, 60 miRNAs were detected in both individuals of couple #1, 32 were up-regulated, and 28 down-regulated in ASD sample. Moreover, 15 miRNAs were detected only in the ASD subject and 18 only in the healthy sibling. The top three down-regulated miRNAs in ASD sample are hsa-mir-937, hsa-mir-3197, and hsa-mir-103a-1, while those up-regulated are hsa-mir-4700, hsa-miR-657, and hsa-mir-2110. A total of 2177 piRNAs were identified in stool samples from couple #1, and 329 were significantly dysregulated. Among these, 112 are up-regulated and 52 down-regulated and present in both subjects belonging to couple #1, whereas 109 are present only in the stool sample from the ASD subject and 56 in sample from the healthy sibling (for details, see Table 2). The top three down-regulated piRNAs in ASD are hsa-piR-6691, hsa-piR-6693, and hsa-piR-29205, while hsa-piR-28269, hsa-piR-32987, and hsa-piR-28059 are significantly up-regulated.   By comparing the results obtained within each couple of siblings, we identified upregulated and down-regulated miRNAs and piRNAs by statistical analyses and defined as statistically significant those ncRNA with log2 fold change ≤ −1 or ≥ 1 and p-value < 0.05 (Fisher's test) (Supplementary Table S4). A total of 750 miRNAs were identified in couple #1, and among these, 93 have a significantly differential expression. Among these, 60 miRNAs were detected in both individuals of couple #1, 32 were up-regulated, and 28 downregulated in ASD sample. Moreover, 15 miRNAs were detected only in the ASD subject and 18 only in the healthy sibling. The top three down-regulated miRNAs in ASD sample are hsa-mir-937, hsa-mir-3197, and hsa-mir-103a-1, while those up-regulated are hsa-mir-4700, hsa-miR-657, and hsa-mir-2110. A total of 2177 piRNAs were identified in stool samples from couple #1, and 329 were significantly dysregulated. Among these, 112 are up-regulated and 52 down-regulated and present in both subjects belonging to couple #1, whereas 109 are present only in the stool sample from the ASD subject and 56 in sample from the healthy sibling (for details, see Table 2). The top three down-regulated piRNAs in ASD are hsa-piR-6691, hsa-piR-6693, and hsa-piR-29205, while hsa-piR-28269, hsa-piR-32987, and hsa-piR-28059 are significantly up-regulated.
In couple #2, we identified 581 miRNAs and 1912 piRNAs. There are 51 significantly differentially expressed miRNAs common to both samples, 11 up-regulated and 40 downregulated in the ASD sample. Moreover, 19 miRNAs were identified only in the ASD sample and 23 only in the control sibling. In couple #2, 222 piRNAs significantly differentially expressed and common to both samples were identified. Among these, 66 were up-and 256 down-regulated in ASD sample, 33 detected only in ASD, and 85 only in sample from the neurotypical sister. The comparison between couple #1 and #2 reveals that there are three common miR-NAs down-regulated in samples from ASD, namely hsa-miR-10b-5p, hsa-miR-22-3p, and hsa-miR-192-5p, and four up-regulated, namely hsa-miR-6760-5p, hsa-mir-6766, hsa-mir-6839, and hsa-mir-3976. The common dysregulated piRNAs are 44, down-regulated 13, and up-regulated 31 (for details see Table 3).

Discussion
The role of gut microbial composition and interactions among the microbial strains as well as the bidirectional communication between the host and gut microbiota are continuously emerging in health and disease [13]. Faeces contain much information useful to shed the light on gut microorganisms and the molecular trafficking involved in the hostmicrobiota cross-talk. For these reasons, faeces samples should be analysed as a whole, at the microbial and molecular level. Nevertheless, most studies about gut microbiota only refer to the prokaryotic component, and the role of eukaryotes and host-microbiota communication is poorly investigated. Thanks to the progress achieved in present years regarding the high-throughput sequencing technologies and in the computational analysis, due to the data generated by these technologies, we produced the eukaryotic and prokaryotic profile as well as the host transcriptional modulator profile (miRNAs and piRNAs) gathered in the faeces collected from a small group of children with ASD or neurotypical development.
Considering bacterial composition, we found a dysbiosis consisting of lower alpha diversity in ASD stools. This is in line with the literature that associates this index with the severity of social impairment, while it is independent of IQ [10,18]. We found that Firmicutes and Bacteroidetes were the most abundant phyla in both ASD and Ctrl faeces samples, and the Bacteroidetes/Firmicutes ratio was higher in ASD samples than in Ctrl, mainly due to decrease of Firmicutes. The imbalance in this ratio is still known in literature [72]; however, some evidence sustains that the Bacteroidetes/Firmicutes ratio was higher in ASD samples than in Ctrl [59,60], while other studies claim otherwise [33,73]. Due to this inconsistency, the role of Bacteroidetes/Firmicutes remains controversial.
At the genera taxonomic level, the Acidaminococcus enrichment and Dialister depletion that we observed in ASD is in line with previous articles [33,74]. Instead, Prevotella_9 and the overall increase in the Prevotella family is arguably reported in ASD children [74,75]. Moreover, an increase in Ruminococcus and Streptococcus and a decrease in Agathobacter were previously associated to ASD [59,76] and negatively associated with sleep and language disorders [77]. In line with our results, the presence of Alistipes in ASD, or more generally in patients with neurodevelopmental disorders, has also been discussed [31,33,78]; Alistipes probably has a role in decreasing serotonin availability, destabilising the gut-brain axis. Finally, literature reports Christensenellaceae as signature of a healthy gut [79,80]. Furthermore, Akkermansia is related to a healthy gut, and it is thought to have anti-obesity and anti-diabetic effects [81,82]; this family promotes gut barrier integrity, modulates immune response, and inhibits inflammation and syntrophy with other microbiota species [83]. Akkermansia is evidenced as mucolytic bacteria, and its lower abundance (or absence) in ASD supports the hypothesis of possible mucus barrier in children with autism [60,84]. Besides, short chain fatty acid (SCFA) produced by several bacteria families are involved in the release of mucus, with the exception for succinate not consumed by Phascolarctobacterium [85] in ASD cohort children. Bacteria families involved in SCFA production were prevalently identified in controls, namely Akkermansia [86] Phascolarctobacterium, Lachnospiraceae, and Agathobacter, a producer of butyrate and beneficial SCFA [77,[87][88][89]. However, recent evidence suggests an involvement of acetate and propionate in various disorders, including obesity [90]. Overall, SCFAs play a modulatory role in the microbiota-gut-brain axis, mediating behaviour and intestinal physiology [91]. In particular, butyrate, the levels of which are known to decreased in ASD [92], has multiple effects on the gut-brain axis, anti-inflammatory, blood-brain barrier, and gut-permeability regulators [93] and on healthy gut physiology, preventing pathogen invasion, modulating the immune system, and reducing cancer progression [94].
Furthermore, for the fungi kingdom, the mycobiome comparative analysis identified only Saccharomycetaceae and Debariomycetaceae families in both ASD and Ctrl groups, while Aspergillaceae (including Penicillium) and Malassenziaceae were found only in neurotypical according to literature [33,95]. The Malassezia role in the healthy gut is still unclear: it is correlated with neurological disease [96], but it represents a large amount of breast milk mycobiota, and it is found in healthy faecal samples [97]. In ASD, instead, only the Debaryomycetace family was predominant and consists of Candida-Lodderomyces_clade and Meyerozyma-Candida_clade that have never been reported before in ASD or in any other gut microbiota study. Although fungi are poorly studied, an important presence of Candida among ASD gut prokaryotes was reported, and anti-Candida albicans IgG antibodies were also detected in the plasma of children with ASD, making genus Candida a new microbial risk factor for the condition [33,98]. Candida, when overgrown, produces root-like structures that penetrate the intestinal wall, causing the leaky-gut syndrome. This allows macromolecules, such as toxins and food antigens, to enter the bloodstream and trigger food intolerance and allergies, as described in a subset of children with ASD [99].
We investigated sncRNA by next-generation sequencing technologies, with the advantage of looking for all sncRNA classes. Indeed, previous studies about sncRNAs in ASD were performed by quantitative real-time PCR or arrays technologies that investigate only a relatively small panel of target miRNAs. Therefore, these studies may return a focused profile, while a complete and unbiased sncRNA view, such as that obtained by sncRNA-seq analysis, is needed. Indeed, as emerged by our results, some sncRNAs that we found over-expressed in ASD target genes and pathways involved in inflammation and intestinal permeability, while others are involved in ASD. Due to the nature of investigated samples, we are unable to define if these latter sncRNAs are limited to stool and act only at gastrointestinal level or pass into the bloodstream and act at systemic level, too. The most represented miRNAs in all ASD are hsa-miR-2110 and hsa-miR-657. The first one was found in EV released from mesenchymal stem cells [100], up-regulated in serum exosome from patients with glioblastoma [101], supporting the hypothesis that hsa-miR-2110 is released into the lumen. The increased expression of hsa-miR-657 is associated to inflammatory response in gestational diabetes mellitus [102]. Interestingly, maternal gestational diabetes mellitus increases the risk factor for autism in offspring [103].
Host miRNAs are able to enter bacteria and regulate their gene expression and growth [34]. However, our results cannot be compared to those of the literature because the study is limited to a small number of miRNAs on Fusobacterium nucleatum and E. coli [34].
In our study, piRNAs are the most represented sncRNA. This is of particular interest because, although the role of piRNAs is still emerging, they have never been investigated so far in the faeces of individuals with ASD. PiRNAs data refer mainly to germline or are related to cancer cells [104,105] or faeces of patients with colorectal cancer [61]. They were detected in EVs in different body fluids [106,107], so we speculate that the piRNAs we found in stools may have been released into the gut lumen and may interact also with micro-and mycobiota.
Considering the miRNA and piRNA target genes common to multiple samples, four were already present in SFARI database: NACC1, SMAD4, TNRC6B, and TTN, and their mutants are implicated in ASD and other neurodevelopmental disorders (See Table 1). Besides, CNV in locus 9p24.3, overlapping CBWD1, was identified in autistic patients [108][109][110]. Numerous genes are also involved in mental retardation, schizophrenia, and other neurodevelopmental disorders (see Table 1 for details). We suppose that the up-regulated miRNAs and piRNAs can negatively influence the expression of these genes and therefore dysregulate the biological pathways involved in neurodevelopment. SLC2A4 is a glucose transporter whose impairment has again been associated with diabetes mellitus [111]. N4BP1 is a potent suppressor of cytokine production, and thus, its inhibition increases innate immune signalling and inflammation [112]. Finally, HSBP1 and IGF1R were found differentially expressed in autism [113][114][115][116][117] and, directly or indirectly, modulated by microbiota [118,119]; namely, HSBP1 is down-regulated after antibiotic administration, and IGF-1 level is affected by gut microbiota composition.
These 26 genes belong to pathways that are, directly or indirectly, associated with ASD [69] or its comorbidities. For instance, IGF1R, SMAD4, and WASF2 genes are involved in "Adherens junction" KEGG pathway and therefore could induce an imbalance in intestinal tight junctions. Interestingly, this pathway could also be implicated in celiac immune response [120], a condition common in ASD. WASF2 is also involved in "bacterial invasion of epithelial cells" and in "regulation of actin cytoskeleton" pathways, controlling cellular actin dynamics [121]. Then again, "Adipocytokine signalling pathway" is involved in cytokine production and inflammation as well as "TGF-β signalling". This latter is emerging as a key regulator of nervous system physiology also involved in nervous system diseases and injury [122]. Moreover, pathways directly connected with ASD are: "Long term depression" that, together with "Long-term potentiation", are the key regulators of long-lasting synaptic plasticity at the basis of learning and memory and whose impairment is involved in many brain disorders, including ASD [123,124]; "Insulin signalling pathway" promotes neuronal circuit development and maturation to an extent that, since its safety/preliminary efficacy for the treatment of Rett syndrome has emerged, clinical trials to ameliorate ASD symptoms are ongoing [125].
Similarly, Hallmark analysis identified "Apical surface", "Fatty acid metabolism", and "IL2 STAT5 signaling" hallmarks. The first one refers to a group of proteins expressed on the apical surface of epithelial cells, including enterocytes, that play a crucial role in cell polarity and thus act together with tight junction proteins in maintenance of epithelial tissue integrity [126]. Instead, the "Fatty acid metabolism" gene set includes genes known to be affected by gut microbiota composition [127]: saturated fatty acids increase bile-tolerant bacteria and reduce microbial diversity, while unsaturated fatty acids, such as omega-3, exert anti-inflammatory activity that contributes to gut health. Thus, an overall diet scheduling fat intake to alleviate gastro-intestinal disorders has been suggested [128]. Finally, IL2/STAT5 signalling pathway has a role in differentiation and homeostasis of both pro-and antiinflammatory T cells, determining the molecular details of immune regulation [129,130].
The limit of this study is the small size of samples; nevertheless, this is the first study that integrates metagenomic data with small ncRNA profile to investigate the host-gut microbiota cross-talk in ASD and that detects and analyses the piRNA profile in stools from individuals with ASD.

Conclusions
Since multiple interactions occur between host and gut microorganisms as well as among different microbial strains, faecal samples represent a source of information that should be studied as a whole. A complex picture emerges from this study due to the abundance and variability of gut microbiota and to the heterogenicity of samples. Based on our results and those of literature, we tried to find a relationship between the dysregulated gut microbial strains and up-regulated sncRNAs in ASD. We identified novel fungi genera and confirmed the bacteria dysbiosis, highlighting the possible role of dysregulated microbiome metabolites in ASD aetiology, too. Moreover, for the first time, we profiled miRNAs and piRNAs in ASD stool samples, which targeted pathways with roles in ASD. Disentangling the host-microbiome cross-talk is crucial to understand the role of dysbiosis in ASD onset and to design diagnostic tools and personalised therapeutic interventions. Host-gut microbiota cross-talk in health and disease, mediated by small ncRNAs, needs further insights, and omic and multi-omic data integration studies play a key role for these purposes [131,132].
The results obtained should be confirmed in large ASD and Ctrl cohorts, and the effect of the up-regulated miRNAs and piRNAs should be studied on microbial cultures as well.    Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.