Gut microbiome evolution from infancy to 8 years of age

Abstract

The human gut microbiome is most dynamic in early life. Although sweeping changes in taxonomic architecture are well described, it remains unknown how, and to what extent, individual strains colonize and persist and how selective pressures define their genomic architecture. In this study, we combined shotgun sequencing of 1,203 stool samples from 26 mothers and their twins (52 infants), sampled from childbirth to 8 years after birth, with culture-enhanced, deep short-read and long-read stool sequencing from a subset of 10 twins (20 infants) to define transmission, persistence and evolutionary trajectories of gut species from infancy to middle childhood. We constructed 3,995 strain-resolved metagenome-assembled genomes across 399 taxa, and we found that 27.4% persist within individuals. We identified 726 strains shared within families, with Bacteroidales, Oscillospiraceae and Lachnospiraceae, but not Bifidobacteriaceae, vertically transferred. Lastly, we identified weaning as a critical inflection point that accelerates bacterial mutation rates and separates functional profiles of genes accruing mutations.

Extended Data Fig. 1 |. Timeline of samples from infants and mothers in this study.

The 26 families are given identifiers between 06 and 48. Twins are distinguished by -1 and -2 following the family identifier, and mothers are distinguished with a C0 before the family identifier. Light blue fill identifies samples that underwent shallow shotgun sequencing. Black fill identifies samples that underwent ‘extended sequencing’, that is, deep shotgun sequencing, pooled long read sequencing, and culture-enhanced sequencing. Four infants 20–2, 29–1, 39–2, and 43–1 were extensively sampled for extended sequencing. Blue diamonds identify the weaning timepoint for all children.

Extended Data Fig. 2 |. Selected species prevalence over time.

Percent of infant cohort at each timepoint that carries a species within Clostridium, Eubacterium, Bacteroides, Streptococcus, Bifidobacterium, Ruminococcus, Veillonella, Roseburia, or Blautia. Measures of center are a smoothed conditional mean (local polynomial regression), with 95% CI shown in gray.

Extended Data Fig. 3 |. Maternal postpartum microbiome dynamics.

(A) Bray-Curtis dissimilarities between stool samples from different mother (“Between”) and stool samples from the same mother over time (“Within”) (n = 87 maternal stool metagenomes, P = 2.2e-16, two-tailed Wilcoxon test). Boxes represent IQR with median line, whiskers extend to 1.5xIQR. (B) Principal coordinate analysis of maternal samples plotted by Bray-Curtis dissimilarity. Samples are colored by months postpartum. Samples from the same mother are connected. (C) Bray-Curtis dissimilarity between mother-infant dyads over time. Distance of child metagenomes to either the average microbiome pro-file for the respective mother (solid orange), the most recent stool sample from the same mother captured before it (dashed orange), and unrelated mothers at the same timepoint (solid navy blue), plotted over time. Measures of center are a smoothed conditional mean (local polynomial regression), with 95% CI shown in gray.

Extended Data Fig. 4 |. Cohort distributions by antibiotic exposure and birth mode.

(A) Percent of pre-weaning months with at least one antibiotic exposure is plotted for each infant, binned by cohort. No significant differences between cohorts are observed (FDR-corrected Kruskal-Wallis test). Median lines are displayed. (B) Total infants delivered via C-section and vaginal birth are plotted and binned by cohort. No significant differences in birth mode ratio are observed (P = 0.71, Chi-square, d.f. = 2). (C) Gestational age at birth of infants, binned by cohort. No significant differences between cohorts are observed (FDR-corrected Kruskal-Wallis test).

Extended Data Fig. 5 |. Enterotype distribution and composition.

(A) Enterotype distribution across infant samples. Principal coordinate analysis plotting Bray-Curtis dissimilarity between all infant samples (n = 1,099 stool), with enterotype indicated by color. (B) Percent pre-weaning exposure to breastmilk of infant stool in adult enterotypes, segmented by microbiome development phase. Percents are defined as the number of pre-weaning months with any breastmilk feeding, divided by total pre-weaning months. Boxes represent IQR with median line, whiskers extend to 1.5xIQR, and * and *** correspond to q < 0.05 and q < 0.001, respectively (FDR-corrected one-way Mann-Whitney test).

