Overview of preterm infant cohort

Our study examined associations between preterm infant PMA, growth, nutrition, clinical factors, and gut microbiota development in a cohort of 95 preterm and 2 full-term infants from PROP and 23 full-term infants from RPRC at the University of Rochester School of Medicine. A total of 719 rectal swab samples were collected weekly from the PROP preterm infants while in the NICU, spanning PMA from 24 to 46 weeks with a good representation across gestational ages. A total of 2 rectal swabs were collected from the 2 PROP full-term infants, and 46 rectal swabs from the 23 RPRC full-term infants: one near birth (≤ 20 day of life [DOL]) and—for the RPRC subjects only—a second at 1 month of age (20 < DOL ≤ 50) (Table 1). The longitudinal analyses included 719 samples from preterm and 48 samples from full-term infants. Relevant available age metrics were gestational age at birth, day of life (DOL), and PMA. To select the age metric for our analyses, we used functional data analysis to fit nutrition and growth variables, as well as the abundance of operational taxonomic units (OTUs) in the microbiota, first using PMA and then DOL as the age variable [21]. The overall fitting variance using PMA was lower for OTU abundance and, for most metrics of growth and nutrition, consistent with previous findings that the temporal dynamics of the preterm infant gut microbiota correspond better to PMA than to DOL [5]. Consequently, we used gestational age at birth and PMA for our analyses.

Table 1 Demographic and clinical variables Full size table

To select nutritional and clinical variables for our study, we first used an initial linear mixed effect regression that associates gestational age at birth, PMA, one main covariate (see Additional file 1: Table S1 for the full list of covariates), and its interaction with PMA, with microbiota taxa abundance. We found that mode of delivery was not significantly associated with microbiota composition after controlling for age, which is similar to recent studies where mode of delivery was not associated with the influence of breast milk in preterm infants or with stool composition in full-term infants [22, 23]. Overall, PMA has the strongest impact on microbiota composition, followed by the ratio of enteral calories, total calories normalized by body weight, proportion of dietary lipids, antimicrobial usage, proportion of dietary protein, and diuretics usage (Additional file 1: Table S1). Corticosteroids, H 2 receptor antagonists, and motility agent usage also have limited associations with some taxa. Based on these results, we applied a full linear mixed-effect regression model to identify the associations between microbiota phase, gestational age at birth, nutrient intake, and medication.

Evaluation of fecal microbiota sampling methods

To develop a sampling protocol that yields highly reproducible and representative fecal microbiota profiles, we first compared fecal microbiota obtained from matched stools or meconium and rectal swabs from five infants. Both sampling methods identify similar composition of operational taxonomic units (OTUs) and alpha diversity or evenness of observed OTUs within each subject (p values > 0.1 between stool-meconium and rectal swabs) and greater diversity between subjects (Additional file 2: Figures S1 and S2). Results from differential abundance testing on a per taxon basis between the three groups (stool, meconium, and rectal swab) are not significant. Differential abundance testing between two groups (stool-meconium and rectal swab) identified one adjusted p value < 0.1 (0.088) for the Clostridiales. In comparison, evaluation of fecal microbiota collected as matched stool and swab samples in other studies demonstrated by composition and diversity metrics that microbiota from both samples is nearly identical [24, 25]. Based on similarity of alpha diversity and the ability of clinical NICU staff to collect and store samples directly from infants at specific times without cross contamination from infant diapers and skin, we selected rectal swabs as the preferred method for sampling gut microbiota.

Three phases of the preterm infant gut microbiota

