Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases

Abstract

Inflammatory bowel diseases, which include Crohn's disease and ulcerative colitis, affect several million individuals worldwide. Crohn's disease and ulcerative colitis are complex diseases that are heterogeneous at the clinical, immunological, molecular, genetic, and microbial levels. Individual contributing factors have been the focus of extensive research. As part of the Integrative Human Microbiome Project (HMP2 or iHMP), we followed 132 subjects for one year each to generate integrated longitudinal molecular profiles of host and microbial activity during disease (up to 24 time points each; in total 2,965 stool, biopsy, and blood specimens). Here we present the results, which provide a comprehensive view of functional dysbiosis in the gut microbiome during inflammatory bowel disease activity. We demonstrate a characteristic increase in facultative anaerobes at the expense of obligate anaerobes, as well as molecular disruptions in microbial transcription (for example, among clostridia), metabolite pools (acylcarnitines, bile acids, and short-chain fatty acids), and levels of antibodies in host serum. Periods of disease activity were also marked by increases in temporal variability, with characteristic taxonomic, functional, and biochemical shifts. Finally, integrative analysis identified microbial, biochemical, and host factors central to this dysregulation. The study's infrastructure resources, results, and data, which are available through the Inflammatory Bowel Disease Multi'omics Database ( http://ibdmdb.org ), provide the most comprehensive description to date of host and microbial activities in inflammatory bowel diseases.

Conflict of interest statement

J. Braun is on the Scientific Advisory Board for Janssen Research & Development, LLC. C.H. is on the Scientific Advisory Board for Seres Therapeutics. J.F.P. and N.J.A. own shares at Diversigen Inc. R.J.X. is a consultant to Novartis and Nestle.

Figures

Fig. 1. Multi-omics of the IBD microbiome in the IBDMDB study.

a, Overview of cohort characteristics. We followed 132 participants (with CD, with UC, or without IBD (control)) for one year each. Principal component analysis (PCA) of SNP profiles shows that the resulting IBDMDB cohort is mostly of European ancestry as compared to the 1000 Genomes (1kG) reference (see Methods). b, Sampling strategy. The study yielded host and microbial data from colon biopsy (baseline), blood (approximately quarterly), and stool (every two weeks), assessing global time points for all subjects and dense time courses for a subset. Raw, non-quality-controlled sample counts are shown. c, Overlap of multi-omic measurements from the same sample (strict) or from near-concordant time points (with differences of up to 2 or 4 weeks; see Methods). d, Principal coordinates analysis (PCoA) based on species-level Bray–Curtis dissimilarity; most variation is driven by a tradeoff between phylum Bacteroidetes versus Firmicutes. Samples from individuals with IBD (CD in particular) had weakly lower Gini–Simpson alpha diversity (Wald test P = 0.26 and 0.014 for UC and CD compared with non-IBD, respectively). e, Mantel tests quantifying variance explained (square of Mantel statistic) between measurement type pairs, with differences across subjects (inter-individual) or within subjects over time (intra-individual; see Methods); results show tight coupling across measurement types. Sample sizes in f. f, PERMANOVA shows that inter-individual variation is largest for all measurement types, with even relatively large effects (for example, antibiotics or IBD phenotype) capturing less variation (see Methods). Stratified tests (CD/UC) consider only samples within the indicated phenotype (note that sample counts decrease for these, resulting in larger expected covariation by chance). Stars show FDR-corrected statistical significance (FDR *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001). Variance is estimated for each feature independently (Methods). ‘All’ refers to a model with all metadata. Total n for each measurement type is shown in square brackets, distributed across up to 132 subjects (Extended Data Fig. 1a, see Methods).

Fig. 2. Metagenomic, metatranscriptomic, and stool metabolomic profiles are disrupted during IBD activity.

a, Relative abundance distributions for ten of the most cross-sectionally significantly differentially abundant metabolites in samples from individuals with IBD, as a ratio to the median relative abundance in individuals without IBD (Wald test; all FDR P < 0.003; see Methods; Supplementary Tables 1–14). Left, fraction of samples below detection limit (see Methods). n = 546 samples from 106 subjects. b, Relationships between two measures of disease activity: patient-reported (Harvey–Bradshaw index (HBI) in CD, n = 680 samples from 65 subjects; simple clinical colitis activity index (SCCAI) in UC, n = 429 samples from 38 subjects) and host molecular (faecal calprotectin (cal), n = 652 samples from 98 subjects). Linear regression shown with 95% confidence bound. c, d, Distribution of microbial dysbiosis scores as a measure of disease activity (c, median Bray–Curtis dissimilarity between a sample and non-IBD samples; see Methods) and its relationship with calprotectin (d, n = 652 samples from 98 subjects). Linear regression with 95% confidence. e, Kaplan–Meier curves for the distributions of the durations of (left) and intervals between (right) dysbiotic episodes in UC and CD. Both are approximately exponential (fits in dashed lines), with means of 4.1 and 17.2 weeks, respectively, for UC, and 7.8 and 12.8 weeks for CD (see Methods). f, Relative abundance distributions of significantly different metagenomic species (n = 1,595 samples from 130 subjects), metabolites (n = 546 samples from 106 subjects), and microbial transcribers (n = 818 samples from 106 subjects) in dysbiotic samples compared to non-dysbiotic samples from the same disease group (Wald test; all FDR P < 0.05; full results in Supplementary Tables 15–28). Also shown are antibody titres for ANCA, ASCA (IgG or IgA), anti-OmpC, and anti-CBir1 antibodies (n = 146 samples from 61 subjects). Boxplots show median and lower/upper quartiles; whiskers show inner fences; sample sizes above boxes.

Fig. 3. Temporal shifts in the microbiome are more frequent and more extreme in IBD.

a, Bray–Curtis dissimilarities within subjects as a function of intervening time difference, as compared to different people or technical replicates; calculated for metagenomic taxonomic profiles (species; n = 1,595 samples from 130 subjects), metabolomics (n = 546 samples from 106 subjects), and functional profiles (KO gene families; n = 818 samples from 106 subjects). Boxplots show median and lower/upper quartiles; whiskers show inner fences. Blue, least-squares power-law fits; orange, thresholds for microbiome shifts (see Methods). Proteomics and species-level transcripts in Extended Data Fig. 5a. Within-subject changes are significantly more extreme in UC and CD than in non-IBD for taxonomic profiles (F-test P = 3.9 × 10−10 and 1.2 × 10−18, respectively) and transcripts (P = 0.00016 and 1.7 × 10−5), with mixed differences for metabolites (P = 0.012 and 0.23). Technical replicates shown (when possible) at 0 weeks. b, Shift frequencies for the top 10 species with greatest change during shifts, ranked by number of shifts as primary contributor, stratified by disease phenotype(s) (full table Supplementary Table 29). c, P. copri is of interest in arthritis and international populations, and it alone retained stable abundances in CD but bloom-relaxation dynamics in controls (two-tailed Wilcoxon test of absolute differences between consecutive time points P = 4.2 × 10−6 between non-IBD and UC, and 1.1 × 10−4 between non-IBD and CD). Plot shows 22 subjects with at least one time point with more than 10% differential abundance (n = 267 samples). d, Ordination of temporally adjacent samples within individual, based on metabolomics (Bray–Curtis principal coordinates on normalized absolute abundance differences). Disease groups separate significantly (n = 440 sample pairs from 106 subjects; PERMANOVA R2 = 2.8%, P < 10−4). Urobilin, urate, and an unidentified untargeted feature that segregates with disease groups in the PCoA are shown (right); HILn_QI1594 (HILIC-neg method m/z = 152.0354, RT = 4.16 min). e, As in c, but for urate (two-tailed Wilcoxon test P = 0.0012 non-IBD–UC, P = 0.044 non-IBD–CD; n = 546 samples from 106 subjects).

Fig. 4. Colonic epithelial molecular processes perturbed during IBD and in tandem with multi-omic host–microbe interactions.

a, Human DEGs (negative binomial FDR P < 0.05, minimum fold change 1.5; Supplementary Table 31) from 81 subjects with paired ileal and rectal biopsies. Ordering by diagnosis, clustering within diagnosis. IL-17 signalling (I) showed strongest enrichment in ileal DEGs (FDR P = 8.2 × 10−12), while the complement cascade (II) was enriched in rectal DEGs from patients with UC (FDR P = 5.2 × 10−8; KEGG gene sets, Supplementary Table 32). Example DEGs shown with I and II. b, Expression of four genes involved in host–microbe interactions–. Inflamed biopsy samples are shown for CD from ileum (left, n = 20, 23, 39 independent samples for non-IBD, UC, CD respectively); for CD and UC in rectum (right column; n = 22, 25, 41 independent samples for non-IBD, UC, CD); non-IBD samples were non-inflamed. Asterisks indicate significant differential expression compared to non-IBD (Fisher’s exact test, FDR P < 0.05; P values in Supplementary Table 31). Boxplots show median and lower/upper quartiles; whiskers show inner fences. c, Significant associations among 10 aspects of host–microbiome interactions: metagenomic species, species-level transcription ratios, functional profiles captured as EC gene families (MGX, MTX and MPX), metabolites, host transcription (rectum and ileum), serology, and calprotectin (sample counts in Fig. 1b, c). Network shows top 300 significant correlations (FDR P < 0.05) between each pair of measurement types (for serology, FDR P < 0.25). Nodes coloured by disease group in which they are ‘high’, edges by sign and strength of association. Spearman correlations use residuals of a mixed-effects model with subjects as random effects (or a simple linear model when only baseline samples were used (biopsies)) after covariate adjustment (see Methods). Time points approximately matched with maximum separation 4 weeks (see Methods). Singletons pruned for visualization (Extended Data Fig. 8). Hubs (nodes with at least 20 connections) emphasized.

Extended Data Fig. 1. Distribution of sample and measurement types and timing.

a, Measurements available over time for each IBDMDB participant. b, Number of identical or closely-aligned time points available for which each measurement type has been generated (see Methods). c, Distributions of the number of processed samples per subject, stratified by disease and measurement type. d, Distributions of time intervals between consecutive samples for each measurement type.

Extended Data Fig. 2. Within-individual stability is a major driver of microbiome differences across measurement types.

a, PCoA and t-SNE embeddings based on Bray–Curtis dissimilarity matrices from stool species abundances, transcripts, proteins, and metabolites. Marginal densities are shown for the PCoAs that show disease separation for some measurements. In the t-SNEs, each subject has been assigned a different hue, showing that small clusters generally represent individuals’ time courses, as inter-individual differences are the greatest driver of microbiome variation across measurement types (Fig. 1f). Sample counts are shown in Fig. 1b. b, Distributions of correlations between functional profiles, captured as UniRef90 gene family abundances, measured from paired metagenomes, metatranscriptomes, and metaproteomes (see Methods). c, Human transcriptional expression was mostly determined by biopsy location rather than IBD phenotype and inflammation (n = 249 samples from 91 subjects). Ordination shows PCA on gene expression levels normalized by library sizes and represented as CPM. Ellipses indicate 95% confidence regions for the indicated sample types. d, Principal coordinates plot (Bray–Curtis on OTU profiles) of community profiles from biopsy samples shows that mucosal microbial communities do not differ significantly by biopsy location (shape), unlike human gene expression in epithelial tissue (Fig. 4b).

Extended Data Fig. 3. Patient-reported, molecular, and microbial disease activity measures.

a, Relationships between six measures of disease activity: the patient-reported HBI in CD, patient-reported SCCAI in UC, faecal calprotectin, the fractions of human reads from stool MGX and MTX, and a dysbiosis score defined here as departures from control population microbiome configurations (Fig. 2c). Rho values are Spearman correlations with ties broken randomly. Linear regression is shown (red line) with 95% confidence bound (shaded). Sample counts are presented in the title bar, though sample counts for a particular correlation may be less as samples must be paired. b, c, PCoA based on metagenomic species-level Bray–Curtis dissimilarities (n = 1,595 samples from 130 subjects), indicating dysbiosis score (b) and whether the sample was defined as dysbiotic (c). d, Number of dysbiotic samples per participant. Colour scheme as in c. e, Relationship between the dysbiosis score, when using the HMP1-II gut data set as reference (n = 553 from 249 subjects), compared to the non-IBD data set. The threshold for the dysbiosis classification is also shown (black lines). The two scores are highly correlated (Pearson ρ = 0.86; two-sided P < 2.2 × 10−16), as are the resulting dysbiosis classifications (odds ratio of 56; Fisher’s exact test P < 2.2 × 10−16). n = 1,595 samples from 130 subjects.

Extended Data Fig. 4. Significant microbial and metabolic perturbations during taxonomic dysbioses.

a, Alpha diversity (Gini–Simpson) as a function of the dysbiosis score for the sample (Pearson correlation –0.60; P < 2.2 × 10−16). n = 1,595 samples from 130 subjects. b, Seven (chosen for space constraints) most differentially abundant species not shown elsewhere in this manuscript (n = 1,595 samples from 130 subjects; Wald test; see Methods; full results in Supplementary Table 15). c, Top 10 differentially abundant acylcarnitines in dysbiosis (n = 546 samples from 106 subjects; Wald test). Carnitine and acylcarnitines are more abundant in dysbiotic CD, whereas C20:4 carnitine is significantly depleted (Supplementary Table 16). d, Top 10 differentially abundant metabolites during dysbiosis not shown elsewhere in this manuscript (n = 546 samples from 106 subjects; Wald test; full results in Supplementary Table 16).

Extended Data Fig. 5. Detecting shifts in longitudinal microbiome multi-omics.

a, Distributions of Bray–Curtis dissimilarities as a function of time difference between samples for protein profiles, species-level transcriptional activity (see Methods), and species-level taxonomy (though excluding subjects with dysbiosis at any time point), otherwise as in Fig. 3a. Removing subjects with dysbiotic samples removes the extreme dissimilarities (near 1) observed in IBD subjects. Boxplots show median and lower/upper quartiles; whiskers show inner fences. b, Distribution of Bray–Curtis dissimilarities between samples from the same subject, two weeks apart versus those from different individuals, allowing us to define a ‘shift’ in the microbiome as a change more likely to have been drawn from the between-subject distribution than within-subject distances (corresponding to Bray–Curtis > 0.54). c, Relative abundance differences of the top ten microorganisms that contributed to each of the 183 detected taxonomic shifts among any two within-subject subsequent time points. Shifts are typically reciprocal (that is, losing a microorganism and regaining it later, or vice versa), and microorganism with frequent high-abundance shifts generally correspond to frequent contributors in Fig. 3b. Sample ordering is from a hierarchical clustering using average linkage followed by optimal leaf ordering. d, As in Fig. 3c, but for E. coli (n = 322 samples from 24 subjects; two-tailed Wilcoxon test of the absolute differences in relative abundances between consecutive time points P = 2.2 × 10−4 for non-IBD to UC, and P = 0.029 for non-IBD to CD), which is frequently implicated in gut inflammation. e, As in b, but showing Bray–Curtis dissimilarities of metabolomic profiles. Here, 22% (96 out of 440) of sample pairs exceed the shift threshold, whereas 13% (183 out of 1,413) exceed the threshold in b. If metagenomic profiles are sub-sampled to match the metabolomics samples, this increases to 14% (57 out of 398) of sample pairs, showing that if we increased the sampling rate, this measurement type would be likely to shift more than the metagenomes. f, As in Fig. 3b, but showing the primary contributors to metabolomic shifts, that is, the metabolite with the largest change in relative abundance during a shift. Note that other metabolites may still experience large changes in abundance (for example, for this reason, urate was not a primary contributor to any non-IBD shifts, though large changes are visible for one non-IBD individual in Fig. 3e). The full table of detected metabolomic shifts is given in Supplementary Table 30. Violin plot shows the density of points around that intake frequency; bandwidth chosen automatically by Silverman’s method.

Extended Data Fig. 6. Mucosal communities and human genetics.

a, Number of biopsy samples available for each biopsy location and inflammation status. b, DEGs (Fig. 4a) with newly identified significant correlations with OTU abundances in biopsies (partial Spearman correlation conditioned on disease status, BMI, age at consent and sex; FDR P < 0.05; n = 54 in ileum and n = 52 independent 16S–RNA-seq pairs; full table in Supplementary Table 33). c, A limited subset of the microbiome trended with genetic variants in targeted testing, including the strongest trend shown here of Parabacteroides distasonis with genotypes of NKX2-3 (a known IBD-associated locus; boxplots show median and lower/upper quartiles; whiskers show inner fences). This is the most significant association by P value among all tested associations between metagenomic taxa and five known IBD loci (nominal significance P = 0.006; no associations passed FDR P < 0.05, mixed effect model with age, sex, antibiotic and immunosuppressant use and first 20 genetic principal components as covariates while specifying subjects as random effects; Wald test; n = 84 subjects of European ancestry with exomes and 960 metagenomes; full results in Supplementary Table 34). d, Association between rs1042712 SNP in the LCT locus and self-reported milk intake from dietary recall. Self-reported short-term milk intake (from dietary recalls accompanying stool samples) was significantly associated with the count of C alleles (29.8% allele frequency) at rs1042712 in the LCT gene locus using a linear mixed effect model accounting for age, sex, first 20 genetic principal components and with subjects as random effects (P = 0.028, linear mixed effect regression with Wald test, see Methods). All available data are plotted for unique subjects of European ancestry with exome data (per-genotype subject count (GG/GC/CC):50/26/8). Differences between IBD and non-IBD groups are not statistically significant (odds ratio 0.27; 95% CI 0.05–1.33; P = 0.10; n = 84 subjects of European ancestry with exomes and 960 dietary surveys; model: IBD (yes|no) ~ intercept + SNP + sex + age + PC1–PC20).

Extended Data Fig. 7. Microbial and host-related subsets of the multi-omic association networks.

a, Subset of the network in Fig. 4c showing metagenomic abundances (octagons) and expression levels (hexagons) of Subdoligranulum, Roseburia spp. and F. prausnitzii and their neighbours (three functionally associated microbial hubs selected for further investigation based on anti-inflammatory associations in the literature, see text). b, The host expression-related subnetwork of the ‘unadjusted’ association network (Extended Data Fig. 9, Supplementary Discussion). Sample counts in Fig. 1b, c.

Extended Data Fig. 8. Significant covariation among multi-omic components of the gut microbiome and host interactors in IBD (adjusted).

Detailed labelling of the association network in Fig. 4c (intended for magnification). The network was constructed from ten data sets: metagenomic species, species-level transcription ratios, functional profiles at the EC levels (MGX, MTX and MPX), metabolites, host transcription (rectal and ileal separately), serology and faecal calprotectin. As in Fig. 4c, measurement types were approximately matched in time with a maximum separation between paired samples of four weeks. The top 300 significant correlations (FDR P < 0.05) among correlations between features that were differentially abundant in dysbiosis were used to construct the network visualized here (for serology, a threshold of FDR P < 0.25 was used). Nodes are coloured by the disease group in which they are ‘high’, and edges are coloured by the sign and strength of the correlation. For this adjusted network, Spearman correlations were calculated using HAllA from the residuals of a mixed-effects model with subjects as random effects (or a simple linear model without the random effects when only baseline samples were used) after adjusting for age, sex, diagnosis, dysbiosis status, recruitment site, and antibiotics (see Methods). Appropriate normalization and/or transformation for each measurement type was performed independently before the model fitting (see Methods). Singleton node pairs were pruned from the network. Source associations are in Supplementary Table 35, sample counts in Fig. 1b, c.

Extended Data Fig. 9. Significant covariation among multi-omic components of the gut microbiome and host interactors in IBD (unadjusted).

The network was constructed from ten data sets: metagenomic species, species-level transcription ratios, functional profiles at the EC levels (MGX, MTX and MPX), metabolites, host transcription (rectal and ileal separately), serology and faecal calprotectin. As in Fig. 4c, measurement types were approximately matched in time with a maximum separation between paired samples of four weeks. The top 300 significant correlations (FDR P < 0.05) among correlations between features that were differentially abundant in dysbiosis were used to construct the network visualized here (for serology, a threshold of FDR P < 0.25 was used). Nodes are coloured by the disease group in which they are ‘high’, and edges are coloured by the sign and strength of the correlation. For this unadjusted network, Spearman correlations were calculated using HAllA from the residuals of the same model as in Extended Data Fig. 8, though without adjusting for dysbiosis (see Methods). Appropriate normalization and/or transformation for each measurement type was performed independently before the model fitting (see Methods). Singleton node pairs were pruned from the network. Source associations are in Supplementary Table 36, sample counts in Fig. 1b, c.