Infant diet and maternal gestational weight gain predict early metabolic maturation of gut microbiomes

Abstract

Commensal gut bacterial communities (microbiomes) are predicted to influence human health and disease^1,2. Neonatal gut microbiomes are colonized with maternal and environmental flora and mature toward a stable composition over 2-3 years^3,4. To study pre- and postnatal determinants of infant microbiome development, we analyzed 402 fecal metagenomes from 60 infants aged 0-8 months, using longitudinal generalized linear mixed models (GLMMs). Distinct microbiome signatures correlated with breastfeeding, formula ingredients, and maternal gestational weight gain (GWG). Amino acid synthesis pathway accretion in breastfed microbiomes complemented normative breastmilk composition. Prebiotic oligosaccharides, designed to promote breastfed-like microflora⁵, predicted functional pathways distinct from breastfed infant microbiomes. Soy formula in six infants was positively associated with Lachnospiraceae and pathways suggesting a short-chain fatty acid (SCFA)-rich environment, including glycerol to 1-butanol fermentation, which is potentially dysbiotic. GWG correlated with altered carbohydrate degradation and enriched vitamin synthesis pathways. Maternal and postnatal antibiotics predicted microbiome alterations, while delivery route had no persistent effects. Domestic water source correlates suggest water may be an underappreciated determinant of microbiome acquisition. Clinically important microbial pathways with statistically significant dietary correlates included dysbiotic markers^6,7, core enterotype features⁸, and synthesis pathways for enteroprotective⁹ and immunomodulatory^10,11 metabolites, epigenetic mediators¹, and developmentally critical vitamins¹², warranting further investigation.

Figure 1:. Taxonomic Composition of Infant Fecal Microbiota

1A: Relative Abundance of Genera, Grouped by Month, Diet, and Delivery Route. Samples are grouped horizontally by month of life, diet (breastfeeding, cow’s milk formula, soy formula), and delivery route. All genera with ≥ 2% relative abundance in any sample are included, sorted vertically by phylum and relative contribution to the aggregate community of all subjects. 1B: Diversity and Major Taxa by Age and Diet. Boxplots (boxes representing interquartile ranges with median shown in black) portray alpha diversity (Shannon index) and relative abundance of Bifidobacteriaceae, Lachnospiraceae, and Enterobacteriaceae over time, separated by diet type: majority breastfeeding (N=75 samples), cow’s milk formula-feeding (N=295), and soy formula-feeding (N=32). Diversity increases with age (p<0.001) and soy (p=0.036). Bifidobacteriaceae positively correlated with breastfeeding (p=0.003), and negatively with soy (p<0.001). Lachnospiraceae increased in association with time (p<0.001) and soy (p<0.001) and decreased with breastfeeding (p=0.014). Enterobacteriaceae decreased with time (p<0.001) and GOS in cow’s milk formula (p=0.003). All p values are from multivariate longitudinal maximum-likelihood GLMMs, Tukey-corrected for multiple comparisons (Table S3). 1C: Principle Coordinate Analysis (PCoA) plot of Taxonomic Families, Colored by Major Taxa. PCoA plots of taxonomic families based on the Bray-Curtis dissimilarity index for all samples (N=402) are shaded from low (purple) to high (green) relative abundance of Bifidobacteriaceae, Lachnospiraceae, and Enterobacteriaceae, highlighting three distinct clusters. Sequential MANOVA (adonis in R, two-tailed) yielded R² values of 0.37 for Bifidobacteriaceae (p=0.001), 0.13 for Lachnospiraceae (p=0.001), 0.11 for Enterobacteriaceae (p=0.001); residual R² from a multivariate model including only these three taxa was 0.38.

Figure 2:. Dynamic Development of Amino Acid Synthesis Pathways

