Bacteria and bacteriophage consortia are associated with protective intestinal metabolites in patients receiving stem cell transplantation

Abstract

The microbiome is a predictor of clinical outcome in patients receiving allogeneic hematopoietic stem cell transplantation (allo-SCT). Microbiota-derived metabolites can modulate these outcomes. How bacteria, fungi and viruses contribute to the production of intestinal metabolites is still unclear. We combined amplicon sequencing, viral metagenomics and targeted metabolomics from stool samples of patients receiving allo-SCT (n = 78) and uncovered a microbiome signature of Lachnospiraceae and Oscillospiraceae and their associated bacteriophages, correlating with the production of immunomodulatory metabolites (IMMs). Moreover, we established the IMM risk index (IMM-RI), which was associated with improved survival and reduced relapse. A high abundance of short-chain fatty acid-biosynthesis pathways, specifically butyric acid via butyryl-coenzyme A (CoA):acetate CoA-transferase (BCoAT, which catalyzes EC 2.8.3.8) was detected in IMM-RI low-risk patients, and virome genome assembly identified two bacteriophages encoding BCoAT as an auxiliary metabolic gene. In conclusion, our study identifies a microbiome signature associated with protective IMMs and provides a rationale for considering metabolite-producing consortia and metabolite formulations as microbiome-based therapies.

Extended Data Fig. 1

Longitudinal dynamics of bacterial, fungal and viral community compositions in patients receiving allo-SCT. a) Quantification of the bacterial and fungal load by 16S and 28S rDNA copy numbers per gram of stool at time-points relative to all-SCT (Day 0). Solid lines represent the smoothed conditional means for the entire cohort (black), Munich (MUC, orange) or Regensburg (REG, red) calculated by locally weighted regression using the locally estimated scatterplot smoothing (LOESS) method. The gray shading indicates the 95% confidence interval for entire cohort. Each individual patient stool sample is plotted as a gray dot superimposed on the graph. The number of samples is indicated. b) Beta diversity analysis illustrating changes in bacteriome, fungome and virome by time-points relative to allo-SCT (Day 0). For bacteriome and fungome, beta diversity was calculated by weighted UniFrac. For virome, beta diversity was calculated by Bray-Curtis dissimilarity. Distances were projected in Principal coordinate analysis (PCoA). Bacteriome: Comparisons between Day +7 (p=0.001), Day + 14 (p=0.001), Day +21 (p=0.001) , Day +28 (p=0.001) vs baseline (Day -7) are significant (pairwise Adonis test adjusted for multiple comparisons). Virome: Comparisons between Day +7 (p=0.013) and Day +28 (p=0.001) vs baseline are significant (PERMANOVA test). The number of samples is indicated. c) Beta diversity analysis illustrating changes in bacteriome, fungome and virome according to study center (MUC or REG) and whether patients received antifungal therapy (”No Antifungals” or “Antifungals”). For bacteriome and virome, beta diversity was calculated by weighted UniFrac. For virome, beta diversity was calculated by Bray-Curtis dissimilarity. Distances were projected in Principal coordinate analysis (PCoA). The number of samples is indicated.

Extended Data Fig. 2

Longitudinal dynamics of SCFAs, BCFAs, IIMs, primary and secondary bile acids in patients receiving in allo-SCT a) Heatmap of normalized Panel 1 metabolite levels in stool samples of allo-SCT patients averaged by time-points relative to allo-SCT. Normalized concentrations are indicated in the adjacent color legend. Clustering based on metabolite expression patterns using the Ward algorithm. Distance measure is Euclidian. The number of samples is indicated. b) Heatmap of normalized Panel 2 metabolite levels in stool samples of allo-SCT patients averaged by time-points relative to allo-SCT. Normalized concentrations are indicated in the adjacent color legend. Clustering based on metabolite expression patterns using the Ward algorithm. Distance measure is Euclidian. The number of samples is indicated. c) Principal component analysis (PCA) of Panel 1 metabolite profiles by time-points relative to allo-SCT. Comparisons at Day 0 (p=0.004), Day +7 (p=0.001), Day +14 (p=0.001), Day +21 (p=0.001) and Day +28 (p=0.001) vs Day -7 are significant (pairwise Adonis test of Euclidean distances adjusted for multiple comparisons). d) PCA of Panel 2 metabolite profiles by time-points relative to allo-SCT. Comparisons at Day +7 (p=0.001), Day +14 (p=0.001), Day +21 (p=0.001) and Day +28 (p=0.001) vs Day -7 are significant (pairwise Adonis test of Euclidean distances adjusted for multiple comparisons).

