Development of Inflammatory Bowel Disease Is Linked to a Longitudinal Restructuring of the Gut Metagenome in Mice.

The gut microbiome is linked to inflammatory bowel disease (IBD) severity and altered in late-stage disease. However, it is unclear how gut microbial communities change over the course of IBD development, especially in regard to function. To investigate microbiome-mediated disease mechanisms and discover early biomarkers of IBD, we conducted a longitudinal metagenomic investigation in an established mouse model of IBD, where damped transforming growth factor β (TGF-β) signaling in T cells leads to peripheral immune activation, weight loss, and severe colitis. IBD development is associated with abnormal gut microbiome temporal dynamics, including damped acquisition of functional diversity and significant differences in abundance trajectories for KEGG modules such as glycosaminoglycan degradation, cellular chemotaxis, and type III and IV secretion systems. Most differences between sick and control mice emerge when mice begin to lose weight and heightened T cell activation is detected in peripheral blood. However, levels of lipooligosaccharide transporter abundance diverge prior to immune activation, indicating that it could be a predisease indicator or microbiome-mediated disease mechanism. Taxonomic structure of the gut microbiome also significantly changes in association with IBD development, and the abundances of particular taxa, including several species of Bacteroides, correlate with immune activation. These discoveries were enabled by our use of generalized linear mixed-effects models to test for differences in longitudinal profiles between healthy and diseased mice while accounting for the distributions of taxon and gene counts in metagenomic data. These findings demonstrate that longitudinal metagenomics is useful for discovering the potential mechanisms through which the gut microbiome becomes altered in IBD. IMPORTANCE IBD patients harbor distinct microbial communities with functional capabilities different from those seen with healthy people. But is this cause or effect? Answering this question requires data on changes in gut microbial communities leading to disease onset. By performing weekly metagenomic sequencing and mixed-effects modeling on an established mouse model of IBD, we identified several functional pathways encoded by the gut microbiome that covary with host immune status. These pathways are novel early biomarkers that may either enable microbes to live inside an inflamed gut or contribute to immune activation in IBD mice. Future work will validate the potential roles of these microbial pathways in host-microbe interactions and human disease. This study was novel in its longitudinal design and focus on microbial pathways, which provided new mechanistic insights into the role of gut microbes in IBD development.

INTRODUCTION

Inflammatory bowel disease (IBD) is an increasingly prevalent chronic autoimmune disease in which the cells of the immune system attack intestinal tissue (1–3). Quality of life deteriorates, and patients die in severe cases. Unfortunately, the etiology of disease remains unclear and is likely complex (4). Discovery of the factors that contribute to IBD onset, development, and severity is needed to ensure accurate and effective health care. Epidemiological studies and animal model experiments have identified genetic factors (5–7) and lifestyle factors that associate with IBD, including diet (8) and exercise (9). But these factors are not precise predictors of disease risk, severity, or response to treatment, and many questions remain regarding disease mechanisms. Elucidating the cryptic etiology of IBD would enable new preventative measures, diagnostics, and therapies.

Recent work has implicated the gut microbiome in the development and severity of IBD (10). Individuals afflicted with Crohn’s disease or ulcerative colitis, the two principal clinical forms of IBD, harbor taxa distinct from those present in healthy controls (11–14). Shotgun metagenomics further revealed that the abundances of several microbial metabolic pathways are significantly altered in IBD guts (13, 15, 16). These associations may be causal, because gut microbes can influence the immune system and intestinal homeostasis. For example, immunosuppressive regulatory T cells (Tregs) are prevalent in the colonic lamina propria (LP) compared to other organs. However, their numbers are reduced in germfree or antibiotic-treated mice, suggesting that microbiota affect colonic differentiation of peripheral Tregs (pTregs) (17, 18). A similar loss of Tregs occurs in people with polymorphisms in IBD susceptibility genes that promote defects in Treg responses (19). Thus, gut microbes have the potential to interact with immune cells and this interaction can be altered due to host genetics and other risk factors in the development of IBD.

We hypothesized that the changes in the immune status of individuals with IBD are associated with temporal alterations in the functional capabilities of their gut microbiota. Understanding how the gut microbiome dynamically changes during IBD and how these changes relate to host symptoms and immune activation could clarify which microbiomic alterations contribute to disease onset and progression and which alterations respond to disease. We are particularly interested in elucidating specific microbial pathways that may induce or exacerbate immune activation and in distinguishing these from pathways required for survival in an inflamed intestinal environment. Addressing these issues requires a prospective, longitudinal study of the microbiome in IBD.

Longitudinal investigations of the microbiome have tended to focus on taxonomic rather than functional changes (20, 21). One study used 16S sequencing in the T-bet−/− RAG2−/− ulcerative colitis (TRUC) mouse model of inflammatory disease to identify how gut microbiome taxonomic composition changes over the course of treatment-induced remission and then investigated how microbial pathway abundances might change over time with ancestral state reconstruction techniques (22). Shotgun metagenomic sequencing provides direct insight into the functions encoded in the microbiome, but it has not been applied to a longitudinal investigation of IBD. As a result, our insight into how the gut microbiome operates dynamically in association with disease development is limited.

Mouse models of disease present an opportunity to quantify the longitudinal covariation between gut microbiome functions and IBD development while overcoming the challenges associated with a prospective human study and reducing the extensive genetic, lifestyle, and microbiome variations among humans. We implemented this approach using a well-documented IBD model (23–29), where the activity of transforming growth factor β (TGF-β) dominant-negative receptor II is driven by the CD4 promoter CD4-dnTβRII (30) (called DNR here). TGF-β is important for inducing pTreg differentiation (31), and its signaling in naive T cells results in activation and nuclear translocation of Smad2/3 molecules and regulation of target genes, including Foxp3 (32–34). Foxp3 then provides a positive-feedback loop by downregulating Smad7, thereby reducing its inhibition of TGF-β signaling (35). Absence of TGF-β signaling in T cells results in loss of Foxp3 expression and defects in the in vivo expansion and immunosuppressive capacity of pTregs (36, 37). However, excess inflammation can also potently inhibit Foxp3 induction by TGF-β (38, 39), and the presence of certain inflammatory cytokines can instead divert differentiation of Tregs into pathogenic Th17 cells (40–45). Thus, due to the involvement of TGF-β in Treg cell differentiation and the requirement for Treg-produced interleukin-10 (IL-10) to maintain intestinal homeostasis, TGF-β signaling in T cells is an important component of intestinal immunity (46–54). Furthermore, mutations in both TGF-β and IL-10 signaling pathways have been implicated in human IBD (55–58). As a result of the blockage of TFGβ signaling on their T cells and of the reduced number of pTregs, DNR animals develop spontaneous colonic inflammation and IBD that is akin to Crohn’s disease (30, 59). In addition to these physiological similarities, the DNR line serves as an effective model of human IBD because (i) human IBD is associated with mutations in SMAD3 (5, 60–62), a direct downstream target of TGF-β RII required for Foxp3 induction in the gut (33), and (ii) DNR mice model the documented effect of Smad7 overexpression in human IBD (63–65).

