Integrative Transkingdom Analysis of the Gut Microbiome in Antibiotic Perturbation and Critical Illness

Abstract

Bacterial microbiota play a critical role in mediating local and systemic immunity, and shifts in these microbial communities have been linked to impaired outcomes in critical illness. Emerging data indicate that other intestinal organisms, including bacteriophages, viruses of eukaryotes, fungi, and protozoa, are closely interlinked with the bacterial microbiota and their host, yet their collective role during antibiotic perturbation and critical illness remains to be elucidated. We employed multi-omics factor analysis (MOFA) to systematically integrate the bacterial (16S rRNA), fungal (intergenic transcribed spacer 1 rRNA), and viral (virus discovery next-generation sequencing) components of the intestinal microbiota of 33 critically ill patients with and without sepsis and 13 healthy volunteers. In addition, we quantified the absolute abundances of bacteria and fungi using 16S and 18S rRNA PCRs and characterized the short-chain fatty acids (SCFAs) butyrate, acetate, and propionate using nuclear magnetic resonance spectroscopy. We observe that a loss of the anaerobic intestinal environment is directly correlated with an overgrowth of aerobic pathobionts and their corresponding bacteriophages as well as an absolute enrichment of opportunistic yeasts capable of causing invasive disease. We also observed a strong depletion of SCFAs in both disease states, which was associated with an increased absolute abundance of fungi with respect to bacteria. Therefore, these findings illustrate the complexity of transkingdom changes following disruption of the intestinal bacterial microbiome.IMPORTANCE While numerous studies have characterized antibiotic-induced disruptions of the bacterial microbiome, few studies describe how these disruptions impact the composition of other kingdoms such as viruses, fungi, and protozoa. To address this knowledge gap, we employed MOFA to systematically integrate viral, fungal, and bacterial sequence data from critically ill patients (with and without sepsis) and healthy volunteers, both prior to and following exposure to broad-spectrum antibiotics. In doing so, we show that modulation of the bacterial component of the microbiome has implications extending beyond this kingdom alone, enabling the overgrowth of potentially invasive fungi and viruses. While numerous preclinical studies have described similar findings in vitro, we confirm these observations in humans using an integrative analytic approach. These findings underscore the potential value of multi-omics data integration tools in interrogating how different components of the microbiota contribute to disease states. In addition, our findings suggest that there is value in further studying potential adjunctive therapies using anaerobic bacteria or SCFAs to reduce fungal expansion after antibiotic exposure, which could ultimately lead to improved outcomes in the intensive care unit (ICU).

Keywords: bacteria; bacteriophages; data integration; fungi; microbiome; multi-omics.

FIG 1

Overview of the composition and diversity of the bacterial, fungal, and viral microbiome. (a) Relative proportions of sequence reads at the genus level assigned to different bacterial and fungal taxa and at the order level for viral taxa. Viral metagenomics of two samples did not pass quality control due to high background levels, and these samples were therefore excluded from further analysis. For bacteria, the Ruminococcaceae, Lachnospiraceae, and Enterobacteriaceae families as well as genera that made up ≥5% of the total microbiota in at least one sample are included; other genera and families are pooled within the category “Other Bacteria.” (b) Alpha diversity metrics of bacteria (top), fungi (middle), and viruses (bottom), using the Shannon diversity index (Shannon) and the observed taxon richness index (Observed). In the box plots, the central rectangle spans the first quartile to the third quartile (the interquartile range [IQR]), the central line inside the rectangle shows the median, and whiskers above and below the box indicate variability outside the upper and lower quartiles. Given the nonparametric nature of the data, P values were calculated using the Wilcoxon rank sum test. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

FIG 2