Characterization of microbiota from all subjects and time points spanning PMA from 24 to 46 weeks identified Bacilli, Gammaproteobacteria, and Clostridia as by far the most abundant taxa, with relative abundances of 41.75, 23.0, and 22.5% respectively, accounting for 87.0% of the total observed abundance (Fig. 1a). The next most abundant classes are Actinobacteria and Bacteroidia, which account for just 6.5 and 5.1% of total observed abundance, respectively. To characterize the apparent developmental phases of the premature infant gut microbiome, we used threshold values for the log ratio of the three predominant bacterial classes—Bacilli, Gammaproteobacteria, and Clostridia—to construct a decision tree that permits objective assignment of individual microbiota samples to one of the three phases based on their composition (Fig. 1a and Methods). We used criteria that distinguished the phases based on their association with prematurity and lower PMA and on the relative dominance of Bacilli, Gammaproteobacteria, and Clostridia in P1, P2, and P3 respectively. The composition of individual samples assigned to each phase at the class level using these criteria is shown in Fig. 1b. The categorical structure of the model, which assumes relative stability within a phase and abrupt shifts in composition between phases, was validated by examining the week-to-week changes of the microbiota within each subject. Quantitative changes in the weekly microbiota samples were determined using weighted UniFrac distance to measure the dissimilarity between consecutive samples. Averaged over all subjects, consecutive samples of the same phase show substantially less dissimilarity week-to-week than consecutive samples of the differing phases (i.e., samples before and after a phase transition; Fig. 1c). Testing the median dissimilarity revealed that it is significantly higher when the phase changed between consecutive samples than when it remained the same (p < 0.0001), suggesting discrete periods of community restructuring corresponding to phase transition.

Fig. 1 Overview of the preterm infant gut microbiota phases and properties. a The decision tree for classifying a microbiota sample into one of the three phases. b A composition bar chart with each sample grouped by phases 1–3 (P1–P3) from left to right. Green, gray, orange, and blue represent Bacilli, Gammaproteobacteria, Clostridia, and Bacteroidia, respectively. c Bar charts representing the average weighted UniFrac distance between consecutive samples of each individual infant. The bars are grouped into three major categories from left to right according to the initial phase of the consecutive samples being assessed. Each bar within a category corresponds to the phase of the second consecutive sample. The height of the bar indicates the average dissimilarity between consecutive samples of the corresponding phases, with exact values included in the table below the graph. d Bar charts indicating the transition probability between consecutive samples within an individual. The groupings and bars within each group indicate the phases of the first and second sample of a pair of consecutive samples, respectively, and are ordered as described in c. Transition probabilities are included in the table below the graph. e The distribution of samples over corrected gestational age in weeks. The dashed line separates the samples into early (< 34 weeks PMA) and late period (≥ 34 weeks PMA), based on functional variance of microbiota composition across all 81 individuals. f Bar charts showing the average composition of the samples in each phase at the genus level, with prominent genera labeled. For two Enterobacteriacaea and one Clostridiacaeae, the genus could not be determined and the family is indicated instead (See Additional file 2: Comment on Figure 3F). Lines connecting segments between phases indicate that the segment represents the same genus in each bar. A complete list of the genera represented here and their relative abundances can be found in Additional file 4: Table S3 Full size image

We quantified the pattern of progression of the gut microbiota with respect to the order of phase transition events using the sequence of phases observed in the consecutive samples from each individual infant to compute the transition probabilities between the phases. For each phase, the probability that the subsequent sample from the same individual will be in the same phase is higher than the probability of transitioning to a different phase (Fig. 1d). Transition from one phase to the next consecutive phase is more likely than the transition from a higher phase to a lower phase (e.g., P2 to P1) or from P1 to P3 directly. Accordingly, we found a strong relationship between PMA and microbiota phases (Fig. 1e). For this preterm cohort, 70% of all P1 samples were observed at PMA of 29 weeks or less; 84% of all P2 samples were observed from 28 weeks to 36 weeks PMA; and 78% of all P3 samples at 33 weeks or later. Eighty six percent of all samples from 37 weeks PMA and later were in P3, suggesting that preterm infant gut physiology and developmental stage influences the microbiota following birth.