Extended Data Fig. 6 |. Evaluating MAG assembly approaches.

Summary statistics for (A) total high- (HQ) and medium-quality (MQ) post-DAS Tool MAGs (Putative Genomes) per stool, (B) average contigs per MAG, and (C) average N50 per MAG. HQ is defined as completeness ≥ 90%, contamination ≤ 5%, strain heterogeneity ≤ 0.5%, and MQ is defined as completeness ≥ 50%, contamination ≤ 5%. Boxes represent IQR with median line, and whiskers extend to 1.5xIQR. For each panel, six approaches were considered against 10 Putative Genomes: (1) timepoint-filtered (TF) long-read meta-assemblies, (2) timepoint-specific (TS) short-read meta-assemblies, (3) metaSPAdes merged with filtered long-read meta-assemblies, (4) metaSPAdes–nanopore merged with filtered long-read meta-assemblies, (5) unfiltered long-read meta-assemblies merged with timepoint-specific metaSPAdes, and (6) OPERA-MS.

Extended Data Fig. 7 |. Putative and Reconstructed Genome capture by time.

(A) Putative Genome count per sample increases until approximately 3 YOL before stabilizing (n = 177 stools from 20 infants). (B) Timepoints 5 and 6, covering the 3 YOL and 7–8 YOL infant samples, harbor more Putative Genomes than the early-life timepoints, taken within the first 2 YOL (n = 16 infants, 6 stools per infant). (C) Putative Genome count per timepoint is stable longitudinally (n = 10 mothers 3–4 stools per mother). (D) Comparison of Putative Genome counts per individual before and after quality filtering (n = 10 mothers, 20 infants). (E) Number of quality-filtered Putative Genomes per dRep secondary cluster during MAG dereplication. (F) Comparison of MAG counts for each sample before and after dereplication. For (A) and (F), measures of center are a smoothed conditional mean (LOESS local polynomial regression), with 95% CI shown in gray. For (B) and (C), boxes represent IQR with median line, and whiskers extend to 1.5xIQR.

Extended Data Fig. 8 |. Generalized mutation rates of persisting RGs in mothers and infants.

Aggregate breadth- and genome size-adjusted popSNPs per persisting RG plotted by years since seeding with (A) linear regression or (B) local regression trendline. Shaded region represents 95% CI.

Extended Data Fig. 9 |. Mutation rates and mutated gene profiles binned by sub-cohort and period relative to weaning.

(A-C) Breadth-and genome-size adjusted popSNP counts by years since seeding (A) for all persisting RGs (n = 255 breastfed-origin persisting MAGs, 359 formula-origin persisting MAGs; P = 0.912, two-tailed Mann-Whitney test), (B) for persisting RGs seeded pre-weaning (n = 125 breastfed-origin persisting MAGs, 192 formula-origin persisting MAGs; P = 0.754, two-tailed Mann-Whitney test), or (C) for persisting RGs seeded post-weaning (n = 130 breastfed-origin persisting MAGs, 167 formula-origin persisting MAGs; P = 0.965, two-tailed Mann-Whitney test), each binned by pre-weaning diet. (D) Mutation rates per persisting RG, binned by seeding in respect to weaning. Non-mutating RGs from each cohort are excluded prior to statistical analysis (n = 225 pre-weaning persisting MAGs, 221 post-weaning persisting MAGs; P < 0.0001, two-tailed Mann-Whitney test). (E) For each RG that persists through weaning, mutation rate calculated from seeding to the last timepoint (TP) before weaning, and mutation rate calculated between the immediate timepoints flanking weaning. RGs that did not accrue any mutations pre-weaning are excluded prior to statistical analysis (n = 45 persisting MAGs, P = 0.0292, two-tailed Wilcoxon matched-pairs test). (F) Pairwise Bray-Curtis dissimilarity between pre-weaning and post-weaning distributions relative to that of each mother, binned by intra- vs. inter-family comparisons (n = 4 infants, 4 mothers; P = 0.63 same vs. different family, pre-weaning, P = 0.72 same vs. different family, post-weaning, two-tailed Kruskal-Wallis test). Boxes represent IQR with median line, and whiskers extend to 1.5xIQR.

