Eight male and eight female Wistar rats aged 10–12 weeks and weighing 270–290 g were purchased from Vital River Laboratory Animal Technology (Beijing, China). The animals were housed at 23 $\pm$ 2 ${}^{\circ }$C and 50% $\pm$ 10% humidity under a 12-h light/dark cycle. The experimental procedures were in line with the USA National Institutes of Health Guide for the Care and Use of Laboratory Animals and approved by the Animal Care and Use Committee of Peking University People’s Hospital (2019PHE007).

We constructed the animal model in accordance with previous studies [44, 45]. Female rats were paired with male rats overnight. The day after a vaginal plug was detected was considered embryonic day 0.5 (E 0.5). Pregnant rats were housed individually and randomly distributed into two groups: the VPA-treated group (n = 4) and the normal saline (NS) group (n = 4). Pregnant rats in the VPA-treated group were intraperitoneally injected with VPA (Sigma: P4543) at a dose of 600 mg/kg on E 12.5. In the NS group, the rats received the same volume of NS. After weaning at postnatal day 21 (PND 21), every 4-5 pups of the same group and sex were raised in the same cage until the end of fecal sample collection.

We randomly selected 1 male and 1 female pup from each litter for our experiment. A total of 16 pups from litters in the VPA-treated group and NS group were included and divided into 4 groups, with 4 animals each: the (1) male control group (mCtrl group), (2) male VPA group (mVPA group), (3) female control group (fCtrl group), and (4) female VPA group (fVPA group).

Pups were individually placed into clean cages. After defecation, we used sterile tweezers to collect the fecal pellets into sterile cryotubes and immediately transferred them into liquid nitrogen. After collection of feces, all the fecal samples were stored at –80 ${}^{\circ }$C until analysis by untargeted metabolomics and 16S rRNA sequencing.

16S rRNA gene sequencing was performed based on the method described previously with some modifications [46]. Total fecal genomic DNA was extracted with an E.Z.N.A Stool DNA Kit (OMEGA, USA) and quantified by a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The V3-V4 region of the 16S rRNA gene was amplified with the primers 341F (5${}^{\prime }$-CCTAYGGGRBGCASCAG-3${}^{\prime }$) and 806R (5${}^{\prime }$-GGACTACNNGGGTATCTAAT-3${}^{\prime }$). Twenty microliter PCR mixtures were set up with 4 $\mu$L 5$\times$FastPfu buffer, 2 $\mu$L deoxynucleoside triphosphates (dNTPs) (2.5 mM), 0.8 $\mu$L each primer, 0.4 $\mu$L FastPfu polymerase, and template DNA (10 ng). PCR was performed in a GeneAmp 9700 (ABI) thermocycler (Thermo Fisher Scientific, Waltham, MA, USA) with Phusion High-Fidelity PCR Master Mix. The parameters were set as follows: initial denaturation at 95 ${}^{\circ }$C for 5 min; 27 cycles of denaturation at 95 ${}^{\circ }$C for 30 s, annealing at 55 ${}^{\circ }$C for 30 s, and elongation at 72 ${}^{\circ }$C for 45 s; and an additional elongation at 72 ${}^{\circ }$C for 10 min.

PCR products were detected by 2% agarose gel electrophoresis and purified by a quick gel extraction kit (Qiagen, Dusseldorf, Germany). Library pools were produced with equal amounts of PCR products with a TruSeq Nano DNA LT sample prep kit (Illumina, San Diego, California, USA); these pools were amplified through paired-end sequencing on the Illumina MiSeq PE300 platform.