To obtain insight into how the longitudinal dynamics of the microbiome associate with IBD onset and progression, we followed DNR and littermate wild-type (WT) controls from weaning through severe disease. We used shotgun metagenomics to quantify how fecal microbiome structure and function change over the course of disease development in DNR mice and identified components of the microbiome that both associate with and predict immune status. We focus on longitudinal changes in biological pathways (i.e., groups of genes performing a coherent function), using estimated abundances of KEGG modules from DNR and WT metagenomes. Our work indicates that the microbiome may contain biomarkers of IBD development, clarifies mechanisms through which the microbiome may contribute to disease development, and reveals how gut microbes operate to succeed in an inflamed intestinal environment.

RESULTS

Age-matched female WT and DNR littermates were monitored longitudinally for IBD development over a period of 9 weeks, starting at 4 weeks of age, upon being weaned from their mother. As this is a T cell-mediated IBD model, we quantified peripheral CD4 and CD8 T cell activation by flow cytometry and measured the longitudinal change in the CD44hi activated fraction, which includes both effector and memory T cells (see Fig. S1 in the supplemental material). We also measured the weight of the animals over time (Fig. 1A). As expected, WT mice gained weight and maintained a constant fraction of activated T cells. DNR mice, conversely, stopped gaining weight and experienced a sharp increase in CD4 T cell activation followed by a gradual increase in CD8 T cell activation starting at 7 weeks of age (Fig. 1). These results indicate that in our facility, the DNR mice developed signs of IBD starting around week 7 and full disease by week 9. DNR mice had to be euthanized by week 15, as they had lost more than 15% of their maximum body weight. Similar to the T cell activation phenotype observed in the blood after week 7 (Fig. 1B and C), the DNR animals had a larger fraction of activated T cells in the spleen and the gut-draining mesenteric lymph node (MLN) at week 15 (Fig. S2).

FIG 1 

IBD development correlates with peripheral T cell activation in DNR mice. (A) Animal weight over time. n = 7 WT (blue) and 8 DNR (orange) mice. (B) Percentages of activated CD4 T cells among peripheral blood mononuclear cells (PBMCs). (C) Percentages of activated CD8 T cells among PBMCs. (B and C) n = 6 WT (blue) and 6 DNR (orange) mice. Error bars are standard errors of the means.

Longitudinal phenotypic monitoring of T cell activation in PBMCs of WT and DNR “bleeder” mice. Representative FACS plots of data from WT bleeder no. 89 and DNR bleeder no. 85 determined at the indicated time points are shown. The gate is set for the CD44hi population to show the percentages of activated effector and memory T cells over time. Download FIG S1, EPS file, 2 MB.

Phenotypic monitoring of CD8 T cell activation in WT and DNR “pooper” mice after disease onset. Representative FACS plots correspond to spleen or mesenteric lymph node (MLN) isolated from WT pooper no. 90 and DNR pooper no. 83 at week 15. The fractions of activated (CD44hi) and effector (CD44hi KLRG-1hi) CD8+ T cells are shown. Download FIG S2, EPS file, 1 MB.

We used shotgun metagenomics to assess how the functional potential of the gut microbiome diversifies over the course of disease progression. Specifically, we collected stool samples from parallel cohorts of DNR and WT mice weekly and performed shotgun metagenomic sequencing of samples obtained from the mice when they were 4, 5, 6, 8, 10, 12, and 13 weeks of age (see Table S1 in the supplemental material). We then quantified the abundance of KEGG modules encoded in each metagenome with ShotMAP (16), which revealed 373 modules present in at least one sample. These module abundances were then used to quantify how the within-sample diversities (alpha-diversities) of microbiome functions differed over time in DNR and WT mice. A. Kruskal-Wallis test performed to analyze the change in KEGG module Shannon entropy over time (Fig. S3) found that the DNR mice were relatively stable in their functional alpha-diversity (P = 0.47) compared to WT mice (P = 0.078). We also observed that functional alpha-diversities differed among individuals within a line over time and that this variation differed between lines in association with disease activation (week 7). Specifically, the coefficient of variation (CV) of KEGG module Shannon entropy data from a generalized linear model was higher among WT mice than DNR mice after disease activation (P = 0.0085). We also found that the CV was higher among DNR mice prior to activation, though this difference is reduced when the disproportionately variable week 5 samples are removed from the analysis (P = 0.21). These results show that the functional diversity of the mouse gut microbiome is relatively constrained early in life but increases over the lifetimes of WT but not DNR individuals.

Levels of functional alpha-diversity of the gut microbiome as measured by Shannon entropy differed over time across cohorts (left). In fact, there was significantly greater variation in the levels of Shannon entropy between the lines corresponding to the periods before and after disease onset (right). Download FIG S3, EPS file, 1.5 MB.

Project metadata. Download TABLE S1, XLSX file, 0.02 MB.

