Report

Growth dynamics of gut microbiota in health and disease inferred from single metagenomic samples

See allHide authors and affiliations

Science  04 Sep 2015:
Vol. 349, Issue 6252, pp. 1101-1106
DOI: 10.1126/science.aac4812

Estimating bacterial growth dynamics

The pattern of sequencing read coverage of bacteria in metagenomic samples reflects the growth rate. This pattern is predictive of growth because bacterial genomes are circular, with a single origin of replication. So during growth, copies of the genome accumulate at the origin. Korem et al. use the ratio of copy number at the origin to the copy number at the terminus to detect the actively growing species in a microbiome (see the Perspective by Segre). They could spot the difference between virulent and avirulent strains, population diurnal oscillations, species that are growing in irritable bowel disease, and what happens when a host's diet changes. Results were consistent in chemostats, in mice, and in human fecal samples.

Science, this issue p. 1101; see also p. 1058

Abstract

Metagenomic sequencing increased our understanding of the role of the microbiome in health and disease, yet it only provides a snapshot of a highly dynamic ecosystem. Here, we show that the pattern of metagenomic sequencing read coverage for different microbial genomes contains a single trough and a single peak, the latter coinciding with the bacterial origin of replication. Furthermore, the ratio of sequencing coverage between the peak and trough provides a quantitative measure of a species’ growth rate. We demonstrate this in vitro and in vivo, under different growth conditions, and in complex bacterial communities. For several bacterial species, peak-to-trough coverage ratios, but not relative abundances, correlated with the manifestation of inflammatory bowel disease and type II diabetes.

Characterization of microbiome composition and function through shotgun sequencing has provided many insights into its roles in health and disease. Gene calling (1, 2), functional/pathway analysis (36), metagenomic-wide association studies (7, 8), genome assembly (9, 10), and metagenomic single-nucleotide polymorphism (SNP) detection (11) have all shown associations between microbiome configurations and susceptibility to several diseases, including obesity (4, 12), type II diabetes (7), auto-inflammatory disorders (1, 13), metabolic disease (12, 14), and cancer (15, 16). However, these approaches, which only examine a static snapshot of the microbiome at the point of collection, cannot be used to observe the highly dynamic nature of the microbiota and the differential activity of its microbial members.

Here, we asked whether microbiota growth dynamics could be probed from a single metagenomic sample by examining the pattern of sequencing read coverage across bacterial genomes. Apart from a few examples (17), most bacteria harbor a single circular chromosome, which replicates bidirectionally from a single fixed origin toward a single terminus (18) (Fig. 1A). Thus, during DNA replication, regions that have already been passed by the replication fork will have two copies, whereas the yet unreplicated regions will have a single copy.

Fig. 1 PTR accurately measures in vitro growth rates of E. coli.

(A) Sequenced reads are mapped to complete bacterial genomes, and the sequencing coverage across the genome is plotted. Each bacteria cell in a growing population (top) will be at a different stage of DNA replication, generating a coverage pattern that peaks near the known replication origin (green vertical line in graph), and thus produces a prototypical sequencing coverage pattern with a single peak and a single trough. Bacteria from a nondividing population (bottom) have a single copy of the genome, producing a flat sequencing coverage pattern across the genome. (B) (Top) Sequencing coverage patterns of a nonreplicating (left) and actively replicating (right) E. coli from two human gut metagenomic samples. (Bottom) Distribution of PTRs of E. coli across 583 different human gut metagenomic samples (3, 7, 9) (histogram) and 58 in vitro samples from four growth experiments (box plot). (C) Genome coverage plots of E. coli, measured at different times during an in vitro cell growth experiment, showing that PTRs are highest during exponential growth (time points 1 to 2.5) and lowest in lag (time point 0) and stationary (time points 3 to 7) phases. (D) PTRs (red line) correlate (R = 0.95, P < 10−4) with the growth rate (black line), measured as the derivative of the logged abundance curve [abundance is measured as optical density of the culture (OD); blue line] in the subsequent 30 min (23). N = 2 repeats. Symbols, mean; error bars, mean ± SEM.

