First insights into microbial changes within an Inflammatory Bowel Disease Family Cohort study

Abstract

The prospective Kiel Inflammatory Bowel Disease (IBD) Family Cohort Study (KINDRED cohort) was initiated in 2013 to systematically and extensively collect data and biosamples from index IBD patients and their relatives, a population at high risk for IBD development. Regular follow-ups were conducted to collect updated health and lifestyle information, to obtain new biosamples, and to capture the incidence of IBD during development. By combining microbial data collected at successive time points with extensive anthropometric, medical, nutritional, and social information, this study aimed to characterize the factors influencing the microbiota in health and disease via detailed ecological analyses. Using a microbial dysbiosis metric based on the German KINDRED cohort, we identified strong and generalizable gradients within and across different external IBD cohorts for validation. These community gradients correspond strongly with IBD pathologies, physiological manifestations of inflammation (e.g. Bristol stool score, ASCA IgA, ASCA IgG), and genetic risk for IBD. Anthropometric and medical factors influencing fecal transit time strongly modify bacterial communities. Various Enterobacteriaceae (e.g. Klebsiella sp.) and opportunistic Clostridia pathogens (e.g. C. XIVa clostridioforme), characterize in combination with ectopically colonizing oral taxa (e.g. Veillonella sp. Cand. Saccharibacteria sp. Fusobacterium nucleatum) the distinct and chaotic IBD-specific communities. Weak community and physiological changes are further traceable in a small number of individuals, who developed IBD in the study's runtime. Our findings demonstrate broad-scale ecological patterns which indicate drastic state transitions of communities in IBD patients. These patterns appear to be universal across cohorts and influence physiological signs of inflammation, display increased resilience, but show only limited heritability/intrafamily transmission.

Keywords: 16S; Family Cohort; IBD; dysbiosis; inflammatory bowel disease; microbiota; oralization.

Figure 1.