We then investigated how the composition of gut microbiome functions varies over time and between cohorts (DNR mice versus WT mice) by using an abundance-weighted beta-diversity metric (Bray-Curtis dissimilarity). At a global level, KEGG module abundances were similar in DNR and WT mice prior to week 6 but then diverged over time as IBD developed in the DNR mice (Fig. 2). Furthermore, the diversity of KEGG modules found in a metagenome was significantly associated with the week that the sample was collected within the cohort (permutational multivariate analysis of variance [PERMANOVA] P = 0.01, R2 = 0.42), as well as with the cohort’s weekly mean activated T cell status (pcCD4tCD44hi) (PERMANOVA P = 0.001, R2 = 0.16). Thus, there exist microbiome-encoded functional modules that differ in abundance in association with IBD progression in DNR mice.

FIG 2 

The taxonomic and functional diversity of the gut microbiome associates with IBD development. (A) NMDS ordination plots of the functional (left) and taxonomic (right) beta-diversity of samples from each line illustrate the significant divergence in levels of beta-diversity between lines over time. Functional beta-diversity was measured as Bray-Curtis dissimilarity based on KEGG module abundances, while taxonomic beta-diversity values represent the UniFrac distances of taxa detected in metagenomes. (B) The longitudinal variation of samples along selected NMDS dimensions similarly reveals how DNR and WT lines significantly diverge over time in terms of both their functional (left) and their taxonomic (right) beta-diversity. Smoothed (locally weighted scatterplot smoothing [LOESS]) trajectories of samples from each line over time are plotted, where gray areas represent 95% confidence intervals.

This temporal divergence in DNR versus WT microbiome functions was mirrored in the taxonomic structure of the microbiome (Fig. 2). The compositions of the gut metagenomes of the WT and DNR lines were relatively similar at early time points and began to diverge at week 6. Additionally, the microbiomes of WT mice remained relatively consistent over time compared to those of DNR mice, though they were not without temporal variation. Indeed, similarly to the functional diversity analysis, the levels of taxonomic beta-diversity of the microbiome significantly differed between the lines over time (PERMANOVA P = 0.004, R2 = 0.46), though not in a manner that corresponded to mean activated T cell status (PERMANOVA P = 0.118, R2 = 0.046). Collectively, these analyses indicate that (i) the diversities and structures of the gut microbiome differ over time in WT and DNR mice that are between 4 and 15 weeks of age; (ii) the WT and DNR microbiomes are generally consistent with each other prior to immune activation in DNR mice but diverge afterward; and (iii) immune activation is associated with changes in the subsequent successional diversification of the gut microbiome.

On the basis of these observations, we assessed how specific components of the microbiome associate with disease development. A key novelty of our approach is the use of Tweedie compound Poisson generalized linear mixed-effects models (GLMMs). These models allow us to test for differences in temporal trends in KEGG module abundance between DNR and WT mice while accounting for baseline differences between mice and genotypes, as well as for DNA extraction kit effects. GLMMs enable accurate modeling of non-normally distributed abundance data and correctly account for multiple sources of variation (66), including the intersubject variation that is present in repeated-measures designs such as the longitudinal sampling of individual mice in our study. The Tweedie compound Poisson distribution, which represents a weighted mixture of Poisson and gamma distributions, has a number of other attractive features. Its exponential relationship between variance and mean captures the overdispersion that is frequently present in environmental DNA sequence data, and its point mass at zero allows one-step fitting of zero-inflated data (versus fitting a model to determine the presence or absence of a feature before modeling nonzero components, as in hurdle models). Additionally, the Tweedie compound Poisson distribution is a continuous distribution, allowing us to use a normalized abundance measure, instead of raw counts, as the dependent variable. We provide a more detailed description of the models used in our analysis in Text S1 in the supplemental material.

Modeling methods. Download TEXT S1, DOCX file, 0.1 MB.

We first looked at overall trends of abundance trajectories for DNR mice versus WT mice as quantified by the interaction between genotype and time in the GLMM. These analyses revealed 29 KEGG modules with significant differences in abundance trends between DNR and WT mice (false-discovery rate [FDR] < 0.05). The interaction coefficient was positive for 26 of the significant modules (Table S2), which indicates that those modules became increasingly abundant in DNR mice versus WT mice over time. This set includes modules associated with UMP biosynthesis (M00051), keratin sulfate degradation (M00079), and the type III secretion system (M00332). The three modules with negative interaction coefficients, indicating decreasing abundance in DNR mice versus WT mice over time (Fig. 3), are lysine biosynthesis (M00031), lipooligosaccharide transport (M00252), and melatonin biosynthesis (M00037).

FIG 3 

Summary of GLMM results from 29 modules with significant time by group interaction. (A) The quantity plotted is the predicted marginal mean (PMM) of the slope coefficients. Significance testing was done by comparing goodness-of-fit values from full and reduced GLMM specifications, and the full model was used to produce the PMM estimates shown here. This quantity was primarily calculated to get a succinct summary of the direction of temporal change and does not always coincide with the interaction coefficient that is the focus of the main analysis. The estimates were obtained by running the lstrends function from the lsmeans R package (134). (B) The underlying KO abundance trajectories of a significant module (M00031; lysine biosynthesis) that decreases in abundance in DNR mice and increases in abundance in WT mice over time, as evidenced by a negative model slope and a positive model slope, respectively. (C) The plot was constructed as described for panel B, except that this significant module (M00330; adhesin transport) significantly increased in abundance over time in DNR mice whereas it did not change in abundance in WT mice. For both panel B and panel C, the shaded ribbons represent LOESS confidence bounds.

KEGG modules that exhibit significant group by time interaction coefficients, indicating that they differentially diversified between the two lines over time. Download TABLE S2, XLSX file, 0.1 MB.

To obtain improved temporal resolution regarding the divergence of module abundance in DNR mice, we extended our GLMMs to include a “hinge” at week 7, which is when immune activation initiates in DNR mice. This segmented regression approach has the potential to reveal modules that diverge in abundance between DNR and WT mice either between weeks 4 and 7 or between weeks 7 and 13. Only 13 of the 29 previously identified modules exhibited a significant effect using segmented regression (Fig. 4), likely due to a loss of power from partitioning the data into two smaller sets of samples. However, for these 13 modules, our results clarify when the DNR and WT abundances began to diverge (Table S3). The predominant pattern consisted of similar module abundances prior to week 7, followed by divergence after immune activation (11/13 modules). This pattern suggests that these modules respond to disease or play a role in disease progression.