This concept was previously used to detect the location of the replication origin in synchronized yeast colonies (19) but also holds true in an asynchronous bacterial population in which every cell may be at a different stage of replication. Summed across the population, the copy number of a DNA region will be higher the closer that region is to the replication origin and, conversely, lower the closer that region is to the terminus (20, 21). Hence, the ratio between DNA copy number near the replication origin and that near the terminus, which we term peak-to-trough ratio (PTR), should reflect the growth rate of the bacterial population. At higher growth rates, a larger fraction of cells undergo DNA replication and more active replication forks are present in each cell (22). This results in a ratio higher than 1:1 between near-origin DNA and near-terminus DNA, thereby providing a quantitative readout of the population growth rate (21).

We grew in vitro cultures of Escherichia coli (K-12 strain) and sequenced them at multiple time points during late lag phase, exponential phase, and early stationary phase (23). During stationary phase, when most of the cells in the culture are not growing and thus have a single copy of their genome, we found uniform coverage across the genome (Fig. 1, A to C). In contrast, during exponential growth, when each bacterial cell may be at a different stage of DNA replication, the coverage pattern exhibited a single trough and a single peak, and the peak coincided with the known (24) replication origin (Fig. 1, A to C).

Similar patterns to those seen in vitro were also found for E. coli in 583 publicly available (3, 7, 9) human metagenomic fecal samples (Fig. 1B). PTRs extracted from these samples varied across individuals, in the range of 1 to 2.4, resembling the 1 to 2.6 range of ratios measured in vitro (Fig. 1B). Ratios higher than 2 are indicative of multifork replication, previously documented for E. coli (18, 22).

To examine whether PTRs provide a quantitative measure of growth rate, we calculated the temporal growth rate of E. coli at different times during its growth experiment as the derivative of its abundance across time (23). PTRs were correlated with the measured growth rate, preceding it by 30 min (R = 0.95, P < 10−4) (Fig. 1D), indicating that PTR predicts the change in abundance. To determine whether PTRs accurately reflect steady-state growth rates (as opposed to temporal growth), we grew E. coli in an aerobic chemostat in which steady growth rates were controlled by changing the dilution rate of the system to induce a 16-fold range (23) and found excellent correlation (R = 0.996, P < 0.001) (fig. S1A) between the calculated PTRs and the measured growth rates. According to theoretical models (21), PTR = 2C/G, where C is the replication time (C-period), and G is the generation time, and thus 1/log2(PTR) is proportional to G/C. This transformation was correlated with measured bacterial generation time (R = 0.96, P < 0.01) (fig. S1B), confirming PTR as its proxy, even when the replication time is unknown.

The relationship between PTR and growth rate extends to other commensal strains, as we found that PTR and temporal measured growth rate were significantly correlated in similar cell growth experiments performed on Lactobacillus gasseri and Enterococcus faecalis under anaerobic conditions (L. gasseri, R = 0.74, P < 0.001; E. faecalis, R = 0.57, P < 0.05) (fig. S2, A and B). PTR also detects changes in growth rate mediated by changes in growth conditions, as PTRs and measured growth rate were correlated in additional cell growth experiments performed on E. faecalis in aerobic conditions and on E. coli in restricted growth conditions (23) (E. coli, R = 0.92, P < 0.001; E. faecalis, R = 0.78, P < 10−4) (fig. S2, C and D). L. gasseri exhibited no growth in aerobic conditions, and accordingly we observed no change in PTRs (fig. S2E). PTRs for L. gasseri and E. faecalis were significantly different between aerobic and anaerobic conditions (E. faecalis, P < 0.05; L. gasseri, P < 0.001, Mann-Whitney U test) (fig. S2E).