Extended Data Fig. 3

MOFA (multi-omics factor analysis) and MEFISTO (a method for the functional integration of spatial omics data) in allo-SCT patients a) Correlation between MOFA-identified Factors and normalized intestinal metabolite concentrations. Associations between Factor values and metabolites were analyzed via Pearson‘s correlation. The correlation coefficient is indicated in the adjacent color legend. The p-values associated with the correlations were corrected for multiple testing with the FDR approach. b) Top 15 Features in bacteriome, virome and metabolites in Factor 4 in descending order according to Feature weight. Larger weights indicate a higher correlation with that Factor, while the positive or negative sign indicates the directionality of that variation, i.e., “+” indicates a positive association, “-“ a negative association. c) Bar plot of time scale parameters assigned to Factors 1 through 10 identified by MEFISTO. MEFISTO assigns a time scale value between 0 and 1 to each Factor, which reflects the degree to which that Factor is dependent on time. A value of 0 implies no time-dependency, a value of 1 strong time-dependency. Of note, results pertaining to the identified Factors, their weight/covariance structure, variance explained across omics entities and Factor values obtained by MEFISTO modelling were almost identical and thus comparable to the output of our original MOFA model. d) Heatmap of normalized abundance of viral contigs assigned to eukaryotic and prokaryotic viruses at time-points relative to allo-SCT. The number of samples is indicated.

Extended Data Fig. 4

Correlation between top 15 bacterial and metabolite as well as bacterial and viral high-weight Features in Factors 1, 3 and 4. a) Heatmaps of pairwise Pearson‘s correlations of top Features across different omics modalities. The Feature values of the top 15 high-weight Features of a given omics modality were correlated with that of another omics modality. (Left) bacterial taxa at genus level and metabolites, (right) bacterial taxa at genus level and bacteriophages at species level. The correlation coefficient is indicated in the adjacent color legend. The p-values associated with Pearson correlation have been corrected for multiple testing by applying the FDR approach to each set of correlations of two omics modalities of a given Factor. b) As in a) for Factor 3. c) As in a) for Factor 4.

Extended Data Fig. 5

Co-abundance of MOFA-identified bacterial and viral Features is associated with high-level IMM expression, which declines progressively after allo-SCT. a) Correlation scatter plots of high-weight Features within Factor 4, comparing normalized abundance of bacterial taxa at genus level (x-axis) with that of bacteriophages at species level (y-axis) together with the metabolite propionic acid. Dots represent samples from individual patients at different time points (91 samples from 45 patients), all of which had both 16S and viral metagenomic sequencing data. Dots are colored by intestinal levels of propionic acid, or in grey if no propionic acid data was available. Normalized concentrations of metabolites are indicated in the adjacent heatmap. Associations between bacterial genera and viral species were analyzed via Pearson correlation and linear regression. The R- and p-values are indicated in each plot. The regression line is drawn in blue and the 95 % confidence interval of the regression line is shaded in grey. b) As in a) for Features in Factor 3 and the metabolite isovaleric acid. c) As in a) for Features in Factor 3 and the metabolite DAT. d) As in a) for Features in Factor 3 and the metabolite ICA. e) Levels of intestinal microbiota-derived metabolites at time-points relative to allo-SCT (Day 0) in μmol per gram of dried stool measured by targeted mass spectrometry. Number of patients per time-point is indicated in Figure 4D. Significance by two-sided Kruskal-Wallis test corrected for multiple testing via Dunn’s test of all time-points against baseline (Day -7). In the box plots, the box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the 10th and up to the 90th percentile. Points below and above the whiskers are drawn as individual points. Points indicate individual patient stool samples sampled at the specified time-points. SBA… secondary bile acids; UDCA…ursodeoxycholic acid; TUDCA…tauroursodeoxycholic acid.