FIG 4 

Modules with slopes significantly differing between groups showed primarily post-disease-onset differences in analyses performed with a segmented GLMM. For each cohort, the segmented GLMM estimate data represent two separate WT slopes (pre-week 7 and post-week 7) and two deviations from those slopes, which represent the time by group interaction that measures how DNR slopes differ from WT slopes. The estimates of these deviations are plotted, with asterisks marking coefficients that were significantly nonzero, with B-H-corrected P values of <0.2.

KEGG modules with significant interactions in the segmented GLMM analysis. Download TABLE S3, XLSX file, 0.03 MB.

Lipooligosaccharide transport (M00252), which is a two-component system with an unknown substrate in the mammalian gut, was the only module that stratified DNR and WT mice both before and after disease onset. To further investigate the potential taxa that may drive this particular signal, we assessed the taxonomic source of the KEGG sequences that recruited metagenomic reads into the module. We also quantified the distance covariance (67) between the longitudinal trajectories of the KEGG orthology groups (KOs) that comprise the module and the trajectory seen with each observed species. The results were mixed, with the former analysis suggesting the presence of primarily Streptococcus contributions, while the latter identified the greatest similarity with Lactobacillus murinus and "Candidatus Arthromitus" group trajectories (Fig. S4). The differences in the taxonomic composition of the reference data underlying these two approaches could account for these inconsistencies, as could the fact that the KEGG analysis relies on amino acid comparisons whereas the species trajectories are determined through nucleotide comparisons. Thus, an uncharacterized lipooligosaccharide transporter encoded in Streptococcus and other gut microbes decreases in abundance over time at a significantly higher rate in DNR mice than in WT mice, starting early in life before weight loss and immune activation.

Two different analyses potentially explain the taxonomic origins of the lipooligosaccharide transport KOs that were observed. (A) Species identities of KEGG ortholog sequences that recruited reads in generating the relevant KO abundances. (B) Values corresponding to distance covariance (dCov) trajectories of K09694 and K09694 and all species, with asterisks marking those dCov values that were significantly nonzero after B-H multiple-testing correction. Download FIG S4, EPS file, 2 MB.

The temporal changes seen with the type III secretion system (M00332) differed between the lines uniquely before disease onset. Specifically, the module decreased in abundance in WT mice over weeks 1 to 7, with KO K03225 primarily driving this effect. On the other hand, this module was relatively stable in DNR mice prior to disease onset, and several of the KOs that comprise the module increased in abundance in DNR mice in the later weeks (Fig. 4). The discovery of a stable, rather than decreasing, abundance of K03225 as an early indicator of IBD in DNR mice is intriguing because type III secretion systems are used by pathogens to invade the gut community and alter the gut environment (68, 69).

We next examined baseline differences in module abundance between DNR and WT mice at weaning. Early differences could result from genotype-specific selection of the gut microbiome or cage effects. Our models revealed 17 modules with significantly different intercepts (FDR < 0.05), which indicates differences in the abundances of the two lines at week 4 (Table S4). Eight of these modules, including several methanogenesis-associated pathways, had positive intercept coefficients, meaning that they were more abundant in DNR mice than in WT mice at week 4. Lipopolysaccharide biosynthesis and eight other modules showed the opposite effect and were higher in abundance in WT mice at weaning. This early-life variation in the microbiome supports hypotheses that suggest that preadolescent development of the microbiome can affect health state later in life. However, these temporal relationships are complex; later changes in abundance, as captured by the time by cohort interaction (which measures the difference in slopes between DNR and WT mice), could reverse the pattern seen at weaning.

KEGG modules with significant intercepts, indicating that they exhibited significantly different abundances between lines at the initial time point. Download TABLE S4, XLSX file, 0.04 MB.

To explore the temporal dynamics of specific gut taxa, we applied the same GLMM analysis to species abundances. This analysis yielded no significant results at an FDR of <0.05, likely due to not having the advantage of grouping components across a higher-order variable. While species could be grouped into higher taxonomic entities, the model assumption that members of the same group tend to covary across samples and over time may be violated because members of the same taxonomy may compete with or ecologically exclude one another (70). We evaluated this possibility by applying nonparametric decomposition of variance components (71) to assess whether within-module or within-genus dispersion decomposition patterns were significantly different from those obtained from random permutations of the underlying data. This auxiliary analysis found that components of functional groups covaried more than components of randomly chosen groups whereas components of taxonomic groups did not (Fig. S5). This observation indicates that grouping taxa would violate GLMM assumptions. Consequently, we instead used a goodness-of-fit test based on functional principal-component analysis (FPCA), which is less rigid in its assumption of linearity and capable of borrowing information across species due to the representation of abundance trajectories as combinations of eigenfunctions derived from the entire data set. This test identified seven species that significantly differed in their levels of variation over time between the DNR and WT cohorts (Fig. 5; Table 1), including greater increases in abundance over time within DNR microbiomes for Escherichia coli and four species from the Bacteroides genus, which are associated with gut inflammation (10).

FIG 5 

Species that showed significantly different trajectory shapes between DNR and WT groups. These results are based on an FPCA-based goodness-of-fit comparison test that identified 7 species that were different at an FDR of <0.05.

TABLE 1 

Species with significantly different trajectory shapes in the FPCA-based goodness-of-fit comparisons

Species ID	P value	FDR	Species name	WT area underLOESS curve	DNR area underLOESS curve
54642	0	0	Bacteroides sartorii	0.01992	0.07699
57185	0	0	Bacteroides xylanisolvens	0.03051	0.02506
57318	0	0	Bacteroides uniformis	0.02297	0.04506
58110	0	0	Escherichia coli O157:H43 strain T22	5.35E−4	0.007523
59684	0.0001	0.0025	Lachnospiraceae bacterium COE1	0.07203	0.04213
59708	0	0	Bacteroides rodentium	0.0136	0.01986
61442	0.0013	0.02786	Lachnospiraceae bacterium A4	0.119	0.1348

LOESS, locally weighted scatterplot smoothing.