Multi-omics factor analysis (MOFA) delineates the sources of transkingdom heterogeneity. (a) Model overview. MOFA takes as the input the three microbiome quantification matrices. MOFA exploits the covariation patterns between the features within and between microbiome modalities to learn a low-dimensional representation of the data in terms of a small number of latent factors (Z matrix) and three different weight matrices (W) (one per kingdom). By maximizing the variance explained under sparsity assumptions, MOFA provides a principled way to discover the global sources of variability in the data. For each latent factor (i.e., each source of variation), the weights provide a measure of feature importance for every feature in each factor, hence enabling the interpretation of the variation captured by every factor. (b) Heat map displaying the percentage of variance explained (R²) by each factor (rows) across the three microbe modalities (columns). Factors 1 and 3 capture coordinated variation across all three microbiome modalities, whereas factor 2 is mostly dominated by heterogeneity in fungal composition. (c) Bar plots showing the fraction of significant associations between the features of each microbiome modality and each factor. P values were obtained using a t test based on Pearson’s product-moment correlation coefficient. Statistical significance is called at a 10% FDR. This plot is useful to interpret whether the variance-explained values displayed in panel b are driven by a strong change in a small number of features or by a moderate effect across a large range of features. (d) Scatterplot of factor 1 (x axis) versus factor 3 (y axis). Each dot represents a sample, colored by condition. Factor 1 captures the gradient in microbiome variation associated with antibiotic treatment and critical illness (from negative to positive factor values), whereas factor 3 captures the variation associated with antibiotic treatment in healthy patients (positive factor 3 values) versus critically ill patients (negative factor 3 values).

FIG 3

Characterization of the transkingdom covariation captured by factor 1 and factor 3. (a) Scatterplots displaying the distribution of bacterial (top), fungal (middle), and viral (bottom) weights for factor 1. A positive value indicates a positive association with factor 1 values, whereas a negative value indicates a negative association with factor 1 values (Fig. 2d). The larger the absolute value of the weight, the stronger the association. For ease of visualization, weights are scaled from −1 to 1. Representative taxa among the top weights are labeled. (b) Heat maps displaying the reconstructed data (see Materials and Methods) based on the MOFA model for the taxa highlighted in panel a. Samples are shown in the columns (sorted based on factor 1 values), and features are shown in the rows. (c) Scatterplots displaying the distribution of bacterial (top), fungal (middle), and viral (bottom) weights for factor 3. A positive value indicates a positive association with factor 3 values, whereas a negative value indicates a negative association with factor 3 values (Fig. 2d). The larger the absolute value of the weight, the stronger the association. For ease of visualization, weights are scaled from −1 to 1. Representative taxa among the top weights are labeled. (d) Heat maps displaying the (denoised) data reconstruction (see Materials and Methods) based on the MOFA model for the taxa highlighted in panel c. Samples are shown in the columns (sorted based on factor 3 values), and features are shown in the rows.

FIG 4

Correlation of total bacterial and fungal loads with fecal levels of short-chain fatty acids in health and critical illness. (a) Association analysis between factor values and SCFA levels. (Left) Pearson correlation coefficients between factor values and the levels of three types of SCFAs: butyrate, acetate, and propionate. (Right) Corresponding FDR-adjusted and log-transformed P values. (b) Box plots showing the SCFA concentrations (y axis) per sample group (x axis). In the box plots, the central rectangle spans the first quartile to the third quartile (the interquartile range [IQR]), the central line inside the rectangle shows the median, and whiskers above and below the box indicate variability outside the upper and lower quartiles. Given the nonparametric nature of the data, P values were calculated using the Wilcoxon rank sum test. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. (c) Scatterplot of factor 1 (x axis) versus factor 3 (y axis) values. Each dot represents a sample, shaped by the sample group and colored by SCFA concentrations (in milligrams per milligram of feces). (d) Scatterplot of fungal-to-bacterial absolute level ratios (after log₁₀ transformation) (x axis) versus SCFA concentrations (after log₂ transformation) (y axis). The line represents the linear regression fit, and the shading represents the corresponding 95% confidence interval. Corresponding Pearson correlation coefficients and P values are also displayed in the top left corner.