Extended Data Fig. 6

Intestinal bacterial but not fungal nor viral diversity predicts outcome after allo-SCT a) 2-year overall survival after Day 21 stratified according to higher and lower fungal (left) and viral (right) alpha diversity. The mean alpha diversity of patient samples at Days +7–21 was calculated and patients were stratified into higher (blue curve) and lower (red curve) diversity groups, defined as above or below the center-specific median Inverse Simpson’s diversity index. For fungome, there were 12 deaths among 32 patients in the lower-diversity group (estimated mean survival time 521 (95% CI 428–614) days) and 10 deaths among 30 patients in the higher-diversity group (estimated mean survival time 597 (95% CI 522–672) days). For virome, there were 4 deaths among 16 patients in the lower-diversity group (estimated mean survival time 587 (95% CI 479–696) days) and 4 deaths among 15 patients in the higher-diversity group (estimated mean survival time 576 (95% CI 450–701) days). Analysis via Kaplan–Meier estimator, significance according to the log-rank test. b) 2-year cumulative incidence of relapse and transplantation-related mortality (TRM) in a competing risks analysis stratified according to higher and lower bacterial (top), fungal (middle) and viral (bottom) alpha diversity, calculated as in a). For bacteriome, there were 12 cases of TRM and 6 relapses among 35 patients in the lower-diversity group and one case of TRM and 8 relapses among 33 patients in the higher-diversity group. For fungome, there were 7 cases of TRM and 6 relapses among 32 patients in the lower-diversity group and 5 cases of TRM and 7 relapses among 30 patients in the lower-diversity. For virome, there was one case of TRM and 3 relapses among 16 patients in the lower-diversity group and 3 cases of TRM and 2 relapses among 15 patients in the higher-diversity. Significance according to Gray’s test. c) 2-year cumulative incidence of GvHD and its competing risk Death in a competing risks analysis stratified by alpha diversity as in a). For bacteriome, there were 12 cases of GvHD and 6 deaths among 35 patients in the lower-diversity group and 1 case of GvHD and 8 deaths among 33 patients in the higher-diversity group. For fungome, there were 7 cases of GvHD and 6 deaths among 32 patients in the lower-diversity group and 5 cases of GvHD and 7 deaths among 30 patients in the lower-diversity group. For virome, there was one case of GvHD and 3 deaths among 16 patients in the lower-diversity group and 3 cases of GvHD and 2 deaths among 15 patients in the higher-diversity group. Statistics as in b).

Extended Data Fig. 7

Characterization of differentially abundant microbial pathways and the species which encode them in IMM-RI low- vs high-risk patients via whole shotgun metagenomic sequencing a) Species-level association with MetaCyc pathways differentially abundant between IMM-RI low and high-risk patients (as shown in Figure 6). The relative abundance and species-level identity of taxa encoding the indicated MetaCyc pathways are shown. b) Box plot of relative abundance of indicated MetaCyc pathways in IMM-RI low- vs high-risk patients. In the box plots, the box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the mean. The whiskers are drawn down to the minimum and up to the maximum. Samples outside the 1.5-fold IQR were regarded as outliers. Patient samples are plotted as a point superimposed on the graph (high IMM-RI: n=10 patients; low IMM-RI: n=7 patients). Significance by 1-sided Wilcoxon Rank Sum Test corrected for multiple comparisons via the Benjamini & Hochberg correction for multiple testing. c) Box plot of relative abundance of acetic acid and propionic acid superclasses in IMM-RI low- vs high-risk patients. Plots, numbers and statistics as in b).

Extended Data Fig. 8