Distributions of F-statistics computed by the DISCO method for KO and species vectors within module and genus groupings, respectively, compared to values from permuted groupings. Results of Kolmogorov-Smirnov tests show significant differences between real and permuted distributions in the functional groupings but no significant differences in the taxonomic groupings. Download FIG S5, EPS file, 1 MB.

DISCUSSION

The results of this study represent the first shotgun metagenomic characterization of IBD development. By using a controlled mouse model, a longitudinal study design, and statistical modeling, we identified novel microbial biomarkers associated with IBD onset and progression. Many of the taxa and functions that we implicated have known roles in immune regulation and pathogenicity, making them plausible candidates for stimulating the disease process, while others likely represent responses of the microbiota to changes in host physiology. Ordination and GLMM analyses enabled us to distinguish these scenarios by identifying significant differences between DNR and WT mice over time (from weaning through severe disease). We discovered that the lipooligosaccharide transport and type III secretion protein abundance trajectories that occur between weaning and immune activation differentiate DNR mice prior to immune activation, making them promising early biomarkers and consistent with a potentially causal role in IBD. Abundances of 17 modules are altered in DNR mice at weaning and could predict IBD risk if they generalize to other mouse models and human disease (see below). Many other modules as well as a few species have altered abundances in DNR mice in later, more-severe stages of disease. Functional and taxonomic diversities also show temporal differences in DNR mice that correlate with immune profiles and/or disease progression. Most of these discoveries would have been missed in a cross-sectional study because the disease association is a longitudinal trend.

By using shotgun metagenomics, we were able to investigate both taxonomic and functional characteristics of the IBD microbiome. Both types of data consistently showed differences between DNR and WT mice. For example, beta-diversity analyses revealed increasing divergence of both taxonomic and functional profiles between DNR and WT microbiomes over the last 4 weeks of the study. In addition, the individual taxa and modules with genotype-specific trajectories predominantly had increased abundance in DNR mice after disease onset. These similarities in the successional diversification of species and genes support the idea that taxonomic changes in IBD have functional consequences that are linked to immune activation. Despite such parallels, our taxonomic and functional results differed in several important ways. Notably, the number of species that stratified lines over time was lower than the number seen with KEGG modules. Furthermore, most of the IBD-associated modules that we discovered were not represented solely in singular species and would have been missed by considering information from taxonomic analyses only. These results could have been due to disease-associated functional redundancy, wherein a gene that is enriched in DNRs might derive from a different species in each mouse. Other potential reasons include (i) higher power due to grouping protein families into modules and (ii) missed taxonomic associations due to the relatively low number of laboratory mouse-associated microbes in the genome tree of life (72). Future work should explore how taxa missed by reference-based quantitation vary in association with IBD in DNR mice.

Despite finding relatively few species that distinguish DNR mice, we can gain insight into the disease process from what is known about how these taxa interact with the host. It is striking that four of the seven species that change in abundance as IBD develops belong to the genus Bacteroides and that three of them are more abundant in DNR mice. Several studies have implicated Bacteroides in intestinal inflammation. For example, the members of a subset of B. fragilis strains carry a proinflammatory metalloprotease toxin that has been identified in 19.3% of patients with active IBD (73), and the inoculation of animals with such strains is associated with severe colitis (10, 74). Subsequent research showed that multiple commensal species of Bacteroides could be incorporated into the gut microbiomes of IBD-susceptible genotypes of mice, including mice with TGF-β susceptibility loci, to induce IBD (75). Supporting the idea that Bacteroides species contribute to IBD, we observed a modest increase (FDR = 0.1898) in the hemophore/metalloprotease transport system module (M00328) in DNR mice as disease progressed. These and other mechanistic hypotheses must be tested, because the species of Bacteroides that we identified are diverse and species within the same genus can exhibit discordant patterns of interaction with host physiology (76).

Cross-sectional and mechanistic investigations of IBD support our finding that disease development is linked to microbiome taxonomy and function (4, 77, 78). The occurrence of progressive divergence of DNR and WT microbiomes as IBD worsens is consistent with a 16S-based study using a different mouse model of IBD in which gut microbes and imputed functions changed in association with disease status and therapeutically induced remission (22). Additionally, studies in germfree mouse models of IBD implicate the gut microbiome in disease development. For example, interleukin-10 (IL-10) knockout mice grown under germfree conditions do not develop colitis, whereas conventionally raised mice do (79). Similar findings have been reported for the TRUC mouse model (80). Furthermore, IL-10 knockout (81) and IL-2-deficient (82) mice manifest differential levels of severity of colitis dependent on the types of taxa that colonize their gut. Human studies of IBD have yet to investigate the disease longitudinally. However, our results are consistent with microbiome case-control studies that found significant differences between the taxonomic (11, 83–87) and functional (13, 16, 22) profiles of IBD patients and those of healthy controls, especially in cases of Crohn’s disease. Additionally, clinical administration of antibiotics shows promise for reducing the intestinal inflammation associated with IBD (88, 89). The longitudinal biomarkers that we identified are promising new candidates for investigation in the context of human disease onset and progression.

Our analyses identified several modules that implicate a pathogenic effect by the DNR microbiome. For example, DNR mouse microbiomes show increases in the abundance of adhesion protein transport modules (M00330) in association with disease, which may help pathobiotic members of the microbiome associate with and metabolize intestinal mucosa (90). Correspondingly, keratan (M00079) and dermatan (M00076) sulfate degradation pathways increase in abundance as disease progresses. Keratan sulfate and dermatan sulfate are glycosaminoglycans (GAGs) that are integral to intestinal mucosa and regulate the permeability of the gut epithelium. These sulfated GAGs are depleted in IBD patients (91), and their metabolism by intestinal bacteria, including Bacteroides thetaiotaomicron, contributes to intestinal colonization (92, 93). Furthermore, Crohn’s metagenomes exhibit an increase in abundance in GAG degradation pathways (16). DNR guts also have elevated levels of type III and type IV secretion systems, which pathogenic organisms leverage to successfully invade the gut microbiome and induce preferable ecological conditions within the gut (68, 69). Curiously, type III secretion abundance shows the opposite effect before immune activation (weeks 4 to 7), perhaps because of broad shifts in community composition after week 7 or alternatively due to microbes with type III secretion systems invading the LP and becoming less abundant in stool over time. Finally, we observed an increase in modules associated with the biosynthesis of isoprenoids, which have been linked to the stimulation of the mammalian immune system (94). Together, these DNR-associated pathways support a pathogenic role of gut microbes in IBD development. Future studies that seek to determine the existence of a microbiome-mediated etiology for IBD should consider these potential mechanisms of disease activation.