To examine whether PTRs can be used to detect clinically relevant changes in culture conditions, we treated an in vitro culture of early log-phase nalidixic acid-resistant Citrobacter rodentium with the bacteriostatic antibiotic erythromycin. Because bacteriostatic antibiotics halt bacterial growth, and thus indirectly inhibit replication, we postulated that erythromycin treatment would decrease PTRs. As a control, cultures were treated with nalidixic acid, with the bactericidal antibiotic kanamycin, or left untreated. Indeed, erythromycin treatment lowered the PTR compared with controls (Mann-Whitney P < 0.001 and P < 10−4 for untreated control or nalidixic acid treated control, respectively) (Fig. 2, A and B). PTR reduction under erythromycin was evident within 30 min after administration and preceded the halt in exponential growth that was detected only 60 min after erythromycin treatment (Fig. 2, A to C). During antibiotic recovery, obtained by washing the cultures after 2.5 hours and removing the antibiotics (23), PTRs increased, consistent with the rise in abundance (Mann-Whitney P < 0.001) (Fig. 2, A and B).

Fig. 2 PTR accurately measures in vitro growth rates in multiple conditions.

(A to C) Absolute abundance levels [colony-forming units per ml (CFU/ml)] (top y axis, blue) and PTR (bottom y axis, red) as a function of time (minutes) of an in vitro culture of C. rodentium treated with erythromycin (bacteriostatic antibiotic in this setting; N = 2 repeats), compared to those of (A) an untreated control culture (N = 3 repeats); (B) a culture treated with nalidixic acid, a drug to which C. rodentium is resistant (N = 3 repeats); and (C) a culture treated with kanamycin, a bactericidal drug in this setting (N = 3 repeats). Background color indicates the treatment period (dark gray, left), recovery period (gray, middle), and early stationary phase (light gray, right). The black vertical line denotes antibiotic washout. PTR changes precede changes in growth. P values are Mann-Whitney U test between abundance (top) or PTR (bottom) of the two different cultures at times 30 to 150 min (left) or times 210 to 300 min (right). (D) Bacterial abundances (CFU/ml; top) and PTR (bottom) of L. gasseri and a mixture of six additional bacterial strains that inhabit the human gut. N = 4 repeats. Symbols, mean; error bars, mean ± SEM.

Next, we grew L. gasseri within a mixture of six commensal bacterial strains (23) and observed a rise in its PTR corresponding with the rise in abundance, with PTR and temporal growth being correlated (R = 0.64, P < 0.001) (Fig. 2D).

Together, the in vitro experiments show that PTRs precede and predict changes in abundance even in the more complex setting of a mixed bacterial community, thereby establishing the link between PTRs and growth rate.

To investigate whether PTRs remained accurate predictors of bacterial activity in a disease setting, we compared the proliferative behavior of virulent and nonvirulent (tir-mutant) strains of C. rodentium, which we used to infect C57BL/6 mice previously depleted of their native microbiota by wide-spectrum antibiotic treatment (23). We compared the in vivo abundance of both strains with PTRs and found that both showed similar behaviors 1 to 5 days post-infection (p.i.), with counts steadily rising from ~104 to 105 CFU/ml at day 1 to ~108 CFU/ml by day 5 (fig. S3). However, at 6 to 9 days p.i., the virulent strain displayed significantly higher counts (Mann-Whitney P < 0.05) (fig. S3) and PTRs (Mann-Whitney P < 0.001) (Fig. 3A) than the nonvirulent strain, likely reflecting preferential mucosal adhesion and proliferation (25). Whereas PTRs of the virulent strain were higher at days 6 to 9 p.i. compared with days 1 to 5 p.i. (Mann-Whitney P < 0.001) (Fig. 3A), PTRs of the nonvirulent strain 6 to 9 days p.i. were even lower than those of in vitro cultures of C. rodentium in stationary phase (Mann-Whitney P < 0.001) (Fig. 3A).

Fig. 3 PTRs reflect the growth dynamics of in vivo microbial communities.

