Study population
Fecal samples were collected from women enrolled in the ENDIA study across five Australian States between February 2013 and October 2017 [22]. ENDIA is a prospective, pregnancy-birth cohort study that follows 1500 Australian children who have a first-degree relative with T1D. Fifteen women (16 pregnancies) without T1D (non-T1D) and 16 with T1D had each provided three fecal samples across pregnancy (total of 96 samples) (Table 1). These were analysed by shotgun WMS (Figure S1A). A second sample set comprising a total of 354 fecal samples collected longitudinally across 145 pregnancies from 51 non-T1D women and 94 T1D women, including the 96 samples analysed by WMS, was analysed by 16S rRNA gene sequencing (Figure S1B). Table 1 summarizes and compares characteristics of the non-T1D and T1D pregnancies.
Whole metagenomic sequencing and 16S sequencing
For the WMS dataset, 43,809,936 ± 9,156,972 (mean ± SD) paired end reads per sample were obtained using an Illumina NovaSeq 6000. Raw reads (SRA accession: PRJNA604850) were pre-processed using KneadData bioBakery tool [23] to eliminate human DNA sequences and filter sequences with poor quality which on average removed 6% of the reads. Reads of 41,013,121 ± 8,739,281 (mean ± SD) were obtained after quality control and read filter steps (Excel file E0). For the 16S rRNA dataset, 42,104 ± 17,296 (mean ± SD) paired end reads per sample were obtained with an Illumina MiSeq. QIIME2 was used to demultiplex and quality-filter the raw sequences (SRA accession: PRJNA604850), generating 24,777 ± 10,749 (mean ± SD) reads per sample (Excel file E0).
Taxonomic diversity and composition of the gut microbiome in non-T1D and T1D women during pregnancy
Metagenomic sequences of the women who provided three fecal samples across trimesters were analysed using MetaPhLan2 implemented within the HUMAnN2 pipeline. Overall, 298 bacterial species were identified with an average of 97 ± 11 (mean ± SD) species per sample.
The top 25 most abundant species accounted for more than 50% of the gut microbiome composition of each subject in any given trimester (Figure 1). The Firmicutes/Bacteroidetes phylum ratio across trimesters was unchanged in non-T1D women but decreased in women with T1D (Figure S2).
Alpha diversity (observed richness or number of species) per sample was calculated (Excel file E1) and generalized estimating equations (GEE) were applied to test for differences between women without and with T1D, between trimesters and to determine if there was an interaction between T1D status and trimester. No significant interactions or differences in richness were found (Figure S3, Excel file E1).
For beta diversity analysis, Bray-Curtis coefficients were calculated between sample pairs, ordinated and plotted by principal coordinate analysis (PCoA; Figure 2). To test for differences in beta diversity, a repeated-measure aware permutational analysis of variance (RMA-PERMANOVA) of the Bray-Curtis coefficients was performed on normalized data. This revealed a significant interaction between T1D status and time at all taxonomic levels. Therefore, differences between T1D and non-T1D women were assessed within trimester. No significant differences were detected in trimester 1. However, differences were significant at the strain level (P=0.02) in trimester 2 and at the strain (P=0.002), species (P=0.001), genus (P=0.015), family (P=0.029) order (P=0.034) and phylum (P=0.032) taxonomic levels in trimester 3 (Excel file E2).
Differences in beta diversity reflect differences in taxonomic composition. To identify differences in taxa between non-T1D and T1D women in pregnancy, differential abundance was analysed in limma; the Benjamini-Hochberg method was used to control the false discovery rate. Taxa were considered differentially abundant if the adjusted P-value expressed as the False Discovery Rate (FDR) was < 0.05; however, borderline significant differences, i.e. FDR from 0.05 to 0.1, were also reported. Note that only taxa for which the prevalence (i.e., proportion of samples with those taxa) was above 50% in at least one group and with a log2 fold-change (logFC) greater than 0.5 or less than -0.5 were considered. T1D status and trimester were combined into a single factor with six levels and comparisons of interests were defined as contrasts. Across all trimesters, the species Bacteroides eggerthii (FDR=0.033), Bacteroides clarus (FDR=0.005), Alistipes sp. AP11 (FDR=0.009), Escherichia coli (FDR=0.002), an unclassified Escherichia bacterium (FDR=0.042), Roseburia hominis (FDR=0.08) and Bacteroides uniformis (FDR=0.088), as well as the family Enterobacteriaceae (FDR=0.09) and the order Enterobacteriales (FDR=0.03), were increased in women with T1D (Figure 3; Excel file E3). On the other hand, the species Bacteroides massiliensis (FDR=0.043) and the family Prevotellaceae (FDR=0.05) were decreased in T1D (Figure 3; Excel file E3).
Differences between non-T1D and T1D women were also assessed within trimesters. In trimester 1, species Escherichia coli (FDR=0.009), Alistipes sp. AP11 (FDR=0.063), Bacteroides clarus (FDR=0.073) and Bacteroides salyersiae (FDR=0.073), and the order Enterobacteriales (FDR=0.0005) were increased in T1D, while Bacteroides massiliensis (FDR=0.016) was decreased in T1D (Figure 3; Excel file E3). In trimester 2, only the family Prevotellaceae (FDR=0.065) and the phylum Firmicutes were decreased in T1D (Figure 3; Excel file E3).
Most differences between non-T1D and T1D women were in trimester 3, in which species Bacteroides eggerthi (FDR=0.012), B. clarus (FDR=0.02), E. coli (FDR=0.012), an unclassified Escherichia bacterium (FDR=0.083), R. hominis (FDR=0.007), B. uniformis (FDR=0.012), B. salyersiae (FDR=0.011), B. caccae (FDR=0.012), Parabacteroides distasonis (FDR=0.012), B. faecis (FDR=0.088), Sutterella wadsworthensis (FDR=0.083) and the family Bacteroidaceae (FDR=0.078) were increased in T1D, and species Bifidobacterium adolescentis (FDR=0.025) and Ruminococcus bromii (FDR=0.09) and the order Bifidobacteriales (FDR=0.002) were decreased (Figure 3; Excel file E3). Five of the 15 differentially abundant species detected were among the 25 most abundant species in the entire dataset (Figure 1). Differential abundance data are summarised in Excel file E3.
In order to expand on these observations, we sequenced DNA of the V4 region of the 16S rRNA marker gene in a larger dataset comprising 354 samples from 51 non-T1D and 94 T1D women, including the original 96 samples. Sequences were processed with QIIME2 to obtain a feature per sample table and further analysis was performed in R (see Methods). Features were agglomerated based on a phylogenetic tree into an overall total of 349 operational taxonomic units (OTUs) with a mean ± SD of 82 ± 24 (mean ± SD) OTUs per sample (Excel file E0). Fourteen of the 25 most abundant species detected by WMS were also within the 25 most abundant OTUs in the larger 16S rRNA dataset. Moreover, one OTU had 100% sequence identity with both Bacteroides vulgatus and B. dorei and therefore these two species could not be distinguished by 16S rRNA gene sequencing.
Consistent with the metagenomic analysis, alpha diversity did not differ by T1D status in the larger 16S rRNA dataset in any trimester (Excel file E1). Also, in accord with the WMS analysis, an interaction between T1D status and time was significant for beta diversity at the genus and family taxonomic levels and borderline significant at the OTU taxonomic level (Excel file E2). Therefore, differences in beta diversity were also assessed within trimester. The gut microbiome composition differed significantly between non-T1D and T1D women within trimester 2 and 3 at the OTU, genus and family taxonomic levels, but not within trimester 1 (Excel file E2). At the order and phylum taxonomic levels, T1D status and the interaction between T1D status and time were not significant.
Only one OTU classified to the species Bacteroides caccae was differentially abundant in trimester 3 (FDR=0.02) (Figure 4A; Excel file E3). Generally, the prevalence of taxa was lower with 16S rRNA sequencing than WMS. To confirm the differentially abundant taxa detected between non-T1D and T1D women by WMS, relative abundances by 16S rRNA sequencing were agglomerated at family, order and phylum taxonomic levels. Most of the taxa identified in the 16S rRNA dataset exhibited the same trend in abundance between non-T1D and T1D as observed with WMS (Figures 4B-F).
Effect of time and other factors on the gut microbiome during pregnancy
WMS revealed no significant differences in alpha diversity in non-T1D or T1D women due either to time analysed as days into gestation or by trimester, i.e., as continuous or categorical variables, respectively (Excel file E1). Due to a significant interaction between T1D status and time, differences in beta diversity between trimesters were assessed in non-T1D and T1D women separately (Excel Table E2). No significant differences in beta diversity due to trimester were observed at any taxonomic level in non-T1D or T1D women (Excel Table E2). However, significant differences were present by time analysed as days into gestation at the strain and genus levels (and borderline differences at the species and family levels) in T1D samples and at the strain, species, genus and family levels in non-T1D samples (Excel file E2).
In T1D women, differences in the abundance of specific taxa between trimesters were identified. Thus, the abundance of strain Collinsella aerofaciens St. GCF000169035 (FDR=0.087) decreased from trimester 1 to 2 (Figure 5; Excel file E4). Lachnospiraceae bacterium 8 1 57FAA St. GCF000185545 (FDR=0.04), Akkermansia muciniphila St. GCF000020225 (FDR=0.04), Ruminococcus bromii St. GCF000209875 (FDR=0.002) and unclassified strains from the species Roseburia intestinalis (FDR=0.04), Oxalobacter formigenes (FDR=0.04), Bifidobacterium adolescentis (FDR=0.04) and Parabacteroides goldsteinii (FDR=0.04) decreased from trimester 1 to 3 (Figure 5; Excel file E4).
Also, in T1D women the abundance of orders Clostridiales (FDR=0.011), Coriobacteriales (FDR=0.025), Pasteurellales (FDR=0.002), Verrucomicrobiales (FDR=0.004), Erysipelotrichales (FDR=0.002), Selenomonadales (FDR=0.002) and Enterobacteriales (FDR=0.002) decreased from trimester 1 to 2 (Figure 5; Excel file E4). The orders Pasteurellales (FDR=0.025), Verrucomicrobiales (FDR=0.004), Enterobacteriales (FDR=0.04) and Bifidobacteriales (FDR=0.06) decreased from trimester 1 to 3 (Figure 5; Excel file E4). Only the abundance of the order Lactobacillales (FDR=0.018) and an unclassified strain from species B. massiliensis (FDR=0.087) increased from trimester 1 to 2 (Figure 5; Excel file E4). In non-T1D women, the abundance of an unclassified strain from the species Haemophilus parainfluenzae (FDR=0.045) and the family Pasteurellaceae (FDR=0.04) decreased from trimester 1 to 3 and from trimester 2 to 3, respectively, while the abundances of the order Bifidobacteriales (FDR=0.05) and the phylum Actinobacteria (FDR=0.03) increased from trimester 1 to 3 (Figure 5; Excel file E4).
Similar to WMS, 16S rRNA sequencing revealed no differences in alpha diversity by time (Excel file E1) or in beta diversity by trimester (Excel file E2). Differences in beta diversity by time were observed in T1D women only at the family level (Excel file E2). In the differential abundance analysis, within T1D women only, the abundance of an OTU classified to the genus Lachnospiraceae (FDR=0.074), the family Clostridiaceae 1 (FDR=0.047) and the phylum Bacteroidetes (FDR=0.005), increased from trimester 1 to 2 (Excel file E4). Family Peptococcaceae (FDR=0.037) and the order Lactobacillales (FDR=0.028) increased from trimester 2 to 3 (Excel file E4). An OTU classified to the Oscillibacter genus (FDR=0.009), the family Bacteroidaceae (FDR=0.006) and the phylum Bacteroidetes (FDR=0.012) increased from trimester 1 to 3 (Excel file E4). Within non-T1D women only, the orders Lactobacillales (FDR=0.011) and Pasteurellales (FDR=0.022) increased from trimester 1 to 3 (Excel file E4).
No significant associations were found between beta diversity and other factors included in the models, viz. age at conception, body mass index (BMI), parity and human leukocyte antigen (HLA) class II genotype, for either metagenomic or 16S rRNA data with the exception of parity which was significantly different only in trimester 1 at the order taxonomic level, for the 16S rRNA data (Excel file E2). As expected, non-T1D and T1D women differed in serum 1,5-AG, a marker of short-term glycemic control [24] (Table 1), but in T1D women serum 1,5-anhydroglucitrol (1,5-AG) was not related to beta diversity (Excel file E2). The frequency of pre-eclampsia was higher in the larger group of 145 women (Table 1) but was not related to beta diversity (Excel file E2). The frequency of elective Caesarean section was higher in T1D women (Table 1) but was not related to beta diversity (Excel file E2). Non-T1D and T1D women also differed in total carbohydrate but not fiber intake (Table 1). In the larger group of 145 pregnancies, total carbohydrate intake was lower in those with T1D (P=0.016), most likely reflecting dietary advice. Based on 16S data, interactions were found between carbohydrate and fiber intake and time; therefore, beta diversity was assessed per trimester. Carbohydrate intake was related (borderline P-values greater than 0.05 but less than 0.1) to beta diversity at family, order and phylum levels only in trimester 2. Fiber intake was significantly related to beta diversity at genus and family levels only (borderline significance for order and phylum taxonomic levels) in trimester 3 (Excel file E2). No relationships between either carbohydrate or fiber intake and beta diversity were found with WMS data.
Functional annotation of gut microbiome taxa
For functional profiling, metagenomic sequences were processed with HUMAnN2. Sequences were annotated, gene abundances calculated and regrouped into KO (Kegg Orthology) and MetaCyc reaction functional categories, and complete metabolic pathways were quantified obtaining a total of three functional profiles. A total of 5,480 KO, 3,014 metaCyc reactions and 420 complete pathways, were obtained. No significant differences in richness were detected between T1D and non-T1D women or across pregnancy in any of the three functional profiles (Excel file E1). On the other hand, in the analysis of beta diversity, the interaction between T1D status and time was significant. Therefore, differences between non-T1D and T1D women were assessed in each trimester and were only significant in trimester 3 for pathways, metaCyc reactions and KO functions (Figures 6; Excel file E2).
Differential abundance analysis revealed significant differences between non-T1D and T1D in specific pathways and enzymes only within trimesters 1 and 3 (Figure 7; Excel files E5-7). A comprehensive list of differentially abundant pathways, KOs and MetaCyc reactions is presented in Supplementary Excel files E5-10). Of interest, a pathway (CMP-3-deoxy-D-manno-octulosonate biosynthesis I) and two enzymes (K02852 and K01791) involved in the synthesis of bacterial lipopolysaccharides (LPSs) were enriched in T1D women (Figure 7; Excel files E5 and E7), with contributions from Alistipes shahii, B. thetaiotaomicron, B. ovatus, B. vulgatus, B. cellulosilyticus, Parabacteroides distasonis, Bacteroides uniformis among others (Figure S4). In addition, a KO and a MetaCyc reaction that were increased in T1D women in trimester 3 were involved in biofilm formation (K04334) and the synthesis of the antibiotic mannopeptimycin (K00815), respectively, both being contributed solely by E. coli (Figures 7 and S5; Excel file E5). Table 2 summarizes these and other key pathways and enzymes differentially abundant in T1D women.
Pathways and enzymes related to the biosynthesis and/or transport of branched chain amino acids (BCAAs; two pathways and K00826), short chain fatty acids (SCFAs) acetate, propionate and butyrate; K00171, K13788 and K00209, vitamins B1 (thiamine; K00941, K04487 and K00878), B5 (pantothenate; K00997), B6 (pyridoxine; K00868, K06215 and K08681), B7 (biotin; K16784) and B9 (folate; K13940), as well as the degradation of starch (K01208), were decreased in T1D women relative to the non-T1D (Figure 7; Excel file E5-E7) due to decreased contributions from Eubacterium rectale, E. eligens, Faecalibacterium prausnitzii, Akkermansia muciniphila, Bifidobacterium adolescentis, Bifidobacterium longum, Bifidobacterium bifidum and Lachnospiraceae bacterium 5 1 63FAA (Table 2; Figures S6-S11).
Functional differences were also present between trimesters (Excel files E8-E10). Most involved a decrease in abundance from trimester 1 to 3. Thus, in T1D women only the abundance of beta-N-acetylhexosaminidase (K01207) involved in the degradation of mucin decreased significantly from trimester 1 to 3, associated with a decrease in A. muciniphila, R. hominis, B. longum and B. adolescentis (Figures 7 and S12; Excel file E8). Enzymes involved in the synthesis of acetate (K01512), butyrate (K00248 and K00209), acetyl-CoA (K13788; a precursor of acetate and butyrate) and vitamins B1 (K00878) and B5 (K00997) also decreased from trimester 1 to 3 in T1D women only, while pyruvate ferredoxin oxidoreductase (K00171), involved in the conversion of pyruvate into acetyl-CoA, increased from trimester 1 to 3 in non-T1D women only (Figures 7 and S12-S14; Excel file E8).