For a point of comparison with our phase-based clustering, we performed Dirichlet multinomial mixture (DMM) modeling, using the composition of each sample at the class level as input. The optimal model fit was achieved with four Dirichlet components (Additional file 2: Figures S3 and S4A–B and Additional file 3: Table S2). DMM component three corresponds to our P1 cluster, component two corresponds to P2, and components one and four correspond to P3. A majority of all samples (89.9%) were classified as representing the Dirichlet component matching the phase of the sample. DMM components one and four within the P3 cluster correspond to more and less mature sub-types. DMM component four was the more mature sub-type, with the average sample occurring 2 weeks after the average component one sample, and exhibited the canonical characteristics of P3 (high Clostridia, high diversity) with little or no recognizable characteristics of P2. DMM component one was the less mature sub-type, with samples exhibiting the distinguishing characteristics of P3 while retaining to some extent features of P2 (relatively high Gammaproteobacteria). The high concordance observed between DMM components and phases provides statistically grounded support for our heuristic model (See Additional file 3: Table S2 for details).

Variance and abundance of taxa across three microbiota phases

Taxonomic analysis of all samples identified 16 phyla, 38 classes, 73 orders, 158 families, and 383 genera. Compositional differences across the phases at all taxonomic levels were characterized and pairwise comparisons were made between phases. The average composition of the samples in each phase at the genus level is shown in Fig. 1f, which represents the abundance of genera relative to bar size. The most significantly differentially abundant taxa between P2 and P3 were a variety of Clostridiales elevated in P3, including the genera Veillonella, Finegoldia, Clostridium, and Anaerococcus. The most significant differences between P1 and P3 were observed among Staphylococcaceae which were elevated in P1, and among Clostridiales genera Finegoldia and Veillonella which were elevated in the P3. A complete list of differentially abundant taxa can be found in Additional files 4 and 5: Tables S3 and S4A–C.

Functional capacity of microbiota phases

The inferred functional capacity of the microbiota was compared between the three phases, revealing differences potentially relevant to nutrient processing and microbiota-derived metabolites that contribute to establishment and maintenance of gut mucosal homeostasis (Fig. 2). P1 exhibited enrichment for bisphenol A (BPA) degradation and carotenoid synthesis pathways with BPA being an environmental contaminant frequently found in preterm infants due to repeated exposure to plastics in medical devices [26,27,28] and carotenoids conferring protection of gut microbiota against oxidative stress [29, 30]. Additional pathways were found to be significantly differentially abundant when comparisons were made between phases, including an increased capacity for synthesis of isoquinoline alkaloids, glycan and lipopolysaccharide (LPS) in P2 and P3. Protein translation, fatty acid biosynthesis and glycolysis and gluconeogenesis were increased in P1. A complete list of differentially abundant putative functions can be found in Additional file 6: Table S5.

Fig. 2 a–i Functional capacity of microbiota phases. The functional capacity of the microbiota present in each sample was inferred using PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) [63]. Each gray panel corresponds to one function and each point within a gray panel represents one sample. The samples are stratified by phase along the x axis, with red circles corresponding to P1, orange triangles corresponding to P2, and green squares corresponding to P3 samples. The sample position on the y axis indicates the relative abundance of the specified KEGG pathway, calculated as the fraction of times functional components of that pathway occurs across all organisms in the sample, with the contribution of each organism weighted by its relative abundance. Within each phase, samples are plotted on top of a box plot, which is centered on the median, with notches indicating an approximately 95% confidence interval, boxes indicating the boundaries of the first and third quartiles, and whiskers extending to the largest and smallest values no further than 1.5*(inter-quartile range) from the boxes. Points beyond the whiskers are outliers. If the notches of two boxes within the same gray panel do not overlap on the y axis, there is strong evidence that the true medians differ [69]. Functional pathways that are differentially enriched among the three phases include those that contribute to the degradation of phthalates on NICU medical devices (bisphenol degradation), protection against oxidative stress (carotenoid biosynthesis), microbiota driven increases in lipopolysaccharide (LPS) concentrations (lipopolysaccharide biosynthesis), short-chain fatty acids (fatty acid biosynthesis), isoquinoloine alkaloid biosynthesis, glycolysis and gluconeogenesis, glycan biosynthesis and metabolism, membrane transport and translation Full size image