The quality control (QC) of original data was carried out using Trimmomatic v0.39 software (http://www.usadellab.org/cms/index.php?page=trimmomatic). Based on the overlap (minimum: 10 bp) between paired-end (PE) reads after QC, PE reads were assembled using Flash v1.2.11 software (McKusick-Nathans Institute of Genetic Medicine, Johns Hopkins University School of Medicine, USA). QIIME v1.9.1 software (Knight Lab in University of California, San Diego, USA) (QIIME allows analysis of high-throughput community sequencing data) was adopted for processing, and VSEARCH v2.14.1 software (Department of Informatics, University of Oslo, Norway) (VSEARCH: a versatile open-source tool for metagenomics) was used to detect chimera sequences. We used the Uclust method in the QIIME v1.9.1 software package to perform operational taxonomic units (OTUs) clustering analysis according to 97% sequence similarity level. Rarefaction was not applied in our research. On the basis of the Silva database Release 138 (http://www.arb-silva.de), taxonomic annotations were made for the OTUs of each sample using a confidence threshold of 70%. LefSe (http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=lefse_upload) was used to evaluate the effect of each taxon on discrimination between different groups. The Shannon, Simpson, Chao1 and ACE indices of microbial communities were calculated by Mothur (https://mothur.org/). Principal coordinates analysis (PCoA), based on Bray-Curtis dissimilarity, was analyzed by permutational multivariate analysis of variance (PERMANOVA) in R language v3.2.1 (University of Auckland, New Zealand) to determine the different bacterial communities of beta-diversity.

Fecal samples (100 mg each sample) were triturated in liquid nitrogen and vortexed with 300 $\mu$L lysis buffer (80% methanol solution). The mixtures were homogenized at 45 Hz for 4 min and ultrasonicated at 4 ${}^{\circ }$C and 40 kHz for 5 min. Then, we placed the mixtures at –20 ${}^{\circ }$C for 30 min and centrifuged at 4 ${}^{\circ }$C and 14000 $\times$ g for 30 min. The supernatant of each sample was collected for UHPLC-MS analysis. The column was a ACQUITY UPLC HSS T3 column (100 mm $\times$ 2.1 mm i.d., 1.8 $\mu$m; Waters). The mobile phase A was water (containing 0.1% formic acid), and the mobile phase B was acetonitrile (contains 0.1% formic acid). The flow rate, injection volume, and column temperature were 0.25 mL/min, 5 $\mu$L, and 35 ${}^{\circ }$C, respectively. A QC sample (prepared after metabolite extraction, mixing equal volumes of all samples) was used to evaluate the stability of the analytical system and ensure the reliability of the results.

UHPLC-MS was performed on a UPLC I-class instrument (Waters, ACQUITY UPLC). It was equipped with a binary solvent delivery manager, an autosampler and a mass spectrometer (Waters, Vion IMS QTof) with an electrospray interface. Sample mass spectrometric data were collected in either positive or negative ionization mode. We set the source temperature, desolvation temperature and desolvation gas flow to 120 ${}^{\circ }$C, 500 ${}^{\circ }$C and 800 L/h, respectively. Centroid data were assembled from 100–1200 m/z with a 0.3 s scan time. The capillary voltages were 3.2 kV and –2.2 kV in positive ionization mode and negative ionization mode, respectively. Before pattern recognition, the original data were subjected to preprocessing. Baseline filtering, peak identification, integration, retention time correction, peak alignment and normalization were performed using Progenesis QI 2.3 software. Then, the METLIN database (https://metlin.scripps.edu), Human Metabolome database (http://www.hmdb.ca/) and ChemSpider database (www.chemspider.com) were applied for metabolite annotation. The fragmentation mass spectrum of metabolites was matched to fragmentation mass spectra in library and scored. The results with higher score were accepted as the identification results.

Principal component analysis (PCA) and orthogonal partial least squares discrimination analysis (OPLS-DA), performed by the SIMCA-P software package (v13.0, Umetric, Umea◦, Sweden), were used to observe the degree of dispersion and distinguish the overall differences in metabolic profiles, respectively. In OPLS-DA model, metabolites with variable importance in projection (VIP) $>$1 were set as differential variables [47, 48]. The method of 100 replacement tests was carried out to examine the fitting effect of OPLS-DA and prevent the model from overfitting. Variables with VIP $>$1.0 and p 0.05 calculated by Student’s t-test were defined as differential metabolites. All differentially abundant metabolites were annotated in pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) [49]. Pathway enrichment analysis was performed using the Python software package (https://docs.scipy.org/doc/scipy/). Receiver operator characteristic (ROC) analysis was carried out using R language (version 3.2.1) to calculate the area under the curve (AUC) and evaluate the degree of influence of metabolites on the distinction between the VPA and Ctrl groups in different sexes.

Quantitative sequencing data with a nonnormal distribution were analyzed by the Wilcoxon rank sum test, and p-values were checked using the Benjamini Hochberg false discovery rate (FDR). Data were presented as the mean $\pm$ standard deviation (SD). The Pearson correlation coefficient was used to express the correlation degree between characteristic genera and characteristic fecal metabolites. GraphPad Prism 8.0.1 software (GraphPad Software, San Diego, California, USA) was used for data analysis and graphing. A double-tailed p 0.05 was accepted as statistically significant.

To detect the taxonomical characterization of the GM in male and female VPA rats, fecal samples were analyzed by 16S rRNA gene sequencing. We obtained 943,972 sequences and 390,572,407 bases from 16 samples with an average of 58,998 sequences per sample. A total of 1501 OTUs were annotated in all samples (1501 at the phylum level; 1499 at the class level; 1496 at the order level; 1448 at the family level; 861 at the genera level; 115 at the species level). The analysis of $\alpha$-diversity showed no significant differences in either sex between Ctrl and VPA rats using the Chao1 estimator (Fig. 1A), ACE estimator (Fig. 1B), Shannon index (Fig. 1C) and Simpson index (Fig. 1D). PCoA showed that the overall composition of the GM was not different between the mCtrl and mVPA groups (R = 0.193, p = 0.151, Fig. 1E), while between the fCtrl and fVPA groups, the overall composition of the GM was significantly different (R = 0.267, p = 0.036, Fig. 1E). At the phylum level, the mVPA group showed an increased Firmicutes/Bacteroides (F/B) ratio compared with that of the mCtrl group (Fig. 1F), while the fVPA group showed a decreased F/B ratio compared with that of the fCtrl group (Fig. 1F). The LefSe results showed that at the genus level, Ruminococcus_2, Candidatus_Soleaferrea, Desulfovibrio, Candidatus_Saccharimonas, Intestinimonas, [Eubacterium]_xylanophilum_group, [Eubacterium]_brachy_group and [Bacteroides]_pectinophilus_group were characteristic genera in the mVPA group (Fig. 1G), while Ruminiclostridium, Acetatifactor, Desulfovibrio, [Eubacterium]_xylanophilum_group and Candidatus_Saccharimonas were characteristic genera in the fVPA group (Fig. 1H). The relative abundances of characteristic genera of either sex of VPA rats showed no significant differences (Tables 1,2). These results revealed the taxonomical characterization of GM constituents in male and female VPA rats, with obvious sex-specific differences.

Fig. 1.

The difference in gut microbiota in different sexes between the Ctrl and VPA rats. (A) Chao1 index. (B) ACE index. (C) Shannon index. (D) Simpson index. (E) PCoA chart. The X and Y axes represent two selected principal axes. The percentage represents the explanatory value of the principal axis for the difference between samples. (F) Column chart of species composition at the phylum level. (G–H) LDA chart. The impact of the representative species abundance on the differences between (G) the mVPA and mCtrl groups or between (H) the fVPA and fCtrl groups. mCtrl, male control group; mVPA, male VPA group; fCtrl, female control group; fVPA, female VPA group.

Untargeted metabolomics was performed to investigate the characteristic fecal metabolites in male and female VPA rats. PCA (Fig. 2A–D), OPLS-DA (Fig. 2E–H) and model validation (Fig. 2I–L) revealed that in both sexes, the distribution of samples between VPA and Ctrl rats was different, implying a significant difference in fecal metabolites. Between the mCtrl and mVPA groups, we detected 51 annotated differentially abundant metabolites (Fig. 2M,N, Supplementary Table 1). A total of 124 annotated differentially abundant metabolites were detected between the fCtrl and fVPA groups (Fig. 2O,P, Supplementary Table 2). Detailed information on these annotated differentially abundant metabolites is shown in Supplementary Table 1 and Supplementary Table 2. These results indicated that fecal metabolite levels were altered in both sexes of adolescent VPA rats.

Fig. 2.

The difference in fecal metabolites in different sexes between the Ctrl and VPA rats. (A–D) Principal component analysis (PCA). (E–H) Orthogonal partial least squares discrimination analysis (OPLS-DA). (I–L) OPLS-DA model validation. (M–P) Volcano map of differentially abundant metabolites. mCtrl, male control group; mVPA, male VPA group; fCtrl, female control group; fVPA, female VPA group.

We further analyzed these annotated differentially abundant metabolites in KEGG pathways. Between the mVPA and mCtrl groups, among the 51 annotated differentially abundant metabolites, 7 metabolites participated in 10 specific KEGG pathways (Table 3, Fig. 3A) and were significantly enriched in 1 KEGG pathway named the pathway for cyanoamino acid metabolism (p = 0.0027, FDR = 0.0159, Fig. 3A). Among the 7 metabolites, we focused on 4 metabolites with VIP $>$2 (Table 3). ROC analysis showed that the AUC values of these 4 metabolites were larger than 0.85 (Fig. 3B–E), suggesting their potential as important biological indicators to distinguish between the mVPA and mCtrl groups. Between the fVPA and fCtrl groups, among the 124 annotated differentially abundant metabolites, 15 metabolites participated in 34 specific KEGG pathways (Table 4, Fig. 4A). However, these metabolites were not significantly enriched in any KEGG pathways (Fig. 4A). Among the 15 metabolites, we focused on 8 metabolites with VIP $>$2 (Table 4). ROC analysis showed that the AUC values of these 8 metabolites were larger than 0.85 (Fig. 4B–I), suggesting their potential as important biological indicators to distinguish between the fVPA and fCtrl groups. These results revealed the potential mechanism of ASD induced by prenatal exposure to VPA as well as the characteristic fecal metabolites of male and female adolescent VPA rats. Sex-specific differences were obvious among these results.

Fig. 3.

KEGG analysis of fecal metabolites in male VPA rats. (A) KEGG pathway enrichment column chart in male VPA rats. The abscissa is the name of KEGG pathway. (B–E) Receiver operator characteristic curves of the four characteristic fecal metabolites. (B) Hesperetin. (C) Linamarin. (D) Melphalan. (E) Protoanemonin.

Fig. 4.

KEGG analysis of fecal metabolites in female VPA rats. (A) KEGG pathway enrichment column chart in female VPA rats. The abscissa is the name of KEGG pathway. (B–I) Receiver operator characteristic curves of the eight characteristic fecal metabolites. (B) 4-(Methylnitrosamino)-1-(3-pyridyl)-1-butanol. (C) Allocholic acid. (D) Androstenedione. (E) Brassicasterol. (F) cis-Acetylacrylate. (G) D-1-Piperideine-2-carboxylic acid. (H) LysoPC(14:1(9Z)). (I) Palmitoyl glucuronide.

We used Pearson correlation coefficients to assess the correlation between characteristic genera and fecal metabolites. Fig. 5A shows the correlation between 4 characteristic metabolites and 8 characteristic genera in male VPA rats. The abundance of hesperetin and protoanemonin negatively correlated with that of Candidatus_Soleaferrea and Desulfovibrio. In addition, the abundance of protoanemonin negatively correlated with that of Ruminococcus_2. Fig. 5B shows the correlation between 8 characteristic metabolites and 5 characteristic genera in female VPA rats. The abundance of Ruminiclostridium, Acetatifactor and [Eubacterium]_xylanophilum_group positively correlated with that of brassicasterol, allocholic acid, palmitoyl glucuronide and D-1-piperideine-2-carboxylic acid but negatively correlated with that of 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanol. The abundance of Acetatifactor and [Eubacterium]_xylanophilum_group positively correlated with that of LysoPC(14:1(9Z)) but negatively correlated with that of androstenedione. In addition, the abundance of cis-acetylacrylate negatively correlated with that of Ruminiclostridium, Acetatifactor and Candidatus_Saccharimonas. Correlation analysis revealed a tight relationship between the GM and fecal metabolites in VPA rats of both sexes. The results of the correlation between characteristic genera and characteristic fecal metabolites in VPA rats also showed sex-specific differences.

Fig. 5.

Correlation analysis of gut microbiota and fecal metabolites in male and female VPA rats. (A) Correlation heatmap between sixteen characteristic fecal metabolites and eight characteristic genera in abundance in male VPA rats. (B) Correlation heatmap between sixteen characteristic fecal metabolites and eight characteristic genera in abundance in female VPA rats. The legend on the right is the color range of different R values. The abscissa represents characteristic genera, and the ordinate represents characteristic fecal metabolites. ${}^{*}$p 0.05, ${}^{**}$p 0.01, ${}^{***}$p 0.001.