BCoAT-coding VC-1 and VC-2 bacteriophages are associated with MOFA Factors 1 and 3 and the IMM-RI a) Gene alignment plot of VC-2. The identity overlap (in percent) is indicated in the adjacent color legend. The BCoAT AMG is highlighted in red. b) Gene alignment plot of VC-1 as in a). c) Box plots of Factor values for Factors 1 and 3 (averaged across time-points Days +7–21) according to whether VC-1 was detected by viral metagenomic sequencing (“Yes”) or not („No“). The center line corresponds to the median, the box ranges from the 25^th to the 75^th percentile. Whisker length corresponds to the largest/lowest data point that does not exceed the 75^th/25^th percentile +/− 1.5-fold IQR. Blue: detected (n=16 patients); red: not detected (n=13 patients). Significance by two-tailed Mann-Whitney-U-test corrected for multiple testing via FDR. d) As in c) for VC-2. Blue: detected (n=13 patients); red: not detected (n=16 patients). e) Detection of the BCoAT-coding VC-1 in patient samples stratified according to IMM-RI. Bar plots show percentage, exact values are provided. The numbers of samples screened vs those in which VC-1 was detected is indicated below. f) As in e) for VC-2.

Extended Data Fig. 9

Onset of acute GI-GvHD shifts intestinal bacterial and viral communities and impacts IMM expression profiles a) Quantification of the bacterial and fungal load by 16S and 28S rDNA copy numbers per gram of stool stratified by patients diagnosed with GI-GvHD (GI-GvHD, red) vs control allo-SCT patients (No GI-GvHD, blue). The box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the minimum and maximum. Samples outside the 1.5-fold IQR were regarded as outliers. Significance by two-tailed Wilcoxon rank sum test adjusted for multiple comparisons via the Benjamini & Hochberg procedure. Each individual patient is plotted as a point superimposed on the graph. Number of patients: for 16S & 28S n=22 vs n=37, corresponding to „No GI-GvHD“ vs „GvHD“, respectively. b) Intestinal fungal alpha diversity (Richness, Inverse Simpson’s diversity index) stratified by patients diagnosed with GI-GvHD (GI-GvHD) vs control allo-SCT patients (No GI-GvHD). Plots, numbers and statistics as in a). c) Beta diversity analysis illustrating changes in bacteriome, fungome and virome in patients with GI-GvHD vs control allo-SCT patients. For bacteriome and fungome, beta diversity was calculated by weighted UniFrac. For virome, beta diversity was calculated by Bray-Curtis dissimilarity. Distances were projected in PCoA. Bacteriome: Comparison between patients with GI-GvHD vs No GI-GvHD is significant (p=0.026 by one-side pairwise Adonis test). Virome: Comparison between patients with GI-GvHD vs No GI-GvHD is significant (p=0.03 by PERMANOVA test). d) Levels of indicated microbiota-derived metabolites stratified by patients diagnosed with GI-GvHD vs control allo-SCT patients. Significance by two-tailed Mann-Whitney test. In the scatter-dot plots, the box is plotted at the mean. Error bars indicate standard deviation. Each individual patient is plotted as a point superimposed on the graph. The number of patients per group are indicated in the legend.

Extended Data Fig. 10