Effect of microbiota phase on infant growth, nutrient intake, and medication

These observations prompted us to explore the potential relationship between microbiota phase, parenteral and enteral nutrient intake, medication and infant growth. Using linear mixed-effect regression models that account for subject-specific variation, with microbiota phase, gestational age at birth, nutrient intake, and medication as explanatory variables, we assessed the association of weight Z-score (standard deviation score (Z-score) of infant’s weight based on weight percentiles of a reference population matched for prematurity and sex; used as the dependent variable) and these covariates. Significant associations with weight Z-score include gestational age at birth, the phase of the microbiota, the ratio of major macronutrients, the proportion of calories administered enterally, the receipt of motility agents, antibiotics, diuretics, corticosteroids, and several interaction terms between medication/nutrition and the phase of the microbiota (Table 2). The significance with respect to weight Z-score of interaction terms between relative lipid and protein intake and the phase of the microbiota indicate that the observed association between these macronutrients and growth depends on the composition of the gut microbiota and differs between microbiota phases.

Table 2 Associations of weight Z-score with microbiome phase, nutrition, and clinical covariates Full size table

The longitudinal patterns of rectal microbiota phase transitions for 95 preterm and 25 full-term subjects are shown in Fig. 3a relative to gestational age at birth. Growth of the subjects is shown as change in weight Z-score from birth to NICU discharge for preterms and birth to 1 month for full terms. Comparison of preterm infants in P1 (N = 42; mean birth GA (gestational age) = 27.43 weeks) with those in P2 or P3 (N = 55; mean birth GA = 30.29 weeks) at the time of their first microbiota sample showed significant difference (p < 0.0001) in mean birth GA between these two groups. Furthermore, the most premature subjects (< 29 weeks birth GA) were significantly more likely to be in P1 than the full-term subjects at their first sample (61.8 vs 32.0%, p = 0.025). The change in weight Z-score is associated with length of time in each phase, with the lowest change (at the red end of the spectrum in Fig. 3a) in subjects (i.e., JE573, J5028, J1B12) who remain in phase 1 or 2 for prolonged periods. The largest negative change in weight Z-scores was associated with delays in transition to a P3 gut microbiota (p = 0.0023). Delayed achievement of P3 was also associated with prematurity, with full-term subjects reaching P3 by 1 month of age much more frequently than preterm subjects (100 vs 53.4%, p = 0.0001). Similarly, greater PMA-adjusted growth by discharge (preterms) or 1 month (full terms) was observed in the full-term subjects than in the preterm subjects (mean change in weight Z-score − 0.033 vs − 1.269, p ≈ 0.0). Eleven infants were treated for necrotizing enterocolitis (NEC) and two of these died of the disease. In those who survived, NEC was frequently followed by more than 2 weeks in P2 (J6B6F, J900B, J00F9, J2B52, and J8648). One infant who required a jejunal ostomy remained in P1 for an extended period of time (J0BE5). Thus, prolonged periods in P1 and P2 may represent the effects of lengthy antibiotic treatment and/or lack of enteral nutrition. Although the number of cases is insufficient for statistical assessment, our data suggest an association between delayed transition to P3 and a long-standing feeding intolerance in the NICU that results in administration of elemental amino acid-based formula (maroon ‘E’ in the right-hand margin of Fig. 3a).