Our identification of pathways that change in association with IBD development generates many novel hypotheses about the mechanisms through which gut microbes contribute or respond to disease development. Future studies can explicitly test these hypotheses to discern the cause-and-effect relationship between the gut microbiome and inflammatory bowel disease. Several KEGG modules with different abundance dynamics in DNR mice versus WT mice appear to be associated with the microbiome’s acclimation to the disease environment. For example, we observed increases in the abundance of two-component systems (M00511, M00482) that may contribute to a cell’s ability to manage the elevated oxidative stress that exists during active IBD (95). We also observed increases in abundance in pathways associated with cellular chemotaxis (M00515, M00507). This result is consistent with observations of increased cell motility pathways in the gut microbiomes of TRUC mice suffering active colitis using imputations from 16S data (22). This result also aligns with prior work that implicated Toll-like receptor recognition of flagellar bacterial antigens in the development of intestinal inflammation (96, 97). On the basis of these observations, we speculate that, given that intestinal permeability frequently increases during IBD flare-ups, chemotaxis pathways help microbiota scavenge the metabolic resources required to survive inside an inflamed gut or to invade the host (98).

We also observed several biosynthetic modules that increased in abundance in association with IBD development. For example, the modules related to the biosynthesis of UMP, leucine, proline, and ammonia changed in association with disease. These results may suggest that the metabolic preferences and needs of the organisms that comprise the microbiome change as disease develops. Alternatively, it may be that more T cells are entering the gut, becoming activated, and consequently consuming the local resources, which in turn results in bacteria activating biosynthetic pathways to survive and compete. Our finding that pathways associated with ammonia production (M00531) increase in abundance in DNR mice is noteworthy because prior studies have found that IBD associates with a lower pH in the intestinal lumen (99), and the production of ammonia by bacteria may serve to buffer such pH changes. Additionally, these pathways are utilized when bacteria metabolize proteins, amino acids, and urea, and the increase in this pathway may indicate a preferential utilization of these substrates by the microbiome (or, as described above, immune cells) during disease.

Furthermore, we observed increases in modules associated with choline metabolism, specifically, in betaine and phosphatidylcholine biosynthesis modules. Recent work has connected the gut microbiome’s production of these metabolites to increased cardiovascular disease risk (100). Our finding is important because a growing number of studies indicate that IBD patients have an increased risk of developing cardiovascular disease, especially during flare-ups (101, 102). The mechanisms underlying this increased risk are not well resolved but may relate to a proposed explanation for the increased cardiovascular disease risk observed in HIV-infected patients (103, 104). In this model, changes in the relative proportions of protective and pathobiotic gut microbiota, especially those capable of translocating across the gut epithelium, activate a chronic systemic inflammation that increases cardiovascular disease risk. It is thus tempting to speculate that, on the basis of our observations in these mouse models of disease, IBD-associated and perhaps HIV-associated changes in the microbial metabolism of choline contribute to or at least indicate the presence of this increased risk of cardiovascular disease.

Another intriguing hypothesis emerges from our observation that levels of heme transport genes are elevated in DNR mice as IBD develops. Bacteria use this module to scavenge iron from the environment. Iron is a crucial component for many cellular processes, but gut microbes seldom have access to free iron and instead sequester it from host sources, such as heme (105, 106). Heme concentrations may be increased in IBD, as a common feature of the disease is intestinal bleeding (107). Hence, we hypothesize that gut microbes that can take advantage of this heme may flourish in DNR mice. It is intriguing to further speculate that microbial sequestration of heme contributes to IBD (e.g., through signaling to the immune system) or to iron deficiency in IBD patients (108).

One surprising discovery was an increase in the abundance of pathways associated with the production of benzoate (M00538) in DNR mice. Benzoate is a carboxylic acid produced by microbial degradation of dietary aromatic compounds and is a precursor of hippurate biosynthesis in mammals (109). Prior work suggested that hippurate may be a useful diagnostic of Crohn’s disease given that it is found at significantly lower levels in the urine of patients (109) and that the gut microbiome’s production of benzoate is responsible for these differences in urinary hippurate levels (110). Our results are inconsistent with this prior work in that they indicate that intestinal benzoate biosynthesis levels are higher in sick animals. This difference may be due to variations in the host species being investigated, including how benzoate is subsequently metabolized in the gut or by the host. Alternatively, the potential of the DNR microbiome to make excess hippurate may not be realized given that we performed DNA sequencing. Future mechanistic studies could measure benzoate and hippurate levels and quantify the benzoate proteins at the RNA or protein level in DNR mice versus WT mice.

Most of the taxonomic and functional IBD biomarkers that we identified become increasingly abundant in DNR mice throughout the disease process. But three modules, namely, melatonin biosynthesis (M00037), lysine biosynthesis (M00031), and lipooligosaccharide transport (M00252), show the opposite trajectory and decrease in abundance over time in DNR mice relative to WT mice. Melatonin has a dual effect on the immune system, acting in a stimulatory manner in early infection and in an immunomodulatory manner in cases of prolonged inflammation (111). The effects of melatonin produced by gut commensals have not been studied as extensively as those of endogenous melatonin. Traditionally, melatonin acts as a potent antioxidant, although additional quorum signaling functions in bacteria have been recently reported (112). The reduction in melatonin biosynthesis capacity observed in the DNR mice could have been caused by the expansion of species that can tolerate a highly oxidative environment (113) or by microbes that utilize other strategies for neutralizing reactive oxygen species. Without metabolite data, it is not possible to definitively say that the final concentrations of melatonin are reduced in the disease state, since the decrease can be offset by host production. With respect to lysine biosynthesis, this module is also depleted in human IBD microbiomes (13), indicating that there may exist similar mechanisms of interaction between disease context and the gut microbiome across species. Future work should empirically test the potential role of these microbiome functions in the development of IBD, especially in individuals that are genetically susceptible to the disease.

