The gut microbiota of preterm infants develops predictably1,2,3,4,5,6,7, with pioneer species colonizing the gut after birth, followed by an ordered succession of microorganisms. The gut microbiota is vital to the health of preterm infants8,9, but the forces that shape these predictable dynamics of microbiome assembly are unknown. The environment, the host and interactions between microorganisms all potentially shape the dynamics of the microbiota, but in such a complex ecosystem, identifying the specific role of any individual factor is challenging10,11,12,13,14. Here we use multi-kingdom absolute abundance quantification, ecological modelling and experimental validation to address this challenge. We quantify the absolute dynamics of bacteria, fungi and archaea in a longitudinal cohort of 178 preterm infants. We uncover microbial blooms and extinctions, and show that there is an inverse correlation between bacterial and fungal loads in the infant gut. We infer computationally and demonstrate experimentally in vitro and in vivo that predictable assembly dynamics may be driven by directed, context-dependent interactions between specific microorganisms. Mirroring the dynamics of macroscopic ecosystems15,16,17, a late-arriving member of the microbiome, Klebsiella, exploits the pioneer microorganism, Staphylococcus, to gain a foothold within the gut. Notably, we find that interactions between different kingdoms can influence assembly, with a single fungal species—Candida albicans—inhibiting multiple dominant genera of gut bacteria. Our work reveals the centrality of simple microbe–microbe interactions in shaping host-associated microbiota, which is critical both for our understanding of microbiota ecology and for targeted microbiota interventions.
This is a preview of subscription content
Subscription info for Chinese customers
We have a dedicated website for our Chinese customers. Please go to naturechina.com to subscribe to this journal.
Get time limited or full article access on ReadCube.
All prices are NET prices.
All Illumina raw sequencing reads, including cohort samples and validation samples, have been deposited at the European Nucleotide Archive (ENA) under study accession number PRJEB36435. Source data for all figures are included in the Supplementary Tables. The public rDNA databases SILVA (119 release; https://www.arb-silva.de/documentation/release-119/) and UNITE (2017-12-01 release; https://unite.ut.ee/repository.php) were used to annotate the operational taxonomic unit (OTU) table.
Custom scripts for microbiome analyses and interaction inference are available at https://github.com/katcoyte/MK-SpikeSeq.
Charbonneau, M. R. et al. A microbial perspective of human developmental biology. Nature 535, 48–55 (2016).
Bäckhed, F. et al. Dynamics and stabilization of the human gut microbiome during the first year of life. Cell Host Microbe 17, 690–703 (2015).
Stewart, C. J. et al. Temporal development of the gut microbiome in early childhood from the TEDDY study. Nature 562, 583–588 (2018).
Derrien, M., Alvarez, A. S. & de Vos, W. M. The gut microbiota in the first decade of life. Trends Microbiol. 27, 997–1010 (2019).
Yatsunenko, T. et al. Human gut microbiome viewed across age and geography. Nature 486, 222–227 (2012).
Lim, E. S. et al. Early life dynamics of the human gut virome and bacterial microbiome in infants. Nat. Med. 21, 1228–1234 (2015).
Palmer, C., Bik, E. M., DiGiulio, D. B., Relman, D. A. & Brown, P. O. Development of the human infant intestinal microbiota. PLoS Biol. 5, e177 (2007).
Lynch, S. V. & Pedersen, O. The human intestinal microbiome in health and disease. N. Engl. J. Med. 375, 2369–2379 (2016).
Honda, K. & Littman, D. R. The microbiota in adaptive immune homeostasis and disease. Nature 535, 75–84 (2016).
Fischbach, M. A. Microbiome: focus on causation and mechanism. Cell 174, 785–790 (2018).
Widder, S. et al. Challenges in microbial ecology: building predictive understanding of community function and dynamics. ISME J. 10, 2557–2568 (2016).
Vrancken, G., Gregory, A. C., Huys, G. R. B., Faust, K. & Raes, J. Synthetic ecology of the human gut microbiota. Nat. Rev. Microbiol. 17, 754–763 (2019).
Walter, J., Armet, A. M., Finlay, B. B. & Shanahan, F. Establishing or exaggerating causality for the gut microbiome: lessons from human microbiota-associated rodents. Cell 180, 221–232 (2020).
Wolfe, B. E., Button, J. E., Santarelli, M. & Dutton, R. J. Cheese rind communities provide tractable systems for in situ and in vitro studies of microbial diversity. Cell 158, 422–433 (2014).
Connell, J. H. & Slatyer, R. O. Mechanisms of succession in natural communities and their role in community stability and organization. Am. Nat. 111, 1119–1144 (1977).
Bertness, M. D. & Callaway, R. Positive interactions in communities. Trends Ecol. Evol. 9, 191–193 (1994).
Shade, A. et al. Macroecology to unite all life, large and small. Trends Ecol. Evol. 33, 731–744 (2018).
Gregory, K. E. et al. Influence of maternal breast milk ingestion on acquisition of the intestinal microbiome in preterm infants. Microbiome 4, 68 (2016).
Gibson, M. K. et al. Developmental dynamics of the preterm infant gut microbiota and antibiotic resistome. Nat. Microbiol. 1, 16024 (2016).
DiBartolomeo, M. E. & Claud, E. C. The developing microbiome of the preterm infant. Clin. Ther. 38, 733–739 (2016).
La Rosa, P. S. et al. Patterned progression of bacterial populations in the premature infant gut. Proc. Natl Acad. Sci. USA 111, 12522–12527 (2014).
Costello, E. K., Carlisle, E. M., Bik, E. M., Morowitz, M. J. & Relman, D. A. Microbiome assembly across multiple body sites in low-birthweight infants. MBio 4, e00782-13 (2013).
Stewart, C. J. et al. Temporal bacterial and metabolic development of the preterm gut reveals specific signatures in health and disease. Microbiome 4, 67 (2016).
Pammi, M. et al. Intestinal dysbiosis in preterm infants preceding necrotizing enterocolitis: a systematic review and meta-analysis. Microbiome 5, 31 (2017).
Gasparrini, A. J. et al. Persistent metagenomic signatures of early-life hospitalization and antibiotic treatment in the infant gut microbiota and resistome. Nat. Microbiol. 4, 2285–2297 (2019).
Reynolds, L. A. & Finlay, B. B. Early life factors that affect allergy development. Nat. Rev. Immunol. 17, 518–528 (2017).
Gensollen, T., Iyer, S. S., Kasper, D. L. & Blumberg, R. S. How colonization by microbiota in early life shapes the immune system. Science 352, 539–544 (2016).
Bokulich, N. A. et al. Antibiotics, birth mode, and diet shape microbiome maturation during early life. Sci. Transl. Med. 8, 343ra82 (2016).
Shao, Y. et al. Stunted microbiota and opportunistic pathogen colonization in caesarean-section birth. Nature 574, 117–121 (2019).
Pasolli, E. et al. Extensive unexplored human microbiome diversity revealed by over 150,000 genomes from metagenomes spanning age, geography, and lifestyle. Cell 176, 649–662 (2019).
The Integrative HMP (iHMP) Research Network Consortium. The Integrative Human Microbiome Project. Nature 569, 641–648 (2019).
Nash, A. K. et al. The gut mycobiome of the Human Microbiome Project healthy cohort. Microbiome 5, 153 (2017).
Limon, J. J., Skalski, J. H. & Underhill, D. M. Commensal fungi in health and disease. Cell Host Microbe 22, 156–165 (2017).
Koskinen, K. et al. First insights into the diverse human archaeome: specific detection of Archaea in the gastrointestinal tract, lung, and nose and on skin. MBio 8, 00824-17 (2017).
Durán, P. et al. Microbial interkingdom interactions in roots promote Arabidopsis survival. Cell 175, 973–983 (2018).
Carr, A., Diener, C., Baliga, N. S. & Gibbons, S. M. Use and abuse of correlation analyses in microbial ecology. ISME J. 13, 2647–2655 (2019).
Contijoch, E. J. et al. Gut microbiota density influences host physiology and is shaped by host and microbial factors. eLife 8, e40553 (2019).
Vandeputte, D. et al. Quantitative microbiome profiling links gut community variation to microbial load. Nature 551, 507–511 (2017).
Stämmler, F. et al. Adjusting microbiome profiles for differences in microbial load by spike-in bacteria. Microbiome 4, 28 (2016).
Ishwaran, H. & Rao, J. S. Spike and slab variable selection: frequentist and Bayesian strategies. Ann. Stat. 33, 730–773 (2005).
Gonze, D., Coyte, K. Z., Lahti, L. & Faust, K. Microbial communities as dynamical systems. Curr. Opin. Microbiol. 44, 41–49 (2018).
Bucci, V. et al. MDSINE: Microbial Dynamical Systems INference Engine for microbiome time-series analyses. Genome Biol. 17, 121 (2016).
Freilich, M. A., Wieters, E., Broitman, B. R., Marquet, P. A. & Navarrete, S. A. Species co-occurrence networks: can they reveal trophic and non-trophic interactions in ecological communities? Ecology 99, 690–699 (2018).
Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey data. PLOS Comput. Biol. 8, e1002687 (2012).
Faust, K. et al. Microbial co-occurrence relationships in the human microbiome. PLOS Comput. Biol. 8, e1002606 (2012).
Watts, S. C., Ritchie, S. C., Inouye, M. & Holt, K. E. FastSpar: rapid and scalable correlation estimation for compositional data. Bioinformatics 35, 1064–1066 (2019).
Stein, R. R. et al. Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLOS Comput. Biol. 9, e1003388 (2013).
Fisher, C. K. & Mehta, P. Identifying keystone species in the human gut microbiome from metagenomic timeseries using sparse linear regression. PLoS ONE 9, e102451 (2014).
Pammi, M., Liang, R., Hicks, J., Mistretta, T. A. & Versalovic, J. Biofilm extracellular DNA enhances mixed species biofilms of Staphylococcus epidermidis and Candida albicans. BMC Microbiol. 13, 257 (2013).
We thank all of the infants and their families who participated in the study; J. Xavier, J. Ordovas-Montanes, O. Cunrath and members of the Rakoff-Nahoum laboratory for discussions; and L. Martin for assistance with sample collection. K.Z.C. is funded by a Sir Henry Wellcome Postdoctoral Research Fellowship (grant 201341/A/16/Z) and a University of Manchester Presidential Fellowship; R.S.G. is supported by grants 1R01AI153257-01 and 5R01AI139633-03; and S.R-N. is supported by a Career Award for Medical Scientists from the Burroughs Wellcome Fund, a Pew Biomedical Scholarship, a Basil O’Connor Starter Scholar Award from the March of Dimes, P30DK040561, K08AI130392-01 and a NIH Director’s New Innovator Award DP2GM136652.
C.R.M. receives grant funding from Mead Johnson Nutrition. C.R.M. also provides consulting services for Mead Johnson Nutrition, Alcresta and Fresenius Kabi, and sits on the Scientific Advisory Boards of Plakous Therapeutics and LUCA Biologics. All other authors declare no competing interests.
Peer review information Nature thanks Jeff Gore and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
a, d, A set of single-kingdom mock communities with a fixed composition of 10 bacterial (a) or 10 fungal (d) species and variable total microbial loads (indicated by the pie chart schematics underneath), were quantified using MK-SpikeSeq for relative composition (coloured bars) and absolute abundance (black or grey bars). b, e, Correlations between expected (based on initial microbial densities and known dilution factors) and MK-SpikeSeq-measured total absolute abundances of bacteria (b) and fungi (e) show that MK-SpikeSeq reliably detects absolute abundances of bacteria and fungi. For e, as exact rDNA copy numbers per fungal cell are undefined, the expected total abundances of internal transcribed spacer 1 (ITS1) are only estimates (here using 200 rDNA copies per fungal cell). c, f, Changes in absolute abundance for individual members (colour-coded as in a, d) of the bacterial (c) and fungal (f) mock communities are largely consistent with known dilution factors. g, A set of serial dilutions of a human faecal sample was quantified using MK-SpikeSeq for relative composition (coloured bars; shown are the phylum-level taxa) and absolute abundance (empty bars). h, Changes in absolute abundance for individual taxa (colour-coded in phyla as in g) across kingdoms are largely consistent with known dilution factors.
Extended Data Fig. 2 MK-SpikeSeq captures key ecological dynamics in multi-kingdom mock communities.
a, b, Two sets of defined multi-kingdom consortia, including ten bacteria and ten fungi (left, colour-coded as in Extended Data Fig. 1a, d), were assembled to model a ‘true’ (a) and a ‘false’ (b) negative interaction between one focal bacterium and one focal fungus, by varying the abundances of either these focal species or other background members. The MK-SpikeSeq quantifications of focal species highlight either consistent (a) or distinct (b) patterns between relative abundance and absolute abundance (middle). Relative abundance data may lead to a false prediction of cross-kingdom interaction between the focal species, whereas absolute abundance data measured by MK-SpikeSeq could disentangle these distinct mock ecological dynamics (right).
Extended Data Fig. 3 MK-SpikeSeq outperforms other quantification methods in cross-kingdom specificity.
A set of 40 test samples including human stools and soil samples were used to compare kingdom-specific quantifications of absolute abundance. a, MK-SpikeSeq compared with total DNA yields. Pearson correlation tests show that total DNA yields mostly only reflect bacterial community abundances. b, MK-SpikeSeq compared with flow cytometry cell enumerations using gating strategies targeted for either prokaryotes or fungi. For prokaryotic enumerations, two soil samples are highlighted owing to their high archaeal loads that cannot be distinguished from bacterial counts by flow cytometry. For fungi enumerations, results using one gating strategy are shown; attempts using two additional gating strategies show a similar overestimation of fungal counts (Supplementary Table 4). c, MK-SpikeSeq compared with kingdom-specific qPCR. Horizontal dashed lines show the limit of detection using qPCR, based on the negative control (DNA extraction of water); vertical dashed line shows the limit of detection using MK-SpikeSeq, based on a minimum of one read of non-spike-in archaeal 16S (arch16S) normalized against the average arch16S sequencing depth. Samples below the limit of detection are excluded from correlational tests. Some samples with an arch16S signal lower than the MK-SpikeSeq limit of detection showed arch16S qPCR signals that were higher than the negative control, probably owing to bacterial signals bleeding into archaea-specific qPCR. For a–c, Pearson correlation r and two-sided P values are shown (no adjustment for multiple comparisons). d, Comparison of 16S genus-level profiles sequenced with (s) or without (ns) spike-in shows largely unaltered community compositions having exogenous spike-in. e, f, Flow cytometry gating strategies used in prokaryotic (e) and fungal (f) cell counting (see b), with green showing bacterial and fungal cells and purple showing microsphere particles provided in the bacteria-counting kit. Higher voltage settings were used in flow cytometry for prokaryote cell counting than for fungi cell counting.
Extended Data Fig. 4 MK-SpikeSeq outperforms qPCR in the sensitivity of detection and robustness to sample background.
a, Comparison of sensitivity between MK-SpikeSeq and qPCR using tenfold serial dilutions of Escherichia coli and C. albicans. MK-SpikeSeq showed an increase in sensitivity of around 100–1,000-fold relative to qPCR in low-bacterial-abundance samples (detecting as few as 10 bacterial cells). For MK-SpikeSeq of E. coli samples, two levels of spike-in were used to cover the whole range of detection under the sequencing depth of around 10,000–100,000 reads per sample. For qPCR, horizontal dashed lines indicate the negative control (DNA extraction of water) and vertical dashed line shows the threshold below which pooled 16S sequencing yielded fewer than 100 reads (sequencing failed, probably owing to too low a signal). b, Comparison of robustness to host-cell background between MK-SpikeSeq and qPCR using test samples with fixed amounts of E. coli and C. albicans and a variable number of Caco-2 colonic cells. MK-SpikeSeq detected consistent (<2-fold variations) microbial abundances in samples with high host-cell background whereas qPCR under-measured microbial abundances by tenfold (ΔCt > 3.3). n = 2 for the 106 host cells group; n = 1 for the other groups.
a, In our first phase of NICU sequencing (Supplementary Information), we identified a number of samples (highlighted in red dots) that failed to yield more than 1,000 ITS1 reads per sample after quality filtering (red dashed line). Many of these samples that failed sequencing showed much lower (ΔCt > 5) ITS1 qPCR signals than the spike-in control (green dashed line), indicating poor DNA extractions of fungi in these samples. Frequency histograms of measurements are shown next to the axes. b, Reprocessing of 10 of these samples that failed sequencing led to increased ITS1 qPCR signals, indicating improved DNA extractions. c, These reprocessed samples also yielded the desired number (over 10,000) of ITS1 reads, passing our rDNA-sequencing criteria. P values by two-tailed paired Student’s t-test (b, c). d, Eight of the reprocessed samples showed non-zero fungal communities, and only two had no detectable fungal signal. The composition (coloured bars) and total abundance (empty bars) of fungal communities in these reprocessed samples are shown.
Extended Data Fig. 6 Bacterial samples cluster according to composition and infant age, but not according to infant diet, delivery mode or sex.
a, PCoA based on Bray–Curtis dissimilarities of bacterial community composition between samples (genus level). Samples are coloured by dominant taxa or are shown in white when diversity is high (IS > 4). b, PCoA with samples coloured by infant age. c, PCoA with samples coloured by infant diet. d, PCoA with samples coloured by the sex of the infant. e, PCoA with samples coloured by delivery mode. f, PCoA with samples coloured by cluster membership, calculated using DBSCAN (density-based spatial clustering of applications with noise). g, Stacked bars represent the distribution of the dominant genus within each cluster and dot plots illustrate the average day of life of samples within each cluster. A Kruskal–Wallis test with Bonferroni correction showed statistically significant differences in day of life of samples between clusters (Χ2 = 254, P << 0.0001, degrees of freedom (df) = 3). Data are mean ± s.d. h, Stacked bars indicate the proportion of genera exhibiting each type of noise per infant. Dark noise indicates increasing temporal dependence; white noise suggests that temporal dynamics are entirely random.
a–c, Scatter plots of rDNA-based total abundances of bacteria (a), fungi (b) and archaea (c) against infant day of life (DOL), measured by MK-SpikeSeq in the first phase. The red lines denote the linear regression fit and the 90% confidence bands of the best-fit line of absolute abundances on a logarithmic scale. Spearman correlations show that bacterial and archaeal—but not fungal—loads are positively associated with infant age. Samples with an undetectable kingdom-specific rDNA signal are not plotted. For archaea that show scarce signal in the cohort (c, left), a separate presence–absence plot and Χ2 test of binned samples (c, right) also show a positive correlation between archaeal loads and infant age. d, Diagnostics for linear mixed-effects model.
Extended Data Fig. 8 Fungal community composition does not map to infant age, diet, sex or delivery mode.
a, PCoA based on Bray–Curtis dissimilarities of fungal community composition between samples (genus level). Samples are coloured by dominant taxa or are shown in white when diversity is high (IS > 4). b, PCoA with samples coloured by infant age. c, PCoA with samples coloured by infant diet. d, PCoA with samples coloured by the sex of the infant. e, PCoA with samples coloured by delivery mode. f, PCoA with samples coloured by cluster membership, calculated using the DBSCAN algorithm. g, Stacked bars represent the distribution of the dominant genus within each cluster and dot plots illustrate the average day of life of samples within each cluster. A Kruskal–Wallis test with Bonferroni correction showed no statistically significant differences in day of life of samples between clusters (Χ2 = 3.06, P = 0.69, df = 5). Data are mean ± s.d. h, Stacked bars indicate the proportion of genera exhibiting each type of noise per infant. Dark noise indicates increasing temporal dependence; white noise suggests that temporal dynamics are entirely random. Notably, the mycobiomes of five infants could not be classified.
Extended Data Fig. 9 Microbe–microbe interactions are predominantly asymmetric, and inferring interactions from data that are based on relative abundance generates misleading results.
a, Heat map plotting interactions inferred by the gLV model. Each row of the heat map illustrates the effect on the target genera by other members of the gut community (left columns) or documented use of antimicrobial agents according to the clinical metadata (right columns). b, Histogram of individual antibacterial (purple) or antifungal (green) interaction strengths, split by kingdom. Antibacterial agents primarily inhibit bacteria, and antifungal agents primarily inhibit fungi; however, there is not a significant bias in the likelihood of either antimicrobial agent inhibiting their target kingdom (exact binomial tests, null hypothesis H0: P(Inhibition) = 0.5, P > 0.05). c, Stacked bar shows the proportion of different interaction types occurring between genera. Over 80% of interactions are asymmetric, being either exploitative (+/−), commensal (+/0) or amensal (−/0). d, To confirm the value of our absolute abundance methods, we inferred inter-genus interactions from relative abundance data alone using the FastSpar46 algorithm. This approach robustly identifies co-occurrence relationships between different microbial taxa in a manner that accounts for the compositional nature of relative abundance data. Notably, correlation networks cannot infer asymmetric interactions, thus this approach cannot detect the exploitation of Staphylococcus by Klebsiella. It also erroneously infers that Staphylococcus increases the growth of Candida, and cannot detect the inhibition of Klebsiella by Candida or Enterococcus. e, Steady-state relative abundances of bacteria of those subcommunities predicted to be feasible and/or linearly asymptotically stable. f, Steady-state relative abundances of fungi of those subcommunities predicted to be feasible and/or linearly asymptotically stable.
a, b, Biological replicate samples of mouse stools characterized by CFU counting of strains of interest in Fig. 3d, f were subjected to MK-SpikeSeq to determine rDNA-based absolute abundances of the specified strains (a, MK-SpikeSeq measurement for Fig. 3d; b, MK-SpikeSeq measurement for Fig. 3f). c, d, Repeat in vivo colonization experiments (c, repeat of Fig. 3d; d, repeat of Fig. 3f). Data are mean ± s.e.m.; *P < 0.05, **P < 0.01 (by two-tailed Student’s t-test). For c, n = 5 per group, t-test of K. pneumoniae CFU between the groups with or without S. epidermidis pre-colonization. For d, n = 4 for the C. albicans + K. pneumoniae group, n = 3 for the other two groups, t-test of K. pneumoniae CFU between C. albicans and C. parapsilosis groups.See Supplementary Table 15 for exact P values.
Included are Methods (SI section 1-6) and Supplementary Figures 1-18 (SI section 7). The Methods section provides details of the MK-SpikeSeq methodology, the NICU cohort, our computational analyses and in vitro and in vivo experiments. Supplementary Figures show multi-kingdom microbiome profiles of individual infants in the NICU cohort.
Included are Supplemental Tables 1-15, referenced in the Main text and SI text, providing the source data for all Figures.
About this article
Cite this article
Rao, C., Coyte, K.Z., Bainter, W. et al. Multi-kingdom ecological drivers of microbiota assembly in preterm infants. Nature 591, 633–638 (2021). https://doi.org/10.1038/s41586-021-03241-8
Nature Microbiology (2022)
The ISME Journal (2022)
Nature Reviews Microbiology (2022)
Nature Microbiology (2021)