Fig. 1 |. Microbiome dynamics from birth through middle childhood.

a, Principal coordinate analysis plotting Bray–Curtis dissimilarity among all samples (n = 1,203 stool; FDR-corrected repeated-measures PERMANOVA P = 0.001 for day of life). Cohort, months since childbirth and Shannon diversity are indicated by shape, color and background shading. b, Bray–Curtis dissimilarity between consecutive stool samples from the same infant (n = 52 individuals), timepoint-matched samples within twin pairs (n = 26 dyads), timepoint-matched samples within mother–infant dyads (n = 52 dyads) and timepoint-matched samples across unrelated children (n = 52 individuals). Measures of center are a smoothed conditional mean (LOESS local polynomial regression), with 95% confidence interval (CI) shown in gray. The black vertical line represents the average weaning timepoint, with 95% CI shown in gray. c, Species diversity by time. For each individual, the number of species per highlighted genus captured in each stool, by time, as a percent of all species ever identified for that individual for that genus is indicated. Measures of center are a smoothed conditional mean (local polynomial regression), with 95% CI shown in gray. d, R² from repeated-measures PERMANOVA between clinical metadata and community composition, segmented by pre-weaning, transitional and stable timepoints. * and · correspond to q ≤ 0.05 and 0.05 < q ≤ 0.10, respectively. Inf., infant; Mat., maternal; Inpt., inpatient; Abx, antibiotics; PC, principal component.

Fig. 2 |. Pre-weaning exposure to breastmilk impacts transition through microbiome enterotypes.

a, Children binned into ‘Breastfed’, ‘Intermediate’ or ‘Formula’ cohorts reflecting each child’s percent of total pre-weaning months exposed to breastmilk. b, Density plots of analyzed samples for each cohort over time. Gray shading represents 95% confidence interval for the weaning month of life across all infants. c, Transition model representing nine enterotypes, with average relative abundance plotted on the right. Circle size represents samples per enterotype; circle color represents the average percent of pre-weaning months that each infant was exposed to milk before sampling; edge thickness represents the number of transitions between enterotypes; and gray shading represents the weaning window (last month of liquid food, ±3 months). All samples not within the weaning window are binned by MOL. d,e, Distribution of breastfed, intermediate and formula-fed cohorts by enterotype during developmental (d) (birth through weaning, n = 352 stools) and stable (e) (>36 MOL, n = 356 stools) phases. Boxes represent IQR with median line; whiskers extend to 1.5× IQR; and *** represents q < 0.001 (one-way ANOVA with Tukey’s post hoc test). f, First post-weaning month in an adult-like enterotype, by cohort. * represents q < 0.05 (Dunn’s corrected Kruskal–Wallis test); BF, I and F correspond to breastfed, intermediate and formula cohorts, respectively. NS, not significant.

Fig. 3 |. 3,995 RGs are recovered from 10 twin pairs and their mothers.

a, Workflow schematic for the rigorous construction of RGs from deep short-read and long-read sequencing of stool and culture-enriched communities (created in BioRender). b, Distribution of the 16,324 PGs by quality status. c, Completeness, contig count, log(N50), contamination and strain heterogeneity comparisons among MQ RGs (n = 1,911), HQ RGs (n = 2,084) and NCBI reference isolate assemblies (n = 1,050). Boxes represent IQR with median line; whiskers extend to 1.5× IQR; and **** corresponds to q < 0.0001 (Dunn’s corrected Kruskal–Wallis test). het., heterogeneity.

