Gut microbiota diversity and composition in children with autism spectrum disorder: associations with symptom severity

Samples collection

Approximately 60 samples were collected from children who had not undergone antibiotic treatment at the Affiliated Hospital of Putian University. For all the ASD patients involved in this study, their samples were gathered during outpatient consultations. The control group samples were sourced from neurotypical children who had not received any medications within one month prior to sample collection. Fresh fecal samples were obtained through rectal swabs. Subsequently, certain fecal samples from ASD patients and neurotypical children were excluded as a result of insufficient sample weight or suboptimal sample quality following DNA amplification.

Consequently, only 46 fecal samples were utilized for further sequencing and analysis, with 16 originating from neurotypical children, 15 from children with moderate/severe autism, and 15 from children with mild autism. The autistic children were admitted to the Affiliated Hospital of Putian University. Subsequently, the diagnosis of ASDs was established in accordance with the Diagnostic and Statistical Manual of Mental Disorders, 5th Edition (First, 2013). Additionally, evaluations were conducted using the Autism Diagnostic Observation Schedule and the Autism Behaviour Checklist. The scores of the Childhood Autism Rating Scale (CARS) (Schopler et al., 1980) were calculated by an experienced child neuropsychiatrist. The samples were stored in five ml tubes and promptly frozen at −80 °C until required for use. This study received approval from the Ethics Committee of the Affiliated Hospital of Putian University (Approval No.: 2023090-1), and the informed consent forms were signed by the parents. The detailed information and scoring criteria of all samples are presented in Table 1.

Main reagents and instruments

This study employed the following reagents and instruments: The FastPure Stool DNA Isolation Kit, provided by Shanghai Majorbio Bio-Pharm Technology Co., Ltd (Shanghai, China). The FastPfu Polymerase, sourced from Beijing TransGen Biotech Co., Ltd (Beijing, China). The NEXTFLEX Rapid DNA-Seq Kit, manufactured by Bioo Scientific (Austin, TX, USA). The T100 Thermal Cycler, produced by Bio-Rad Laboratories, Inc. (Hercules, CA, USA). The JY600C Double Stable Time Electrophoresis Apparatus, supplied by Beijing Junyi Dongfang Electrophoresis Equipment Co., Ltd (Beijing, China). The Synergy HTX, made by Biotek (Winooski, VT, USA). The Illumina’s NextSeq 2000 PE300 platform, which was utilized for high-throughput sequencing and the sequencing work was completed by Shanghai Majorbio Bio-Pharm Technology Co., Ltd (Shanghai, China).

Sample DNA extraction

The total genomic DNA of the microbial community present in the fecal samples of children with ASD (designated as the disease group) and those of neurotypical children (termed the neurotypical group) was extracted in strict accordance with the operating instructions provided by the FastPure Stool DNA Isolation Kit (Vazyme, Nanjing, China). Subsequently, the integrity of the extracted genomic DNA was examined through 1% agarose gel electrophoresis. Meanwhile, the DNA concentration and purity were quantified using the NanoDrop2000 instrument (manufactured by Thermo Fisher Scientific, Waltham, MA, USA).

PCR amplification and sequencing library construction

Taking the DNA extracted from the fecal samples as described above as a template, the V3–V4 variable region of the 16S rRNA gene was amplified through polymerase chain reaction (PCR). The upstream primer 338F (5′-CTCCTACGGGAGGCAGCAG-3′) and the downstream primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′), which carried the Barcode sequence, were utilized for this amplification. For amplification of the ITS gene, the upstream primer ITS3F (5′-GCATCGATGAAGAACGCAGC-3′) and the downstream primer ITS4R (5′-TCCTCCGCTTATTGATATGC-3′) were used. The PCR amplification program was set as follows: an initial pre-denaturation step at 95 °C for 3 min, followed by 27 cycles. Each cycle comprised denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s. Subsequently, a stable extension was carried out at 72 °C for 10 min, and finally, the samples were stored at 4 °C. The PCR instrument used in this process was the T100 Thermal Cycler.