(A) the flow diagram displays the change of the cohort size and its subsets between time points. The highlighted transitions are onset cases detected during the study’s runtime (N_BL→F1=4 (CD:2, UC:2), N_F1→F2=1 (CD:1), N_F2→F3=2 (uIBD:2)). (B) Analyses of selected physiological inflammation markers (ASCA IgA, ASCA IgG, calprotectin, Bristol stool score) with respect to IBD condition and relevant covariates via linear models at baseline (see Table 1, Table S5 for additional time points). Optimal model results (min. AIC) are shown and average pairwise differences with respect disease groups and/or differences in group slopes/direction of association are highlighted in the individual plots (# p≤0.1000, * p≤0.0500, ** p≤0.0100, *** p≤0.0010). (C) The violin plot displays the average differences of genetic predisposition to CD, UC, and IBD in general, as based on LDpred2 derived polygenic risk scores (PRS). The average differences in PRS were tested via Wilcoxon rank tests contrasting healthy individuals (N=785), with healthy future onset cases (grey; N_BL=4, N_F1=1, N_F2=2), CD patients (orange; N=551), UC patients (red; N=438), and patients with unclassified IBD (uIBD; darkred; N=32). In particular, compared with healthy controls, the average risk of IBD is significantly and consistently higher in patients with CD and UC, as well as in future onset cases. (D) Scatterplots show the significant relationships between CD-PRS and selected physiological inflammation markers at the baseline time point, considering anthropometric covariates (linear model on residuals) and general correlation without additional covariates (Spearman correlation). Additional analyses of IBD-PRS and UC-PRS with physiological markers can be found in Figure S1 and Table S6.

Figure 2.

Phylum and ASV level abundances with respect to IBD status, association to clinical measures of inflammation and analyses of alpha diversity. (A) Overview of individual phylum abundances across time points and health conditions in the KINDRED cohort. (B) Differential abundance analyses of ASVs based on the baseline time point (BL). Displayed are the log Fold changes for each significantly different ASV clustered by genus classification, including standard errors. Color coding indicates the phylum membership of the ASV. The plot only displays taxa with significant differential abundance in the respective comparisons. (C) Partial correlation of CLR transformed taxon abundances with core physiological measures via ppcor, combining the p values of Spearman-, Kendall-, and Pearson correlations via Brown’s method and corrected via FDR. Correlations were adjusted for age, gender, and BMI and Spearman ρ is used to visualize correlation strength between taxa and clinical measures. Overlapping patterns of differential abundance for the respective taxa in the KINDRED cohort, Maltese-, and Swedish cohort is indicated in the bottom color bars. (D) Differences in Chao1 species richness between healthy individuals, future onset cases, CD, UC, and uIBD cases at baseline with additional onset cases (N_onset = 7, N_CD = 551, N_UC = 438, N_uIBD = 32, N_contr. = 785; Wilcoxon test; see Materials and Methods). (E) Correlation of species richness with the microbial dysbiosis index (MD-index), which results in a negative, but non-linear relationship (optimal AIC based fit, Table S8) between diversity and dysbiosis (BL (polynomial): F_2,1809=226.97, p<2.2×10⁻¹⁶, adj.R² = 0.1997, N = 1812; F1 (polynomial): F_2,645 = 83.778, p<2.2×10⁻¹⁶, adj.R² = 0.2037, N = 648; F2 (polynomial): F_2,536 = 55.745, p<2.2×10⁻¹⁶, adj.R² = 0.1691, N = 539; linear model). (F) Relationship between Chao1 species richness measured at baseline and several IBD relevant clinical measures, which shows significant relationships between alpha diversity and host physiology, in a disease-specific manner (see Table S9, excluding uIBD). Average pairwise differences between disease groups, or differences in slope/direction of association are highlighted in the plots (# p≤0.1000, * p ≤ 0.0500, ** p ≤ 0.0100, *** p ≤ 0.0010).

Figure 3.

(A) Non-metric multidimensional scaling (NMDS) of Bray-Curtis distances among baseline samples, displaying the significant clustering by health conditions and significant correlations of clinical inflammation measures with community distance (BL, see Table S16 & Table S17). (B) NMDS displaying the gradient of community dysbiosis as expressed by MD-index, in parallel with clinical measures of inflammation and healthy onset cases highlighted in red (*, develop IBD in F1). (C) Correlation of MD-index and the first NMDS axis showing a clear gradient of dysbiosis in the community. Onset cases are distributed within the range of standard deviation around the mean of the community distribution (NMDS1) and the MD-index. (D) Community variability between health conditions as measured by the distance to the group centroid, is overall significantly different between health conditions (F_3,1808 = 46.0315, p = 0.00001, PERMANOVA), and significantly increased in CD and UC patients as compared to healthy controls (Figure S11; Table S18). (E) Principle coordinate analysis of German-, Swedish-, and Maltese samples, highlighting the transferability of the dysbiosis gradient across cohorts (MD-index derived from German samples), (F) as well as a common disease wise clustering of communities irrespective of sample origin (Table S20). (G) Community variability between health/IBD conditions within and between the German-, Swedish-, and Maltese cohorts showing an increased variability in IBD cases. (H) Mean differences of dysbiosis (MD-index derived from German samples) within and across cohorts, with the strongest differences between healthy and CD individuals.(I) Line plots visualize the correlation of selected physiological inflammation markers with the microbial dysbiosis index (MD-index), and show disease-specific differences as compared to healthy control individuals. Significant pairwise comparisons with respect to average differences between pathologies, or differences in slope are highlighted in the plots (see Figure S12, Table S20). (J) Visualization of the explained variation of pathologies significantly associated to community composition based on serial PERMANOVA of Bray-Curtis distances in all three time points (Table S21). Variables are displayed if they show significant clustering in at least one time point (# P_FDR≤0.1000, * P_FDR ≤ 0.0500, ** P_FDR≤0.0100, *** P_FDR ≤ 0.0010). for beta diversity analyses including physiological variables, pharmaceutical treatments, and nutrient intake see Table S16, Table S21–23 and Figure S9 & S10.

Figure 4.

(A) Community clusters of the microbial community at the baseline time point, determined by dirichlet multinomial mixture modelling (DMM) and optimal clustering was determined via Laplace goodness of fit optimization. (B) Overlay of the microbial dysbiosis index gradient (Gevers et al.) and community clusters (outlines), including healthy onset patients at baseline (indicated by *). (C) Community clusters display a significantly elevated level of dysbiosis in clusters 2 and 3 (pairwise Wilcoxon tests), (D) as well as elevated levels of community variability in cluster-3 (PERMANOVA). (E) Violinplots visualize the differences of inflammation related physiological variables between community clusters, highlighting elevated levels of inflammatory biomarkers in cluster-3 (pairwise Wilcoxon test).

Figure 5.

(A) Analyses of selected physiological fecal inflammation markers (fecal calprotectin, TNFα, IL6, IL1β, Zonulin) with respect to pre- and post-diagnosis in onset patients (ANOVA-P value, paired t-test-P_paired value) show no significant differences between disease states. Individuals with paired samples (pre- and post- onset/diagnosis) are connected by lines (N=4). (B) the violin plots display nominally significant differentially abundant genera, between pre- and post-diagnosis onset patients (ANOVA-P value, paired t-test-P_paired value), as well as (C) increasing Shannon diversity in post-diagnosis onset patients, insignificantly elevated Chao1 levels and (D) dysbiosis (MD-index). (E) Principle coordinate analysis of Bray-Curtis distances between onset patients pre- and post-diagnosis (N_Contr.=7, N_CD=3, N_UC=1). Individuals with paired samples (pre- and post- onset/diagnosis) are connected by lines (N=4). (F) non-metric multidimensional scaling (NMDS) of Bray-Curtis distances between all samples across all time points (BL, F1, F2), displaying the significant clustering by health conditions (F_2,2951=6.77216, p<0.00010, R²=0.00469, adj. R²=0.00389, PERMANOVA, excluding uIBD) and (G) the gradient of community dysbiosis as expressed by MD-index. The NMDS further highlights onset cases before and after disease manifestations and shows directly paired samples as connected. The marginal boxplots display differences between NMDS scores of these paired onset samples (paired t-test). (H) Correlation of the MD-index and the first NMDS axis showing a clear gradient of dysbiosis in the community composition (ρ=-0.7569, p<2.20 × 10⁻¹⁶; Spearman rank correlation) and the distribution of onset cases around the mean of the community distribution (NMDS1). P-values were not adjusted for multiple testing.

Figure 6.

(A) Spiec-Easi networks of baseline samples (N_BL=1812). Bacterial nodes highlight significant differentially abundant ASVs in the network. Bacteria not showing any differential abundance patterns between IBD patients and healthy controls are signified via (●), bacteria overabundant in IBD (combined CD, UC, uIBD) via (■) and bacteria more abundant in controls are signified via (►). The barplot visualizes node centrality based on the number of connections (node degree) at the baseline time point (see Figure S14-S16 and Table S24 for BL, F1, and F2). Colored boxes highlight corresponding ASV differential abundance patterns in KINDRED and the external cohorts. Significance of centralities is derived from Z-tests against a randomized networks (10’000) and nodes with significantly higher degrees than expected by chance are. (B) Global network characteristics were derived from networks constructed from the disease condition-specific networks of each sampling timepoint, which are informative for stability and structure of the respective networks (Wilcoxon-test healthy vs. CD/UC). Network metrics include centrality based assortativity, network diameter, radius, and size as well as density/clustering,^, and natural connectivity. (C) Network similarity of disease and time point-specific subnetworks, as well as networks derived from the external cohorts (Malta, Sweden) based on graphlet distance and displayed via NMDS (see Materials and Methods section). Networks display compositional differences between healthy and diseased networks (Control vs. IBD (incl. IBD networks): F_1,18=2.4144, p=0.0412, R²=0.1183, adj. R²=0.0693; PERMANOVA). (D) Heritability estimates derived from the likelihood based method lme4qtl using either only kinship information with or without additional environmental and anthropometric covariates (▲-incl. environmental covariates (h²_full), ▼-no covariates (h²_null)). Only the upper quartile of taxa are highlighted (based on h²_full estimate including environmental covariates). Additional information like differential abundance in IBD accross cohorts are depicted for each taxon (Figure S19, Table S25).