(A) Shown are PTRs of virulent [icc169; wild-type (WT); N = 3 mice] and nonvirulent (tir mutant; N = 3 mice) C. rodentium 1 to 5 or 6 to 9 days p.i. of C57BL/6 mice previously depleted of their native microbiota. PTRs of stationary and exponential in vitro C. rodentium cultures are shown for reference. P values are Mann-Whitney U test. See fig. S3 for the corresponding measured abundances. (B) Relative abundance (left y axis, blue) and PTRs (right y axis, red) of Parabacteroides distasonis from fecal metagenomic samples obtained approximately every 6 hours from one human individual on 4 consecutive days (26). Plotted lines are spline interpolations using the displayed data points. Time is with respect to light cycles (Zeitgeber time, horizontal axis). P value is for 24-hour oscillations (23). (C and D) Shown are standardized PTRs (top graphs, mean ± SEM) and specific PTRs for species present in the sample (bottom heat maps), belonging to two human subjects that underwent a radical dietary change. Compared are days in which only white boiled rice was consumed (gray area) and days of normal diet (white area). A global change in bacterial growth dynamics was observed between dietary regimens (**Mann-Whitney P < 0.005, ***P < 0.001).

To explore the utility of the PTR measure within the complex metagenomic setting, we devised a computational pipeline that extracts PTRs for multiple samples within large metagenomic cohorts (supplementary text and fig. S4). In devising this pipeline, we took care to address (i) genomic differences between strains of a certain species; (ii) copy number variation of different genomic regions, and (iii) variable coverage levels stemming from sequencing depth. We show that our method is robust to these drivers of noise (fig. S5 to S8 and supplementary text), attributed to our examination of coverage across the entire genome, as opposed to comparing the coverage of origin and terminus regions directly.

Examining the full length of the genome also allows us to predict the replication origin location in different bacteria. We verified that coverage peak locations coincided with known locations of origins of replication. To this end, we applied our pipeline to 759 metagenomic stool samples from Chinese and European cohorts (7, 9) and predicted the location of the origin and terminus of replication for 187 different microbial strains. Indeed, these predictions, computed solely based on our analysis of the bacterial genome coverage patterns, agreed with the known replication origins of 132 different strains (24) (R2 = 0.98, P < 10−30) (fig. S8), and for 55 strains whose replication origin location is unknown our method generated novel predictions (fig. S9).

To determine whether PTRs can uncover a possible interplay between host genetics and gut microbiome growth dynamics we collected fecal samples from three mouse strains (Swiss Webster, BALB/c, and C57BL/6) (23) grown under identical environmental conditions. Microbiota growth dynamics, as estimated by PTRs, differed significantly across the different mouse strains. In BALB/c mice, PTRs were lower overall compared with C57BL/6 and Swiss Webster mice (P < 0.05) (fig. S10A). The reduction in growth in the BALB/c mice was driven by Parabacteroides distasonis, which displayed consistently lower PTR than in other mouse strains (fig. S10, B to D), indicating that host genetics may affect the growth dynamics of this bacterium.

Another example of the utility of PTRs in assessing physiological microbiome growth patterns is provided in microbiome diurnal oscillation patterns, which we recently linked to host susceptibility to obesity and glucose intolerance (26). Examining PTRs of fecal microbiomes collected from a human volunteer every 6 hours for four consecutive days, we identified two species, out of four that passed our PTR pipeline filters, that showed abundance levels cycling with a 24-hour periodicity (23). For both species, the PTRs also exhibited 24-hour oscillatory patterns (P < 0.05) (23) (Fig. 3B and fig. S11), suggesting that diurnal changes in the abundance of some bacteria were reflected in their PTRs.

To investigate the effect of an extreme dietary change on bacterial growth rates, two healthy human volunteers underwent an acute dietary change, in which they shifted their normal diet to one that contained only boiled white rice for 1 week, after which they reverted back to their regular dietary habits. In both participants, we observed a global change in gut bacterial growth dynamics between dietary regimens, as reflected in statistically significant differences in the PTRs of all bacteria between the days in which rice was consumed and the days in which the participant’s regular diet was followed (P < 0.005 for each participant) (Fig. 3, C and D).