Fig. 3 Temporal distribution of gut microbiota phases, change in infant weight and meconium clearance. a All rectal samples from 95 preterm and 25 full-term infants are plotted against post menstrual age, stratified by subject and sorted by gestational age at birth. Samples for preterm infants include those collected weekly from birth through discharge. Samples for full-term infants include the first sample after birth (collected at ≤ 20 DOL) and a second sample, collected ≤ 50 DOL. Microbiota phases (P1, red circle; P2, orange circle; P3, green circle), birth (gray diamond), stool transition (blue arrowhead), and NEC diagnosis (black square) at discharge are also indicated. Change in weight Z-score from birth to discharge, and elemental feeding requirements (maroon E) at discharge for preterm infants are indicated in the right margin. The lowest to greatest change in weight Z-score from birth to discharge spans the spectrum from red to green. In all infants, except for J94F4, the total weight change in weight Z-score from birth to discharge was negative. Weight Z-score changes in full-term infants were both positive and negative, and negative changes tended to be smaller than those observed in preterms. b Day of life of stool transition and phase transition for 38 preterm subjects in phase one (P1) at the time of their first microbiota sample. The relationship between day of life (DOL) for the initial transition out of P1 and from meconium to normal infant stool was modeled by linear regression. These results demonstrate a highly significant association between the transition out of P1 and from meconium to normal infant stool that is independent of PMA or prematurity, suggesting that the P1 and meconium microbiota are closely associated Full size image

Effect of nutrition and medication on microbiota taxa in each phase

We next examined the effects of nutrient intake and medications on the microbiota within each phase. Changes in taxa abundance would suggest adaptation of the microbiota in response to factors that the preterm infant encounters while in the NICU and may affect development of a mature, functional gut microbiota. Using a multivariate mixed-effects regression model that accounts for subject-specific variations, we assessed changes in taxa abundance at each phase, with PMA, gestational age at birth, total calories per kilogram in the past week, proportion of enteral calories in the past week, ratio of lipids, carbohydrates and protein in the past week, antibiotics and corticosteroid received during the past week, and additional medications as covariates. We found phase-specific changes in the microbiota significantly associated with the ratio of lipids, proteins, or carbohydrates in nutrition (total enteral and parenteral), with the dominant effect of all three nutrients in P3. At the phylum level, Actinobacteria and Proteobacteria are significantly associated with lipid intake, Firmicutes with protein and Actinobacteria, Proteobacteria, and Firmicutes with carbohydrates. Abundance of Bifidobacterium, an Actinobacterium most commonly linked with development and maintenance of the healthy infant gut microbiota [31, 32], is significantly associated with lipid and protein intake in P3, with increased Bifidobacterium abundance associated with increased lipid in the diet and decreased abundance with greater amounts of protein. Among the commonly used NICU medications, increased abundance of Bifidobacterium was significantly associated with use of corticosteroids and H 2 receptor antagonists in the past week in P3. A complete list of nutrition and medication variables significantly associated with changes in microbiota taxa in all the three phases can be found in in Additional file 7: Table S6.

Early and late periods of the preterm infant microbiome

To demonstrate the utility of modeling the microbiome as three compositionally defined phases, we compared gut microbiota development using two constant time periods based solely on PMA. Specifically, all preterm longitudinal samples (n = 705) were divided into two groups or periods of equal functional variance based on fitted microbiota taxa abundance, an early period (< 34 weeks PMA; n = 362) and a late period (≥ 34 weeks PMA; n = 343) (Fig. 1e and Additional file 2: Figure S5). This separation into early and late periods was used as an unbiased point of comparison to assess the utility of the phase-based approach relative to a purely temporal approach in the context of a nutrition-medication-microbiota-growth model. Using linear mixed-effects regression models as described for the phase-based nutrition analysis (Methods-Model A) with the weight Z-score as an outcome variable, we identified more significant associations in the phase based model than the period based model. The complete results of the period-based models can be found in Additional file 8: Table S7. Overall, our results demonstrate that the phase-based model of gut microbiota development in the preterm infant provides a more robust explanation of the data than the period-based model.

Potential associations of nutrition and medications with phase