Lipooligosaccharide transport is the only module to show significant differences in abundance trajectories both preactivation and postactivation. Intriguingly, its abundance was consistently lower in DNR mice than in WT mice throughout our study, with the largest difference occurring during weeks 4 to 7, prior to immune activation and the manifestation of disease symptoms. This finding seems surprising initially, because lipooligosaccharides are the major glycolipids that are produced by mucosal Gram-negative bacteria and are known to have proinflammatory effects (114). However, the two genes (NodI and NodJ) in the lipooligosaccharide transport system are present across diverse prokaryotes, and the substrates of this two-component ABC transporter have not been characterized beyond lipochitin oligosaccharide export in rhizobial bacteria (115, 116). Determining what this system transports in the mammalian gut and how its function changes in IBD is an exciting prospect. Regardless of the mechanism, the consistent and presymptomatic depletion of lipooligosaccharide transport genes in DNR mice makes this module a promising candidate biomarker for predicting and diagnosing IBD.

We relied on a mouse model to quantify the longitudinal interaction between the gut microbiome and disease because the extensive interindividual variations in human genetics, lifestyle, microbiome composition, and disease status and severity can complicate study design, analysis, and interpretation. We used the DNR mouse model because it is relevant to our understanding of the mucosal immunological dysregulation that occurs during human IBD and, consequently, of its interaction with the gut microbiome. Indeed, we observed immune activation in the blood of the DNR mice that was consistent with what has been observed in human IBD (117). The phenotype observed in DNR mice is akin to severe Crohn’s disease, with relatively substantial immunological activation and weight loss by week 12. Interpretations of the microbiome-disease interaction in this model should take into consideration this relatively severe disease status. Alternative mouse lines may be better models for other forms of IBD. Another consideration is that we found some baseline differences in microbiome protein abundances in DNR mice at weaning that may be specific to this genetic model of IBD. Ultimately, comparisons between our results and those obtained by the integrated Human Microbiome Project (iHMP) (118), which is longitudinally evaluating the microbiome and immune status of IBD patients, will clarify the relevance of the findings produced by the DNR model to human populations. Additionally, future research should use this model and build upon our findings to clarify how TGF-β-induced differentiation and function of T cells interact with the taxonomic structure and function of the gut microbiome.

Overall, our results indicate that the development of IBD is associated with corresponding changes in the operation of the gut microbiome. Microbial taxa and KEGG module abundances vary over time and in association with immune activation. Furthermore, our results suggest that the gut microbiome may contribute to disease by activating inflammation through metabolism of mucosa and by expressing proinflammatory and downregulating anti-inflammatory metabolites. Because our study relied on the imputation of microbiome function from DNA sequences, we cannot definitively conclude that the observed differences in the microbiome’s functional profiles manifest as differences in the metabolites produced by the microbiome. Future research that applies direct measurements of microbiome function should be used to validate and expand the results presented here. Regardless, our results hold promise for our understanding of microbiome-mediated IBD disease mechanisms and the potential of using microbiome sequencing of patient stool to classify and potentially even predict disease.

MATERIALS AND METHODS

Growth of mice and microbiome sampling.

We bred two cohorts of DNR and WT littermate control animals in the Gladstone Institutes mouse facility as follows. CD4-dnTβRII (DNR) animals were crossed to the RAG1−/− background to eliminate the T cell-mediated IBD and were transferred from Yale University to Gladstone Institutes in 2010. To initiate the experiments described in this study, DNR-RAG1−/− males were bred with C57BL/6N female animals, and DNR-RAG1−/+ progeny were again crossed to C57BL/6N females to generate a combination of RAG1−/+ and RAG1+/+ DNR and WT age-matched littermate controls. Animals were given regular chow consisting of irradiated PicoLab Rodent diet 20 (LabDiet). Only female animals were used in this study. Four cohoused WT-RAG1+/+ and five cohoused DNR-RAG1+/+ littermates were followed longitudinally for 15 weeks, and fresh fecal samples were collected weekly and stored at −30°C until they were subjected to microbiome processing. All mice from both cohorts were weighed weekly. All animal experiments were conducted in accordance with guidelines set by the Institutional Animal Care and Use Committee of the University of California, San Francisco.

Immune sampling.

Tail vein blood samples were collected weekly from a parallel cohort of “bleeder” mice (n = 6 WT, n = 6 DNR) to quantify how their immune status changed over time. These were distinct individuals from the “pooper” mice cohort (same colony and time period) subjected to stool metagenomics in order to prevent repeated tail vein blood sampling from affecting the health or microbiota of the cohort of pooper mice. Specifically, ∼100 μl (2 to 3 drops) of blood from tail vein was added to 30 μl of 1× heparin (500 units/ml). A 500-μl volume of 1× ACK (ammonium-chloride-potassium) lysis buffer (Lonza) was added directly to the cells, and the mixture was incubated at room temperature for 2 to 3 min. Cells were centrifuged at 4,000 rpm for 5 min. The top layer was aspirated, and another 500 μl of 1× ACK lysis buffer was added followed by centrifugation. Cells were resuspended in fluorescence-activated cell sorter (FACS) buffer (phosphate-buffered saline [PBS]–0.5% fetal bovine serum [FBS]), and, after blocking was performed, Fc receptors with anti-CD16/CD32, single-cell suspensions were incubated with fluorescein isothiocyanate (FITC) CD4 (GK1.5), phycoerythrin (PE) CD62L (MEL14), peridinin chlorophyll protein (PerCP)-Cy5.5 CD8a (53.6.72), and allophycocyanin (APC) CD44 (IM7) mouse antibodies for 30 min at 4°C. Stained cells were washed and acquired on an Accuri C6 cytometer (BD). Blood lymphocytes were gated on CD4+ or CD8+ fractions, and percentages of activated/memory (CD44hi) cells among CD4+ and CD8+ T cells were determined using FlowJo software (Tree Star Inc.). This cohort was separated from those subjected to microbiome sampling to eliminate the effect that repeated bloodletting might have on the microbiome. At 15 weeks of age, two WT mice and three DNR mice from the group of nonbleeding animals were euthanized. Spleen and mesenteric lymph nodes were then processed into single-cell suspensions and subjected to ACK lysis and cell surface staining as described for peripheral blood mononuclear cells (PBMCs). The level of T cell activation was quantified and found to highly correlate with the blood immune status of their “bleeder” littermates (see Fig. S2 and Table S1 in the supplemental material).