The PCR reaction system using the TransStart FastPfu DNA polymerase purchased from TransGen (2* TransStart FastPfu PC, AS221-02, Beijing, China) and carried out according to the manufacture’s instruction, which includes following components: DNA polymerase (0.4 µL), template DNA (10 ng), the upstream and downstream primer were both 0.8 µL (5 µM), 2.5 mM dNTPs (2 µL), 5XTransStart FastPfu buffer (4 µL), and then adjust the total volume to 20 µL using the proper solvent. Each sample was replicated three times to ensure reproducibility. The PCR products of the same sample were combined and then recovered by means of 2% agarose gel electrophoresis. Subsequently, the recovered products were purified. The size of the resulting band fragments was determined by 2% agarose gel electrophoresis, and the recovered products were detected and quantified using the Synergy HTX device. The purified PCR products were then utilized to construct a sequencing library with the help of the NEXTFLEX Rapid DNA-Seq Kit. This construction process involved several key steps: first, linker ligation was carried out; then, magnetic beads were employed to screen and eliminate linker self-ligation fragments; next, the library template was enriched through PCR amplification; and finally, the PCR products were recovered using magnetic beads to obtain the final, ready-to-use library. Sequencing was performed on the Illumina’s NextSeq 2000 PE300 platform (undertaken by Shanghai Majorbio Bio-Pharm Technology Co., Ltd, Shanghai, China). After sequencing, the original data was uploaded to the NCBI SRA database for further analysis and sharing within the scientific community (PRJNA1207206) and the data was also deposited in the FigShare under the DOI: https://doi.org/10.6084/m9.figshare.28303508.

High-throughput sequencing data analysis