We also examined body site–specific microbial growth rates in a metagenomic cohort (3) and found significantly higher PTRs in 229 tongue dorsum (1.36 ± 0.006) and 193 buccal mucosa (1.37 ± 0.007) samples, as compared with 325 stool samples (1.16 ± 0.002; P < 10–50 in both cases) (fig. S12A). Six species were present in both the oral cavity and stool with sufficient coverage to calculate PTRs, with three featuring significantly different PTRs between sites [false discovery rate (FDR)–corrected Mann-Whitney P < 0.1] (fig. S12, B to G), indicating that intersite differences stem not only from distinct bacterial compositions but also from site-specific differences in growth dynamics.

Overall, these results provide examples of functional insights that are not achievable using traditional metagenomics analysis methods and indicate that microbiome growth dynamics vary across diverse physiological conditions and locations.

To determine whether bacterial PTRs are associated with disease and different clinical parameters, we generated PTRs for every species in samples from European (N = 396) (9) and Chinese (N = 363) (7) cohorts. In both data sets, we found large variation in PTRs across samples (Fig. 4). Notably, we found statistically significant associations between the PTRs of 20 different bacteria and multiple clinical parameters, including significant correlations between the PTR of Bifidobacterium longum and occurrence of Crohn’s disease in the Spanish nationals of the European cohort (9) (FDR-corrected Mann-Whitney P < 0.005), and between the PTRs of 12 different bacteria and the occurrence of type II diabetes in the Chinese cohort (Fig. 4) (7). We also found significant correlations between PTRs and the occurrence of ulcerative colitis, body-mass index, the fraction of glycated hemoglobin (HbA1c%, a common marker of long-term glycemic control) (27), fasting serum insulin, and fasting blood glucose levels (Fig. 4).

Fig. 4 Bacterial dynamics correlate with several diseases and metabolic disorders.

PTRs of species from Chinese (7) [N = 363 samples (Q)] and European (9) [N = 396 samples (M)] cohorts are shown (box plots, left; red, median; boundaries, 25th to 75th percentiles) if their relative abundances or PTRs were significantly associated with clinical parameters. Shown are phylum membership; the number of samples for which PTRs were calculated; and a row with colored entries for each statistically significant (FDR-corrected P < 0.05) association between clinical parameters and its PTR (left column block) or relative abundance (right column block). Mann-Whitney U test and Spearman correlations were used for binary and continuous clinical parameters, respectively. (Top block) Species with significant associations between PTR and clinical parameters. (Bottom block) Species with significant associations only between relative abundance and clinical parameters. A, Actinobacteria; B, Bacteroidetes; F, Firmicutes; P, Proteobacteria; V, Verrucomicrobia; BMI, body-mass index.

These associations are independent of—and unobtainable by examining—bacterial abundances, because (i) in correlating PTRs with clinical parameters, we only used samples in which that bacteria was present, thereby withholding information about the presence or absence of the examined bacteria (23); (ii) in only 5 of the 38 statistically significant correlations were the abundance levels of the species also correlated with the same clinical parameter; and (iii) 36 of the 38 significant associations of PTR remained significant after correcting them for relative abundance levels. The PTRs of some species were correlated with clinical parameters only after correction for relative abundance, including Eubacterium rectale and the occurrence of Crohn’s disease (FDR-corrected Mann-Whitney P < 10–4).

As a global measure of the growth dynamics of the entire microbiota, for every sample we calculated both the mean and the median of the PTRs of all the bacteria present. This global measure correlated with fasting blood glucose and HbA1c% levels and with the occurrence of Crohn’s disease and type II diabetes, indicating that global microbiome growth dynamics also associate with disease (Fig. 4).

A preliminary analysis of 40 samples from the Prospective Registry in Inflammatory Bowel Disease (IBD) Study at Massachusetts General Hospital (MGH) (PRISM) cohort (23) showed that only four bacteria passed our stringent pipeline filters for PTR calculation in more than half of the samples. Notwithstanding, Eggerthella lenta presented significantly different PTRs between patients with active Crohn’s disease and patients in remission (FDR-corrected Mann-Whitney P < 0.1). Neither the abundance of E. lenta nor of the other three species differed between active and quiescent Crohn’s patients, highlighting the fact that PTRs reflect an independent feature of the effect of the gut microbiome on its host.