Impact of antibiotics on bacterial abundance, bacterial and viral community composition and IMM expression profiles a) Quantification of the bacterial and fungal load by 16S and 28S rDNA copy numbers per gram of stool stratified by antibiotic exposure: No Antibiotics (“No ABX”, blue) (blue) vs Antibiotics (“ABX”, red). Once a patient was treated with antibiotics, the current and all subsequent samples were classified as “ABX”. The box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the minimum and maximum. Samples outside the 1.5-fold IQR were regarded as outliers. Significance by two-tailed Wilcoxon rank sum test adjusted for multiple comparisons via the Benjamini & Hochberg procedure. Each individual patient is plotted as a point superimposed on the graph. Number of patients: for 16S n=59 vs n=70, for 28S n=56 vs n=70, corresponding to „No ABX“ vs „ABX“, respectively. b) Intestinal fungal alpha diversity (Richness, Inverse Simpson’s diversity index) in paired patient samples according to antibiotic status as in a). Plots, numbers and statistics as in a). c) Beta diversity analysis illustrating the impact of antibiotics on the intestinal bacterial, fungal and viral communities. Each point represents individual patient samples annotated with metadata regarding concomitant antibiotic therapy. For bacteriome and virome, beta diversity was calculated by weighted UniFrac. For virome, beta diversity was calculated by Bray-Curtis dissimilarity. Distances were projected in PCoA. Bacteriome: Comparison between patients with “No ABX” vs “ABX” is significant (p=0.001 by one-side pairwise Adonis test). Virome: Comparison between patients with “No ABX” vs “ABX” is significant (p=0.003 by PERMANOVA test). d) Levels of indicated microbiota-derived metabolites in paired patient samples before and after exposure to ABX. Significance by two-tailed Wilcoxon matched-pairs signed rank test. In the scatter-dot plots, the box is plotted at the mean. Error bars indicate standard deviation. Each individual patient is plotted as a point superimposed on the graph. The number of patients per group are indicated in the legend.

Figure 1:

Study design and sampling scheme a) Graphical representation of allo-SCT centers including patient numbers per center, sampling time-points including number of patients per time-point, sampling scheme including processing and storage, sample numbers per multi-omics modality, and downstream analytics. Created with BioRender.

Figure 2:

Longitudinal dynamics of intestinal bacteriome, fungome, virome and metabolome in patients receiving allo-SCT a) Intestinal microbiome composition in bacteriome (16S amplicon sequencing, n=266 samples), fungome (ITS1 amplicon sequencing, n=245 samples) and virome (viral metagenomic shotgun sequencing, n=138 samples) at time-points relative to allo-SCT (Day 0) by alpha diversity indices (Richness by number of observed zOTUs and Inverse Simpson’s diversity index, respectively). Solid lines represent the smoothed conditional means for the entire cohort (black), Munich (MUC, orange) or Regensburg (REG, red) calculated by locally weighted regression using the locally estimated scatterplot smoothing (LOESS) method. The gray shading indicates the 95% confidence interval of the respective alpha diversity index for the entire cohort. Each individual patient stool sample is plotted as a gray dot superimposed on the graph. The number of patients at each time-point is indicated in the legend. b) Heatmap of intestinal concentration of microbiota-derived metabolites (in μmol/g dry feces) normalized by log transformation (base 10) in stool samples (n=269 samples) of allo-SCT patients ordered by time-points relative to allo-SCT (Day 0). Panel 1 displays short-chain fatty acids (SCFA), branched-chain fatty acids (BCFA), type I-interferon inducing metabolites (IIMs), i.e., desaminotyrosine (DAT) and indole-3-carboxylaldehyde (ICA) as well as lactic acid. Panel 2 shows secondary bile acids (SBA). In each panel, then top 10 metabolites are shown, selected based on highest significance in ANOVA group comparisons between time-points. Clustering based on metabolite expression patterns using the Ward algorithm. Distance measure is Euclidian.

Figure 3:

Multi-omics factor analysis (MOFA) identifies bacterial and bacteriophage consortia associated with intestinal immuno-modulatory metabolites (IMMs) a) Bar plot of the total variance explained by MOFA (in percent) within each omics modality. Exact percentage values are superimposed on the bars. b) Heatmap of variance explained by individual MOFA-identified Factors, which capture major sources of variability within and across omics datasets. The rows show the percentage of variance explained by a given Factor. The columns indicate the omics modality. c) Bar plot showing the fraction of significant associations between the Features of each omics modality and the Factors 1, 3 and 4. Correlations between Features and Factors were analyzed via Pearson’s correlation coefficient. P values were corrected for multiple testing via FDR. Statistical significance is called at p≤0.05. Exact fraction values are superimposed on the bars. d) List of high-weight Factor Features within each omics modality in descending order according their corresponding Feature weight. For each Factor, MOFA learns a “weight” as a measure of importance of every Feature contained within that Factor hence enabling the interpretation of the variation captured by that Factor . A larger weight indicate a higher correlation with that Factor, while the positive or negative sign indicates the directionality of that variation, i.e., “+” indicates a positive association, “-“ a negative association. Shown are the top 15 Features in bacteriome, virome and metabolome for Factor 1. Viral taxa are listed by viral cluster name. Blue or pink font coloring indicates membership within the Oscillospiraceae and Lachnospiraceae families. e) As in d) for Factor 3. Orange coloring indicates bacteriophages that were correlated with Oscillospiraceae and Lachnospiraceae. f) Scatter plots of top Features within Factor 3, plotting the normalized abundance of bacteriophages at species level (y-axis) against that of bacterial taxa at genus level (x-axis). Dots represent samples from individual patients at different time points (91 samples from 45 patients). Dots are colored by intestinal levels of butyric acid, or in grey if no butyric acid data was available. Normalized concentrations of metabolites are indicated in the adjacent color legend. Associations between bacterial genera and viral species were analyzed via Pearson correlation and linear regression. The correlation coefficient R as well as the associated p-value is provided in each plot. The regression line is drawn in blue and the 95 % confidence interval of the regression line is shaded in grey.

Figure 4:

MOFA Factors are associated with outcome and levels of MOFA-identified IMMs decline progressively after allo-SCT a) Factor values of Factor 3 plotted against values of Factor 1 (averaged across Days +7–21). Points are colored and shaped according to the clinical outcome TRM (transplant-related mortality). Blue dots: no TRM at 2-year follow-up (n=58 patients), red crosses: deceased due to TRM at 2-year follow-up (n=15 patients). b) As in a), but Factor values of Factor 4 plotted against values of Factor 1. c) As in a), but Factor values of Factor 3 plotted against values of Factor 4. d) Levels of intestinal microbiota-derived metabolites at time-points relative to allo-SCT (Day 0) in μmol per gram of dried stool samples measured by targeted mass spectrometry. Metabolites identified as high-weight Features in MOFA Factors 1 and 3 are shown. Shown are (1) SCFA; (2) BCFA; (3) IIM and PBA. The number of patients at each time-point is indicated in the legend. Abbreviations are explained in the text box. Metabolites that comprise the Immuno-modulatory Risk Index (IMM-RI) are highlighted by background shading. Significance by two-sided Kruskal-Wallis test corrected for multiple testing via Dunn’s test of all time-points against baseline (Day −7). In the box plots, the box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the 10th and up to the 90th percentile. Points below and above the whiskers are drawn as individual points. Points indicate individual patient stool samples sampled at the specified time-points.

Figure 5:

An Immuno-modulatory Metabolite Risk Index (IMM-RI) is associated with overall survival, relapse, TRM and GvHD a) 2-year overall survival (OS) after Day +21 stratified according to bacterial alpha diversity. The mean alpha diversity of patient samples at Days +7–21 was calculated and patients were stratified into high and low diversity groups, defined as above or below the center-specific median Inverse Simpson’s diversity index. There were 18 deaths among 35 patients in the low-diversity group (estimated mean survival time 470 (95% confidence interval (CI) 382–557) days) and 6 deaths among 33 patients in the high-diversity group (estimated mean survival time 658 (95% CI 599–716) days). The number of patients at risk after 0, 200, 400 and 600 days after Day +21 is indicated. Analysis via Kaplan–Meier estimator, significance according to log-rank test. All hazard ratios are reported in Supp. Table 3. b) OS after a follow-up of 2 years after Day +21 stratified according to high-risk and low-risk IMM-RI, which considers optimal intestinal levels at Days +7–21 of the metabolites indicated in the text box. Optimized thresholds for metabolite concentrations were determined using the Youden Index (Supp. Table 9). Analysis via Kaplan–Meier estimator, significance by log-rank test. c) Cumulative incidence of relapse and transplant-related mortality (TRM) after 2 years after Day +21 in a competing risks analysis stratified according to high-risk and low-risk IMM-RI. Cumulative incidence functions between groups were tested for equality using the Gray’s test, the p-value is indicated. d) Cumulative incidence of GvHD and its competing risk death after 2 years after Day +21 in a competing risks analysis stratified according to high-risk and low-risk IMM-RI. Statistics as in c). e) Receiver operating characteristics (ROC) curve indicating sensitivity and specificity for the IMM-RI to predict mortality. The Area under the curve (AUC) is indicated, positive predictive value (PPV) = 0.516 and negative predictive value (NPV) = 0.895. f) Box plots of Factor values for Factors 1 and 3 (averaged across Days +7–21) in IMM-RI low- vs high-risk patients. The center line corresponds to the median, the box ranges from the 25^th to the 75^th percentile. Whisker length corresponds to the largest/lowest data point that does not exceed the 75^th/25^th percentile +/− 1.5-fold interquartile range (IQR). Blue: low IMM-RI (n=19 patients), red: high IMM-RI (n=30 patients). Statistics by two-tailed Mann-Whitney-U-test corrected for multiple testing via FDR.