The paired-end raw sequencing reads underwent quality control procedures using fastp (available at https://github.com/OpenGene/fastp), and were subsequently assembled with FLASH (accessible at http://www.cbcb.umd.edu/software/flash, version 1.2.11). The detailed operations were carried out as follows: bases at the read termini with a quality value lower than 20 were filtered out. A window of 50 base pairs (bp) was set, and if the average quality value within this window fell below 20, bases were truncated from the end of the window. Reads with a length less than 50 bp after quality control were filtered, and those containing N bases were removed. Based on the overlap relationship between the paired-end reads, paired reads were merged into a single sequence, with a minimum overlap length of 10 bp. The maximum allowable mismatch ratio within the overlap region of the assembled sequence was set at 0.2, and sequences not meeting this criterion were screened out. Samples were distinguished according to the barcode and primers at both ends of the sequence, and the sequence direction was adjusted. The permitted number of barcode mismatches was 0, while the maximum number of primer mismatches was 2.

The optimized sequences underwent quality control and assembly before being processed with the DADA2 plugin for denoising (Callahan et al., 2016), which resulting in ASVs. To minimize bias from uneven sequencing depth in downstream alpha and beta diversity analyses, all samples were rarefied to 26,359 reads per sample. This normalization achieved an average Good’s coverage of 99.09%, indicating robust representation of microbial diversity. Based on the Silva 16S rRNA gene database (version 138) and unite9.0/its_fungi, the ASVs were taxonomically classified using the naive Bayes in Qiime2 (Bolyen et al., 2019). 16S rRNA function prediction analysis was conducted using PICRUSt2 (version 2.2.0), and the Funguild (version 1.0) was employed to evaluate the differential fungal microecology among different groups.

Statistical analysis

Data processing and statistical analyses were conducted using the Majorbio Cloud Platform (https://cloud.majorbio.com). Alpha diversity metrics, such as the Shannon and Ace indices, were calculated using mothur (http://www.mothur.org/wiki/Calculators, v1.30.2) (Schloss et al., 2009). Differences in alpha diversity across groups were evaluated using the Wilcoxon rank-sum test. To examine beta diversity, microbial community dissimilarities were visualized via principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS), both based on Bray-Curtis distances. Additionally, PERMANOVA was applied to determine the statistical significance of structural differences between sample groups.

Differences in gut bacterial and fungal diversity between ASD and neurotypical children

In total, 46 fecal specimens were successfully sequenced. Following the sequence quality control procedures, 1,713,362 sequences were acquired from all samples (Table S1). A combined total of 2,021 amplicon sequence variants (ASV) were obtained across all samples, with 715 ASV originating from the neurotypical group, 666 ASV from the severe ASD group, and 640 ASV from the mild to moderate ASD group. Notably, a higher number of ASV were obtained in the neurotypical group compared to the ASD patients. The rarefaction curve demonstrated that the sequencing depth was sufficient as it plateaued for all samples (Fig. S1), signifying that the majority of microbiota species had been captured.

Alpha diversity was computed to appraise the microbial diversity within the samples. The Kruskal–Wallis H test methodology was employed to calculate and comparatively analyze the inter-group differences in Alpha diversity of both bacteria and fungi among the distinct groups. In the 16S analysis, it was revealed that there were pronounced differences in bacterial diversity among the different groups. Specifically, the Ace index of children with severe ASD did not exhibit a significant difference from that of the neurotypical group (P > 0.05), yet it did display a significant difference from the Ace index of children with mild ASD (P < 0.05) (Fig. 1). Likewise, in the ITS analysis, no significant differences were observed among the three groups (P < 0.05).

However, upon analyzing the gut microbiota health index (Fig. 2), it was uncovered that significant differences existed in the microbiota health index among different groups, and this disparity bore no significant correlation with diversity. When comparing group G1 (the neurotypical group) to group G2 (the group of children with severe ASD), it was observed that the majority of samples in group G1 exhibited a higher gut microbiota health index (GMHI). Likewise, group G3 (the group of children with mild ASD) had a comparatively higher GMHI than group G2. Additionally, we conducted a further analysis of the differences in the microbiota dysbiosis index among the three groups of samples. It was found that, regardless of whether considering bacteria or fungi, there were pronounced differences between group G1 and the other two groups of children with ASD. Nevertheless, this difference appeared to be diminishing as the microbiota dysbiosis index (MDI) of group G1 and group G3 was closer. These results signify that the gut microbiota health index of children with ASD has undergone substantial alterations. In particular, the degree of microbiota dysbiosis is closely associated with the severity of ASD symptoms.

Beta diversity was utilized to assess the diversity of microbial community structure among different groups. As depicted in Fig. 3, the results of the 16S analysis indicate that there are significant differences in bacterial beta diversity between group G2 and both group G1 and group G3 (P < 0.05), whereas there is no significant difference between group G1 and group G3. Similarly, the ITS analysis reveals that there are significant differences in fungal beta diversity between group G2 and both group G1 and group G3 (P < 0.05), yet the difference between group G1 and group G3 becomes less pronounced. Such findings suggest that the gut microbial beta diversity of children with ASD might be related to the severity of their symptoms.

Furthermore, we employed NMDS and PCoA to dissect the differences in beta diversity of gut bacteria and fungi among the three groups of children (Fig. 4). Firstly, both the NMDS and PCoA analyses demonstrated that it was challenging to distinctly discriminate the gut bacteria (16S rRNA) of the three groups of children on the coordinate axes. However, relatively significant differences still prevailed (P < 0.05). Meanwhile, in contrast to group G2, more samples from group G3 and group G1 (the neurotypical group) clustered within the same region. Unlike bacteria, the fungal communities exhibited conspicuous differences among the three groups of children, and these distinctions could be visualized on the coordinate axes. The gut fungal communities from children of different groups were evidently clustered within their respective regions. It is noteworthy that the results of the PCoA analysis indicated that the gut fungi of children in group G2 were markedly different from those in the other two groups, implying that the gut microbiota of children with severe ASD symptoms might have experienced more profound changes in comparison to the neurotypical group.

Differences in gut bacterial and fungal communities between ASD and neurotypical children

The compositions of both bacteria and fungi were analyzed, ranging from the phylum level (Fig. S2) down to the genus level (Fig. 5), and then visualized in the form of a bar graph. In terms of bacteria, the five principal phyla present in all groups were Firmicutes, Bacteroidota, Proteobacteria, and Actinobacteria. However, the distribution patterns of the abundances of these phyla differed between children with ASD and neurotypical children. Moreover, significant differences were detected when comparing children with severe ASD to those with mild ASD showed in the pie graph (Fig. S3). Among the neurotypical children, Firmicutes (63.30%) constituted the predominant phylum, trailed by Proteobacteria (17.63%), Bacteroidota (9.13%), and Actinobacteria (9.06%). Similarly, in children with severe ASD, Firmicutes (79.86%) was also the leading phylum, yet its relative abundance was elevated in comparison to that of neurotypical children. Concurrently, the relative abundance of Proteobacteria (1.40%) was markedly reduced. Likewise, children with mild ASD also exhibited a relatively high abundance of Firmicutes (56.56%), followed by Actinobacteria (18.46%), Proteobacteria (13.06%), and Bacteroidota (10.87%). These findings suggest that, at the phylum level, group G3 was more akin to group G1 in terms of community composition.

Similarly, when examining the phylum level of fungi, Ascomycota, Basidiomycota, and unclassified_k_fungi were identified as the three main fungal groups in all samples. In group G1, Ascomycota accounted for 65.78%, trailed by unclassified_k_fungi (26.32%) and Basidiomycota (6.30%). In group G2, the proportion of Ascomycota increased to 87.59%, accompanied by a slight increment in Basidiomycota (8.27%), while the unclassified_k_fungi exhibited a significant decrease (3.84%). Notably, when compared to group G2, the relative abundance of Ascomycota in group G3 decreased to 73.18%, and that of unclassified_k_fungi increased to 15.92%. These results imply that, in terms of fungal composition at the phylum level, group G3 was more similar to group G1 than to group G2.

At the genus level, notable differences emerge in the gut microbiota of the three groups of children with varying severities of ASD (Fig. 5 and Fig. S4). Firstly, the outcomes of 16S rRNA gene sequencing analysis for bacteria reveal that among neurotypical children, the top five most prevalent genera within the microbiome are Escherichia-Shigella (accounting for 11.7%), Bifidobacterium (8.8%), Faecalibacterium (8.61%), Blautia (7.33%), and Bacteroides (7.26%). In children with severe ASD, the predominant microbial genera are Holdemanella (11.88%), Blautia (10.6%), Bifidobacterium (7.53%), and Faecalibacterium (6.29%). For children with mild ASD, the leading microbial genera comprise Bifidobacterium (17.91%), Escherichia-Shigella (10.86%), Bacteroides (8.88%), Blautia (7.58%), and Faecalibacterium (7.40%).

The results obtained from the ITS analysis of fungi indicate that among neurotypical children, the principal genera consist of Saccharomyces (accounting for 43.19%), Unclassified_k_Fungi (26.36%), Aspergillus (6.23%), Unclassified_o_Saccharomyces (4.97%), and a minor proportion of Candida (2.01%). In children with severe ASD, Candida exhibits a substantial increase, reaching up to 67.52%, whereas the relative abundance of Saccharomyces declines significantly (from 43.19% in neurotypical controls to 4.45% in severe ASD patients). Additionally, the fungal composition in these children also encompasses Aspergillus (7.22%), Cladosporium (4.22%), and Unclassified_k_Fungi (3.84%). Similarly, when compared with children having severe ASD, the main fungal constituents in children with mild ASD display a certain degree of restoration. This is manifested by the fact that the quantity of Candida is markedly lower than that in children with severe ASD, and the amount of Saccharomyces is greater than that in children with severe ASD. Meanwhile, regardless of whether in neurotypical children or those with mild ASD, the abundance of Unclassified_k_Fungi in the gut surpasses that in children with severe ASD. Consequently, from the vantage point of microbial composition, it appears that the principal gut microbiota of children with mild ASD bears a closer resemblance to that of neurotypical children, while the main gut microbiota of children with severe ASD diverges considerably from that of neurotypical children. This further suggests a potential association between the gut microbiota of children and relevant ASD symptoms.

To gain a more profound understanding of the situation, we decided to carry out a hierarchical clustering analysis on the top 50 principal microbial genera of both bacteria and fungi, which led to the generation of a species heatmap (Fig. 6). The heatmap indicates that, regardless of whether we are looking at bacteria or fungi, the gut microbiota of the three groups of children can generally be grouped together, with the clustering of fungi being more conspicuous. Interestingly, in terms of bacteria, group G1 appears to be relatively scattered, while groups G2 and G3 are relatively concentrated. It is rather difficult to determine whether group G1 is nearer to G2 or G3. However, upon conducting a detailed analysis of the ITS results, we would notice that the samples of G2 and G3 essentially form their own separate clusters. Moreover, some of the G1 samples are relatively close to the G3 samples and thus fall into the G3 cluster, such as H18, H4, H25, and so on. Notably, none of the G1 samples are included in the G2 cluster. These findings, taken together, suggest that the main composition of the microbiota in group G2 differs considerably from that of the gut microbiota of neurotypical children, and this difference might potentially lessen as the symptoms improve.

Analysis of functional in gut bacteria and fungi between ASD and neurotypical children

Subsequently, the normal biological functions of bacteria and fungi in both neurotypical individuals and those with ASD were analyzed using PICRUSt2 and FUNGuild, with a one-way ANOVA test being employed. To further explore the differences among each group, the Tukey–Kramer test was utilized. The 16S rRNA functional analysis of bacteria unveiled that a total of 263 KEGG Level 3 modules were present across all samples. At Kyoto Encyclopedia of Genes and Genomes (KEGG) level 1, no remarkable differences were spotted among the three groups. However, when closely examining KEGG level 2 (Fig. 7A), it was found that significant differences existed between Group G2 and both Group G3 and Group G1 (with P < 0.05). In contrast, the difference between Group G1 and Group G3 was not statistically significant (as P > 0.05). Specifically, the amino acid metabolism in G2 (P < 0.001) was substantially higher than that in G1 and G3, while no difference was detected between G1 and G3 in this regard. Conversely, the Metabolism of cofactors and vitamins in both G2 and G3 (P < 0.001) was markedly lower than that in G1.

The functional prediction analysis carried out by FUNGuild demonstrated that the taxonomic variations of fungi in the intestines of children changed in line with different ASD symptoms (Fig. 7B). It also indicated that significant differences were present between G2 and G1 (P < 0.05 for Global and overview maps), whereas no significant difference was noted between G1 and G3. Among them, children with severe ASD had a higher proportion of undefined saprotroph (P < 0.05). Although other fungal taxonomic differences were not statistically prominent, a tendency could still be observed, suggesting that G1 was closer to G3. Collectively, these predictive analyses of bacterial functions and fungal functions provided solid evidence that the intestinal microbiota of children with severe ASD symptoms deviated from that of neurotypical children, while the intestinal microbiota of children with mild ASD seemed to bear a closer resemblance to that of neurotypical children. Hence, it can be concluded that the composition and functions of intestinal bacteria and fungi exerted a certain degree of influence on ASD symptoms.