2A: Selected Amino Acid Synthesis Pathways, By Age and Diet (Breastfeeding vs. Formula), 2B: Selected Amino Acid Synthesis Pathways, By Age and Diet (Breastfed, Cow’s Milk Formula-Fed, Soy Formula-Fed). Scatterplots of normalized abundance (counts per million) of selected amino acid synthesis pathways vs. infant age (days) are shaded according to diet type. In plot 2A, mostly breastfeeding (N=75 samples) is compared with mostly formula feeding (N=327); in plot 2B, current majority breastfeeding (N=75), cow’s milk formula-feeding (N=295), and soy formula-feeding (N=32) are compared. Regression lines with 95% confidence interval shading are drawn using the loess method in R. All p values are two-tailed, from multivariate longitudinal maximum-likelihood GLMMs, Tukey-adjusted for multiple comparisons (Table S3). 2C: Known Reference Ranges for Human Milk Total Amino Acid (TAA) Content. Published TAA reference ranges in term breastmilk are plotted in comparison with USDA standards for infant formula to contextualize panels 2A and 2B. The line graph plots normative human milk TAA content (mg/100mL) content over time for colostrum (origin), transitional milk (0.5 months), two months, and four months post-delivery (Zhang et. al, Table 4). The bar plot shows predicted differences in total amino acid content (mg/100g total nitrogen) between USDA 2009 standards for infant formula and mature human milk (Zhang et. al, Table 8), divided by normative values for human milk; a +0.36 value for methionine indicates that formula has 36% more methionine (mg/total N) than human milk.

Figure 3:. Taxonomic and Functional Changes Associated with Soy Formula

3A: Taxonomic Structure of Soy-Exposed Infants’ Gut Microbiota. Stacked bar plots show relative abundances of taxonomic families over time from four twin pairs with at least one soy-exposed sibling. All families with ≥ 3% relative abundance are included (rare taxa aggregated as “Other”). Pre- and post-soy samples were available for three infants. Bifidobacteriaceae abundance was low pre-soy and throughout the study in all soy-exposed infants, except infant T0186_A, whose Bifidobacteriaceae recovered following soy formula cessation. Soy-discordant twin microbiomes were visibly dissimilar. 3B: Metabolic Pathways by Age and Soy Exposure Status. These boxplots (boxes representing interquartile ranges with median shown in black) show normalized abundance (normalized counts per million CPM) versus age (months) of three functional pathways: chorismate synthesis (PWY-6163), riboflavin synthesis (RIBOSYN2-PWY), and the aggregate methionine synthesis variable METCOMB. Longitudinal plots are separated into soy-naïve (N=364 current, N=359 lifetime), pre-soy (N=6), and soy-exposed samples (N=32 current, N=37 lifetime). Chorismate and riboflavin synthesis pathways increased post-soy, while methionine synthesis pathways decreased. All p values are two-tailed, from multivariate longitudinal maximum-likelihood GLMMs, Tukey-adjusted for multiple comparisons (Table S3). 3C: Changes in Chorismate Synthesis Pathway Homology Following Soy Exposure. The column graph shows numerical differences in total chorismate synthesis pathway (PWY-6163) abundance pre- and post-soy in three soy-exposed infants and one control (N=4 each timepoint). Total PWY-6163 abundance qualitatively increased more in soy-exposed infants than in the control. The heatmap plots numerical pre-post soy difference (normalized CPM) in chorismate synthesis pathway-identified genera, which qualitatively shifted towards Blautia pathway homology with soy exposure.

Figure 4:. Altered Development of Vitamin Synthesis and Carbohydrate Utilization Pathways in Association with GWG (kg)

Boxplot boxes in all panels represent interquartile ranges with the median line shown in black; total N=402 all plots. All p values in all panels are two-tailed, from multivariate longitudinal maximum-likelihood GLMMs Tukey-adjusted for multiple comparisons (Table S3). 4A: Selected Metabolic Pathways by GWG. These panels plot normalized abundance (counts per million) of the glycogen degradation pathway GLYCOCAT-PWY, the aggregate glucose degradation variable GLUCCOMB and the aggregate pyridoxine synthesis variable THISYNCOMB versus GWG (kg). The plots are colored according to GWG quartile in this population. All pathways plotted have a significant positive association with GWG. Sample size by quartile: Q1 (N=93 samples), Q2 (N=93), Q3 (N=108), Q4 (N=108). 4B: Selected Metabolic Pathways by Age, GWG. These panels plot normalized abundance (counts per million) of the same pathways by month of life, stratified by maternal GWG (kg) quartile in this population. The differences between the lowest and highest GWG become more apparent over time. 4C: Selected Metabolic Pathways by Maternal Pre-Pregnancy Body Mass Index. Normalized abundance (cpm) of the same pathways plotted versus maternal pre-pregnancy BMI, colored by GWG (kg) quartile. 4D: Selected Metabolic Pathways by Gestational Age (GA), GWG. These panels plot normalized abundance (counts per million) of the same three pathways versus GA at delivery, stratified by GWG (kg) quartile. In the lowest GWG quartile, the abundance vs GA curve slopes in the opposite direction of the abundance vs GWG curve, suggesting that lower GWG is more important with increasing GA.