Fig. 4 |. Strain persistence, diversity and co-occurrence from birth through middle childhood.

a, Phylogenetic tree of the 399 taxa represented by RGs. Concentric rings from innermost represent phylum, family, total unique RGs from maternal guts and total unique RGs from infant guts. Taxa that are enriched in infant (yellow) or maternal (black) RGs are denoted by asterisks. b, Breakdown of total stool samples analyzed, total RGs and total persisting RGs by maternal or infant origin. c, Total count and transient/persisting split for all taxa with ≥20 intra-infant RGs and/or ≥10 intra-infant persisting RGs. Taxa that are enriched for persisting (light brown) or transient (dark brown) RGs are denoted by asterisks. d, Median first (green) and last (red) appearances with 90% confidence intervals (CIs) plotted for RGs of the 10 taxa with greatest persisting RG counts. Taxa are sorted by phylogeny, in order of Eubacteriales (n = 3), Bifidobacteriaceae (n = 3) and Bacteroidales (n = 4). e, Strain co-occurrence rate—the percentage of all taxon-positive timepoints (TPs) with multiple strains—for all taxa plotted in c. f, Number of taxa with multiple co-occurring strains per infant and maternal stool by MOL. Measures of center are a smoothed conditional mean (LOESS local polynomial regression), with 95% CI shown in red and blue. For a and c, *, **, *** and **** denote q < 0.05, q < 0.01, q < 0.001 and q < 0.0001, respectively (Bonferroni-corrected binomial test).

Fig. 5 |. Intra-family strain-sharing events vary by phylogeny.

a, (i) distribution of intra-family representative RGs by total individuals who carry them; (ii) distribution of RGs found in multiple individuals by sharing type. b, Incidence rate of new strain-sharing events per month per twin pair. Bins represent the months compared for each rate calculation. The blue vertical line represents average weaning timepoint across all infants. c, Breakdown of sharing type for each shared RG per taxon. Taxa are displayed if they meet one or more of the following criteria: ≥7 shared RGs, ≥6 sharing events between infant–infant dyads, ≥3 sharing events between mother–infant dyads and/or ≥4 sharing events between family triads. d, First appearance of a shared strain in mothers and infants by months after birth. Vertical lines represent average first appearance for each group. ** and **** represent q < 0.01 and q < 0.0001, respectively (FDR-corrected pairwise Wilcoxon test). NS, not significant.

Fig. 6 |. The weaning period is a mutation-generating hotspot that triggers a shift in mutated gene functions.

a, Aggregate breadth-size and genome-size adjusted popSNPs per persisting RG plotted by years since seeding in the infant gut. Local regression trendlines (LOESS) are drawn for persisting RGs that colonized before and after weaning, with 95% confidence interval shown in gray. b, Mutation rates per persisting RG, binned by seeding with respect to weaning (n = 317 pre-weaning MAGs and n = 297 post-weaning MAGs; breastfed-origin and formula-origin only; P = 0.0004, two-tailed Mann–Whitney test). c, For each RG that persists through weaning, mutation rate calculated from seeding to the last timepoint (TP) before weaning and mutation rate calculated between the immediate TPs flanking weaning (n = 91 MAGs, P < 0.0001, two-tailed Wilcoxon matched-pairs test). d, Distributions of gene families representing mutated ORFs, binned by occurrence relative to weaning and diet before weaning. Ellipses were estimated using the Khachiyan algorithm at default tolerance. e, Pairwise Bray–Curtis dissimilarity between pre-weaning and post-weaning distributions relative to that of each mother (n = 4 infants and n = 4 mothers; P = 0.0002, two-tailed Mann–Whitney test). Boxes represent IQR with median line; whiskers extend to 1.5× IQR. For b and e, *** corresponds to P < 0.001 (Mann–Whitney test); for c, **** corresponds to P < 0.0001 (Wilcoxon matched-pairs signed-rank test). PCA, principal component analysis.