Metagenome sequencing and analysis.

QIAamp DNA stool minikits (Qiagen, Valencia, CA) were used to extract DNA from stool samples collected at weeks 4, 6, 8, 10, and 12. Samples were incubated in a water bath at an elevated temperature of 95 C to increase the lysis of bacterial cells per manufacturer instructions. A MoBio PowerFecal DNA isolation kit (MoBio, Carlsbad, CA, USA) was used per manufacturer instructions to process stool samples collected at weeks 5 and 13. The kit type was accounted for in statistical modeling to adjust for any potential differences in extraction bias between the two methods.

Purified DNA was prepared for shotgun metagenomic sequencing using the Nextera XT library preparation method (Illumina, San Diego, CA, USA). Libraries were quality assessed using quantitative PCR (qPCR) and a Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and subsequently sequenced using an Illumina HiSeq 2000 sequencing system. This produced an average of 74,427,303 100-bp paired-end sequences per sample. Metagenomic reads were quality controlled using the standard operating procedure defined by the Human Microbiome Project Consortium (119) as implemented in shotcleaner (120). Briefly, reads were quality trimmed using prinseq (121) and mapped against the mouse reference genome sequence (GRCm38) using bmtagger (122). Exact-duplicate reads were collapsed, and the subsequent high-quality data were subject to taxonomic and functional annotation. Functional annotation of metagenomes was conducted using ShotMAP as described in reference 16, with Prodigal (123) to call genes and RAPsearch2 (124) to identify metagenomic homologs of the KEGG database (downloaded February 2015). Reads mapping to mammalian sequences in the KEGG database were discarded, and the subsequent data were used to quantify the abundance of each KEGG orthology group (KO) using the reads per kilobase of genome equivalents (RPKG) abundance statistic (125). Metagenomes were taxonomically annotated using MIDAS as described in reference 72.

Statistical analyses and modeling.

The functional and taxonomic similarities between metagenomic samples were assessed using nonmetric multidimensional scaling (NMDS) as implemented through the nmds function in the labdsv R package (126). Ordinations were visualized using the ordiplot function in the vegan R package (127). For the functional similarity analysis, the vegdist function from the vegan R package quantified the Bray-Curtis dissimilarity based on KEGG module abundances. The taxonomic analysis used the generalized Unifrac (128) distances (alpha = 0.5), which were obtained by using the taxonomic tree from the Living Tree Project (129) and matching the genus and species components of tree leaf labels to the corresponding components of the MIDAS species labels in our data. Assessment of the significance of the clustering of samples in these ordination plots was conducted using PERMANOVA as implemented in the Adonis function in R.

The compound Poisson generalized linear mixed-effects model implemented in the cplm package in R (130) was used to find KEGG modules with significantly different time trends between groups while controlling for static differences between the lines and DNA extraction procedures (Qiagen versus MoBio). Random intercepts and slopes for both subjects and contributing KOs were used to capture variations between subjects and between genes while focusing on the large-scale shifts over the whole collection of abundance profiles contributing to a module. As described more thoroughly in Text S1 in the supplemental material, the general computational procedure consisted of forming subsets of the data with respect to each module’s relevant KO abundances and fitting a full model that described the RPKG abundance as a function of time, group, time by group interaction, sequencing kit, and random effects of each KO and individual. We then use two reduced models, eliminating first the interaction term and then the group term, to obtain P values via likelihood ratio tests. This is one of the recommended significance testing approaches for mixed models since it avoids using the approximations for the residual degrees of freedom that would be necessary to test significance via the t-statistic (66). To limit the number of modules tested, the input data were run through the MinPath Algorithm (131) to select a parsimonious set of modules based on the KOs present. The union of the individual parsimonious sets of all of the samples was used as the final set of tested modules. The approach of testing the dynamics of an entire module by fitting a single GLMM to a set of temporal abundances of multiple genes is modeled on the time course gene set analysis (TcGSA) method of Hejblum et al. (132), with the modification of using a different response distribution (the Tweedie compound Poisson distribution). Significant modules were selected at the 0.05 FDR threshold after controlling for multiple testing via the Benjamini-Hochberg (B-H) procedure. Species time trend differences were tested with the same approach, minus the grouping of multiple trajectories. Additional details of our modeling approach can be found in Text S1. All of the code used in this analysis is available at the following URL: https://github.com/slyalina/Mouse_IBD_2017_paper_supporting_code.

To differentiate functional changes occurring prior to immune activation, we fitted a second hinge regression to the abundances of modules that were found to have a significant time by group interaction in the main GLMM analysis. This second regression placed a break point at week 7, which represents the point at which immune activation initiated (Fig. 1). This allowed for two sets of slopes (before and after disease onset) and two sets of time by group interactions (representing deviations of DNR slopes from WT before-onset and after-onset slopes).

Alterations in the species trajectory curves were additionally tested with an alternate method aimed at highlighting differences in shape rather than slope. This method was an implementation of the FPCA-based difference in goodness-of-fit approach described previously in reference 133. The permutation-based P values from this analysis were B-H corrected, and species passing the 0.05 FDR threshold were retained.

To test the hypothesis that the distribution of the between-KO/within-KO dispersion decomposition statistics for all modules was significantly different from random grouping of functional trajectories (KOs into modules) but was not significantly different from random grouping of taxonomic trajectories (species into genera), we used the DISCO (71) nonparametric test, as well as the simulated null distributions that arise when generating random groupings of KOs and species, to obtain the real distributions of the test statistics in the two scenarios. We then performed a Kolmogorov-Smirnov test to compare the true distributions with their simulated counterparts.

Accession number(s).

The metagenomic data that were generated and analyzed in this study are available in GenBank under BioProject number PRJNA397886.