Overall, we present a new type of metagenomic data analysis that provides an accurate quantitative estimate of the growth dynamics of the microbiota from a single snapshot sample. These estimates have clinical relevance and correspond to changes in absolute abundances, which are masked by and unobtainable through relative abundances.

Using PTRs to “fish out” microbial kinetic behavior in a complex microbiome population could extend our understanding of how flexibly the microbiota responds functionally to environmental signals. We may be able to identify active “driver” and “modulator” species, distinguish them from bystander commensal species, and pinpoint disease-causing or disease-modulating microbes that contribute to multifactorial diseases whose activities may be masked by other bacteria. Furthermore, our method may be able to detect, follow, and assess therapeutic responsiveness of pathogenic or probiotic species introduced into the microbiome.

Supplementary Materials

www.sciencemag.org/content/349/6252/1101/suppl/DC1

Supplementary Text

Materials and Methods

Figs. S1 to S13

Table S1

References (2837)

References and Notes

  1. See the supplementary materials on Science Online for details.
  2. Acknowledgments: We thank the members of the Segal and Elinav laboratories for fruitful discussions. T.K. and D.Z. are supported by the Ministry of Science, Technology, and Space, Israel. T.K. is the recipient of the Strauss Institute fellowship for nutritional research. C.A.T. is the recipient of a Boehringer Ingelheim Fonds Ph.D. Fellowship. E.E. is supported by Yael and Rami Ungar, Israel; Leona M. and Harry B. Helmsley Charitable Trust; the Gurwin Family Fund for Scientific Research; Crown Endowment Fund for Immunological Research; estate of Jack Gitlitz; estate of Lydia Hershkovich; the Benoziyo Endowment Fund for the Advancement of Science; Adelis Foundation; John L. and Vera Schwartz, Pacific Palisades; Alan Markovitz, Canada; Cynthia Adelson, Canada; CNRS (Centre National de la Recherche Scientifique); estate of Samuel and Alwyn J. Weber; Mr. and Mrs. Donald L. Schwarz, Sherman Oaks; grants funded by the European Research Council; the German-Israel Binational foundation; the Israel Science Foundation (ISF); the Minerva Foundation; and the Rising Tide foundation. E.E. is the incumbent of the Rina Gudinski Career Development Chair. This work was supported by grants from the European Research Council (ERC) and the ISF to E.S. The European Nucleotide Archive (ENA) accession number for the microbial shotgun sequences is PRJEB9718. The pipeline tool will be made publicly available at http://genie.weizmann.ac.il/software/bac_growth.html. All human studies were approved by the Tel Aviv Sourasky Medical Center Institutional Review Board, approval numbers TLV-0658-12, TLV-0050-13, and TLV-0522-10, and the Weizmann Institute of Science Bioethics and Embryonic Stem Cell Research oversight (ESCRO) committee. The trial was reported to https://clinicaltrials.gov, identifier NCT01892956. The study did not necessitate or involve randomization. Author contributions: T.K. and D.Z. conceived the project, designed and conducted the analyses, interpreted the results, and wrote the manuscript. J.S. devised and performed all experiments in mice, anaerobic, and pathogenic bacteria. A.W. directed some of the experiments and the production of libraries and sequencing from samples collected in this study and provided critical insights to the manuscript. T.K., D.Z., J.S., and A.W. equally contributed to this work. T.A.-S., M.P.-L., and E.M. performed experiments and extraction and sequencing of samples. M.P.-F. performed experiments. G.J. performed chemostat experiments. A.H. supervised germ-free mouse experiments. N.C. aided in the design of the computational pipeline and in conducting the analyses. A.S.-M. and R.J.X. directed the PRISM cohort. C.A.T. and E.E. directed the circadian experiments. R.S. and E.S. directed the dietary intervention study. E.E. and E.S. conceived and directed the project and analyses, designed experiments, interpreted the results, and wrote the manuscript.
View Abstract

Navigate This Article