After identifying the microbiota phase as a significant factor in infant growth, we sought to identify potential associations of nutrition and medications with phase by including them as explanatory variables and microbiota phases P1 and P2 as the outcome (with P3 as the baseline phase) in a multivariate mixed-effects logistic regression model. We completed separate analyses to identify associations during the early and late periods with nutrition and medications. Postmenstrual age, nutrient ratios, and the proportion of calories from enteral feeding are significantly associated with the phase of the microbiota in both periods. A higher proportion of nutritional lipids is consistently positively associated with the infant gut microbiota being in P2 and negatively associated with P1, while a higher proportion of proteins is positively associated with a P1 microbiota at an earlier PMA, negatively associated with being in P1 at later PMA, and negatively associated with being in P2 irrespective of PMA. Antibiotics are positively associated with a P2 gut microbiota, significant in the later PMA period (p = 0.0015), and nearly significant in the earlier period (p = 0.0778). The variables used in these analyses, as well as their p values and beta estimates, are provided in Table 3A–B and Additional file 9: Table S8A–B.

Table 3 Significant results of mixed effects logistic regression for nutrition and medication Full size table

Association between the meconium microbiota and transition out of phase 1

In addition to nutrition and other external factors that may influence the phase of the gut microbiome, we sought to identify dynamic aspects of host biology that correspond to phase transition. Emerging evidence suggests that the initial newborn infant gut microbiota is partially acquired by maternal transmission from the amniotic fluid and placenta before birth [10,11,12,13]. In utero, the fetus swallows large quantities of amniotic fluid that is colonized with bacteria in those mothers who deliver prematurely [11, 33]. Genera in common between amniotic fluid and the meconium, the earliest fecal material passed by infants, suggests that pioneer colonizers of the infant gut are from this maternal source. In addition to amniotic fluid that has been consumed, meconium is formed from sloughed off gastrointestinal epithelial cells which are generated as debris during periods of rapid digestive tract development and convolution of the intestinal epithelial surface. It has been established that in preterm infants, passage of meconium as stool is both delayed and prolonged and is observed well beyond the first stool, with final clearance occurring up to several weeks after birth [34].

To identify potential associations between the presence of meconium and P1 of the microbiota, we examined the relationship between clearance of the meconium from the stool and the initial transition of the microbiota out of P1 (Fig. 3b). Two infants remained in P1 (J7F5C and J8560), but did not survive beyond the first weeks in the NICU. Two infants cleared their meconium by discharge but their last microbiota sample was still in P1 (J0BE5 and J5633). The remaining 38 infants that were observed to be in P1 at their first rectal sample were included in a linear regression model using gestational age at birth and the DOL of their last P1 sample before their initial transition to another phase as explanatory variables, with the DOL of meconium clearance as the dependent variable. This model explained approximately half of the variation in the day of life of stool transition from meconium to normal infant stool (R-square = 0.51). Phase transition was found to be highly significant in this model (p < 0.0001), while gestational age at birth and the intercept did not exhibit significant associations (p values = 0.35 and 0.23, respectively), indicating that the time of stool transition was not associated with prematurity or PMA once microbiota phase transition was controlled for. On average, the last P1 sample before the initial phase transition occurred 4.7 days before stool transition was observed (Fig. 3b). To assess the similarity between the meconium and P1 microbiota samples, we first categorized all 721 samples as meconium or not, depending on whether the sample was collected from an infant that had not transitioned to normal stool, and as P1 or not, according to the decision tree. A majority of meconium samples were in P1 (59.8%) and P1 samples in meconium (64.4%) (Additional file 2: Figure S6A). We next used linear regression analysis to identify the taxa significantly associated with meconium and then again to identify those associated with P1 (Additional file 2: Figure S6B and Additional file 10: Table S9). The taxa that differ significantly between meconium and non-meconium samples are nearly identical to the taxa that differ significantly between P1 and P2-P3. These results demonstrate a highly significant association between the transition out of P1 and transition between meconium and normal infant stool, and that P1 and meconium share similar microbiota.