Figure 6:

Microbial SCFA and butyric acid biosynthesis pathways are more abundant in IMM-RI low-risk patients a) Volcano plot of differentially abundant microbial pathways identified by metagenomic sequencing in IMM-RI low- vs high-risk patients. Blue dots indicate pathways which are significantly more abundant in IMM-RI low-risk patients, their BioCyc ID is shown. The x-axis shows fold-change (log2), the y-axis p-values (log10) adjusted for multiple comparisons via the Benjamini & Hochberg procedure. b) Waterfall plot of differentially abundant microbial pathways in IMM-RI low- vs high-risk patients identified with Maaslin2. Only features in the IMM-RI low-risk group were significant after applying abundance and p-value cutoffs (p≤0.01) and are shown. BioCyc IDs and MetaCyc pathway names are indicated. The x-axis denotes the effect size. c) Relative abundance of the specified MetaCyc pathways in IMM-RI low- vs high-risk patients. d) Summarized relative abundances within the fermentation superclass in IMM-RI low- vs high-risk patients. e) Summarized relative abundances within the SCFA and butyric acid superclasses and the indicated MetaCyc pathway in IMM-RI low- vs high-risk patients. f) Relative abundance of butyrate kinase in IMM-RI low- vs high-risk patients. g) BCoAT copy numbers (in log10) per gram of stool quantified by qPCR according to low- vs high-risk IMM-RI. h) Factor values for Factors 1 and 3 (averaged across Days +7–21) according to whether BCoAT copies were detected by qPCR (“Yes”) or not detectable (“No”). Blue: detected (n=30 patients), red: not detectable (n=36 patients). i) Detection of BCoAT copies in patients according to time-point. Bar plots show percentage, exact values are provided. The numbers of patients screened vs those in which BCoAT was detected is indicated below. j) BCoAT copy numbers (in log10) in paired patient samples at Day −7 compared to the peri-engraftment period (samples averaged across Days +7–21). k) BCoAT copy numbers in paired patients samples before (“No ABX”) and after exposure to ABX (“ABX”). For all box plots above: the box ranges from the 25th to 75th percentile. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the minimum and up to the maximum. Samples outside the 1.5-fold IQR were regarded as outliers. In h) the whisker length corresponds to the largest/lowest data point that does not exceed the 75^th/25^th percentile +/− 1.5-fold IQR. Individual patient samples are plotted as a point superimposed on the graph. The number of patients in each group is indicated in the legend. Statistics: c) to e): 1-sided Wilcoxon test and Benjamini & Hochberg correction for multiple testing; g) 2-sided Mann-Whitney test; h) two-tailed Mann-Whitney-U-test and FDR correction for multiple testing; j) 2-sided Wilcoxon test; k) 2-sided Mann-Whitney test.

Figure 7:

Detection of BCoAT auxiliary metabolic gene encoded in VC-1 and VC-2 bacteriophages a) Cladogram of bacteriophage-coded BCoAT compared against bacterial reference sequences. BCoAT was detected in 2 unique viral contigs, VC-1 and VC-2, as well as a in a Myoviridae sp. bacteriophage. Bacterial taxonomy is indicated in the legend. b) Gene alignment plot of the VC-2-coded BCoAT against a closely related reference sequence in Lawsonibacter (an Oscillospiraceae family member). The dark-grey shading indicates identity overlap in percent. A detailed plot is provided in the Extended Data Figure 8. c) Circle diagram of the VC-1 and VC-2 genomes by viral genome assembly. Genomic regions and their function are indicated in the legend. The BCoAT AMG is shown in red. d) From left to right: Relative abundance of VC-1 according to low- vs high-risk IMM-RI; detection of VC-1 in patients according to time-point; relative abundance of VC-1 in paired patients samples before (“No ABX”) and after exposure to ABX (“ABX”). Bar plots show percentage, exact values are provided. The numbers of patients screened and in which VC-1 was detected is indicated below. In the box plots, the box ranges from the 10th to 90th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the minimum and up to the maximum. Samples outside the 1.5-fold IQR were regarded as outliers. High IMM-RI: n=72 samples, low IMM-RI: n=63 samples; “No ABX”: n=40 patients, “ABX”: n=39 patients. Significance by two-tailed Mann Whitney test. e) As in c) for VC-2 f) Correlation of relative abundance of VC-1 or VC-2 with intestinal butyric acid levels (in μmol/g dry feces) by Spearman’s correlation. Axis is on a log scale. The R- and p-values are indicated in each plot. Number of pairs=130 samples.

Figure 8:

Onset of GI-GvHD and initiation of antibiotics deplete IMMs a) Intestinal bacterial and viral alpha diversity (Richness, Inverse Simpson’s diversity index) stratified by patients diagnosed with GI-GvHD (GI-GvHD) vs control allo-SCT patients (No GI-GvHD). Significance by two-tailed Wilcoxon rank sum test adjusted for multiple comparisons via the Benjamini & Hochberg procedure. The box ranges from the 25th to 75th percentiles. The line in the middle of the box is plotted at the median. The whiskers are drawn down to the minimum and maximum. Samples outside the 1.5-fold IQR were regarded as outliers. Each individual patient is plotted as a point superimposed on the graph. The number of patients per group are indicated in the legend. b) Levels of IMM-RI metabolites (in μmol per gram of dry stool) stratified by patients diagnosed with GI-GvHD vs control allo-SCT patients. Significance by two-tailed Mann-Whitney test. In the scatter-dot plots, the box is plotted at the mean. Error bars indicate standard deviation. Each individual patient is plotted as a point superimposed on the graph. The number of patients per group are indicated in the legend. c) Intestinal bacterial and viral diversity alpha diversity as in a) in paired patient samples according to antibiotic status: “No Antibiotics (No ABX)” (blue) or “Antibiotics (ABX)” (red). Once a patient was treated with antibiotics, the current and all subsequent samples were classified as “ABX”. Significance by two-tailed Wilcoxon rank sum test adjusted for multiple comparisons via the Benjamini & Hochberg procedure. Box plots as described in a). The number of patients per group are indicated in the legend. d) Level of IMM-RI metabolite levels in patient samples according to antibiotic status as in e). Significance by two-tailed Wilcoxon matched-pairs signed rank test. Scatted-dot plots as described in b). e) Box plots of Factor values for Factors 1 and 3 according to ABX (averaged across Days +7–21). Box plot as described in a), except that whisker length corresponds to the largest/lowest data point that does not exceed the 75^th/25^th percentile +/− 1.5-fold the IQR. Blue: “No ABX” (n=6 patients), red: “ABX” (n=67 patients). Significance by two-tailed Mann-Whitney-U-test corrected for multiple testing via FDR. f) Waterfall plot of differentially abundant bacterial taxa in patient samples according to antibiotic status. Only significant taxa with a p-value cutoff (p≤0.01) are shown. Blue bars indicate species that are more abundant in patients not exposed to ABX, red indicates species that are more abundant after ABX exposure. The x-axis denotes the effect size.