An integrative understanding of the large metabolic shifts induced by antibiotics in critical illness

ABSTRACT

Antibiotics are commonly used in the Intensive Care Unit (ICU); however, several studies showed that the impact of antibiotics to prevent infection, multi-organ failure, and death in the ICU is less clear than their benefit on course of infection in the absence of organ dysfunction. We characterized here the compositional and metabolic changes of the gut microbiome induced by critical illness and antibiotics in a cohort of 75 individuals in conjunction with 2,180 gut microbiome samples representing 16 different diseases. We revealed an “infection-vulnerable” gut microbiome environment present only in critically ill treated with antibiotics (ICU⁺). Feeding of Caenorhabditis elegans with Bifidobacterium animalis and Lactobacillus crispatus, species that expanded in ICU⁺ patients, revealed a significant negative impact of these microbes on host viability and developmental homeostasis. These results suggest that antibiotic administration can dramatically impact essential functional activities in the gut related to immune responses more than critical illness itself, which might explain in part untoward effects of antibiotics in the critically ill.

KEYWORDS:

• Gut microbiota
• antibiotics
• critical illness
• intensive care unit
• metagenomics
• metabolomics
• ITS2

Introduction

Critical illness leads to the admission of more than 5 million patients per year to intensive care units (ICUs) in the United States alone. Intensive or invasive monitoring of ICU patients accounts for approximately 20% of the total US hospital cost, while the worldwide death rates for critically ill patients are increasing at a higher rate than any other common cause of death ¹. Almost half of ICU patients show symptoms related to an initial systemic inflammatory response syndrome (SIRS).² However, besides inflammation, signs of immune exhaustion or ‘paralysis’ might occur simultaneously.³ A disbalance of pro- and anti-inflammatory responses can lead to an increased risk of infection ⁴ and related sepsis, which are responsible for nearly 60% of deaths in ICUs and account for approximately 40% of ICU costs.⁵

Vincent et al. ⁶ reported that while only 54% of ICU patients had suspected or proven infection, as many as 70% received at least one antibiotic, reflecting a rather “liberal use” within contemporary ICUs. This reflects an early antibiotic treatment – in order to avoid the deleterious impact of a missed or delayed antibiotic therapy if infection triggers organ dysfunction.⁷ However, several studies showed that the benefits of antibiotics for prevention of infection, multi-organ failure, and death in the ICU are unclear.⁸ Evidence suggests that approximately 37% of antibiotic treatments are unnecessary or not compliant with guidelines.⁹ Infections and antibiotics can cause a rapid loss of commensal gastrointestinal microbiota,¹⁰ which can result in metabolic and immune disturbances in the critically ill. An important role of gut bacteria is the fermentation of dietary fiber into short-chain fatty acids (SCFAs), which play an important role in preserving gut integrity.¹¹ A decrease in SCFAs concentrations may result in colonization by pathogenic species. For example, several Bacteroides and Bifidobacteria species secrete SCFAs inhibiting pathogen growth.^12,13 Similarly, Clostridium scindens and Ruminococcus obeum produce secondary bile acids (BAs) that prevent the growth of Clostridium difficile ¹⁴ and Vibrio cholerae,¹⁵ respectively. Gut bacteria also play an important role in the activation of host immunity against infections, both through innate and adaptive mechanisms.¹⁶

With advent of the concept of the human being as a “holobiont” and the perception of the gut microbiome as being highly relevant in the regulation of the immune system, attention to interventions affecting the microbiome is now given also in critical illness. Previous studies have focused predominantly on taxonomic information using 16S rRNA gene sequencing to identify differences between health and disease.^17,18

Fungal constituents of the microbiome represent an overlooked but very important kingdom. Research is beginning to show that fungi are critical for maintaining systemic immunity and intestinal homeostasis.¹⁹ The mycobiome of skin,²⁰ gut,²¹ oral cavity,²² and lungs,²³ among other anatomical sites, in healthy individuals has been characterized in previous studies. Based on these studies, it seems that between individuals and anatomic sites there is high variability in the human mycobiome diversity,²⁴ which is consistent with what we know from the Human Microbiome Project (HMP, 2012) for the bacterial microbiome. Most of the anatomic sites in humans are dominated by members of the Basidiomycota phylum, such as Cryptococcus spp., Malassezia spp. and Filobasidium spp., and the Ascomycota phylum including Saccharomyces cerevisiae, Candida spp., and Cladosporium spp.^25,26 Despite the recent findings showcasing the fungal role on host health, host–microbe and microbe–microbe interactions, only less than 0.5% of the published microbiome papers investigate or refer to the fungal population.¹⁹ The mycobiota is increasingly recognized as a critical player in the development of human diseases, including inflammatory bowel disease, allergic airway diseases, skin disease, alcoholic liver disease, autoimmunity, neurological disorders, and metabolic syndrome.^27–30 In relation to microbe–microbe interactions, existing studies indicate that a competitive association between bacteria and fungi exists in the human gut. This was shown in antibiotic-treated subjects and germ-free mice, where an overgrowth of particular fungi in the gut and/or susceptibility to fungal infection was observed.³¹ In addition, overgrowth of fungi due to antibiotics treatment has been associated with the development of allergic airway responses to Aspergillus fumigatus mold spores.³² Part of the chemical warfare between bacteria and fungi is also the secretion of antifungal peptides from epithelial cells, which can be induced by commensal bacteria such as Blautia producta and Bacteroides thetaiotaomicron.³³ Previous research related to critically ill patients and mycobiome is very limited.^34–39 One recent study characterized the mycobiome of the lower respiratory tract of patients in ICU showing that Candida spp. dominated the fungal community in both with and without antibiotic therapy patient groups.³⁴ Another prospective pilot study showed an increase of C. albicans in the oral mycobiome after an admission to the ICU.³⁸ However, systematic investigation of the interactions between the gut mycobiome of critically ill patients and the bacteria functional activity is currently lacking.

Here, by a comprehensive characterization of the microbiome, mycobiome, and functional potential of the gut community and individual species, we demonstrate that even though antibiotics do not significantly disturb the bacterial and fungal composition of critically ill patients, as observed in healthy individuals,⁴⁰ they cause abundance changes in a handful of species that are highly connected with the production of SCFAs and BAs, allowing the expansion of pathogenic species.

Results

Highly distinct microbiome in ICU patients

To investigate the gut microbiome composition of critically ill patients, we initially included 70 critically ill patients. Of these, 54 were diagnosed with probable or microbiologically confirmed infection as defined by Calandra et al. ⁴¹ (respiratory tract (n = 37), abdominal (n = 6), bones/soft tissue (n = 3), chest (n = 2), catheter associated (n = 1), urogenital (n = 1), and unknown (n = 4)). These 54 received antibiotic treatment, whereas 16 did not presented an infection and did not receive antibiotic treatment. Of those, we collected stool samples from 49 patients receiving antibiotic treatment (broad spectrum beta-lactam antibiotics; n = 19 meropenem and n = 30 piperacillin/tazobactam) and 14 patients without antibiotic treatment, as well as from 12 healthy human volunteers. Basic anthropometric and clinical characteristics of the participants are displayed in Supplementary File 1. There were no significant differences between the two ICU groups in gender, age, BMI, type of admission or surgery and length of ICU or hospital stay (continuous data were compared by the Student’s t-test, dichotomous variables by the chi-squared test, a p-value < 0.05 was considered significant). Similarly, there were no significant differences between the healthy volunteers and ICU patients in the basic demographic characteristics.

We assessed the structure of the gut microbiome via shotgun metagenomics generating 39 million high-quality reads per sample on average (Supplementary Table 1). We used MetaPhlan2 ⁴² for taxonomic profiling, and we identified 428 species in total. Alpha-diversity measured as Shannon, Simpson, and Chao1 index dropped significantly in both antibiotics treated (ICU⁺) (Wilcoxon rank-sum test, P = 8.6e-6, P = 5.3e-6, P = 3.8e-5, Shannon, Simpson and Chao1 indices, respectively) and untreated subjects (ICU^–) (Wilcoxon rank-sum test, P = .0037, P = .0025, P = .0054, Shannon, Simpson and Chao1 indices, respectively) compared to healthy individuals (Figure 1(a)). In contrast to what has been observed in healthy volunteers,^44,45 antibiotics administration had no significant impact on the alpha diversity when comparing the ICU⁺ against the ICU^– patients (Wilcoxon rank-sum test, P = .79, P = .79, P = .64, Shannon, Simpson and Chao1 indices, respectively). The type of beta-lactam used had no significant impact either (Supplementary Figure 1).

Figure 1. Distinct gut microbiota signatures in ICU patients

(a – b, d) Box plots showing the median (centerlines), first and third quartiles (box limits) and 1.5x interquartile range (whiskers) measurements. A comparison was considered significant if P < .05.(a) Alpha diversity of bacterial species using Shannon (left), Simpson (middle), and Chao1 (right) indices. Significant differences were determined using Wilcoxon rank-sum test.(b) Beta dispersion of bacterial species measured as the distance of the samples from one group to the group centroid in multivariate space. Significant differences were determined using Wilcoxon rank-sum test.(c) Principal component analysis (PCoA) of Bray-Curtis dissimilarity between bacterial species abundance profiles. Significant differences were determined using PERMANOVA and were considered significant if P < .05.(d) FEAST ⁴³ estimation of microbial source contribution for each “sink”. Here, sinks are species level relative abundances from samples from a publicly available ICU cohort.¹⁰ For sources, we used species level relative abundances of our Healthy, ICU^– and ICU⁺ groups. Significant differences were determined using Kruskal-Wallis test (*P < .05; **P < .01, ***P < .001, ****P < .0001).(e) Receiver operating characteristic curve (ROC) of a cross-validated random forest classifier. The model was trained on taxonomic and functional profiles from 63 samples from the ICU^– and ICU⁺ groups. Model performance was summarized as area under ROC (AUC). The average AUC value and confusion matrix (threshold 0.5) are calculated based on a 5-fold cross validation results.

We subsequently investigated the variation in the microbiome structure of the three groups by calculating the distance to centroid (Figure 1(b)). The gut microbiome of the ICU⁺ group had significantly higher structural variation compared to both healthy individuals (Wilcoxon rank-sum test, P = 1.1e-12) and the ICU^– group even though it did not reach statistical significance (Wilcoxon rank-sum test, P = .056). Next, we calculated the beta-diversity (Bray–Curtis distance) of the three groups, which showed that although there was a high community dissimilarity between ICU⁺ and Healthy (PERMANOVA, P = .001) and between ICU^– and Healthy (PERMANOVA, P = .001), the ICU⁺ and ICU^– groups did not show significant differences (PERMANOVA, P = .104) (Figure 1(c), Supplementary Figure 9(a)). Similarly, to alpha diversity, when we examined the two beta-lactam groups (meropenem and piperacillin/tazobactam) separately, we did not observe any significant difference (Supplementary Figure 1, Supplementary Figure 9(b)). To further evaluate our findings, we repeated the taxonomic analysis using the Metagenomic Species method (MGS) ⁴⁶ and we observed the same patterns as above (Supplementary Figure 2, Supplementary Figure 9(c)).

As evident from both, the alpha- and beta-diversity comparisons, the similarities in the structure of the gut microbiome of the ICU⁺ compared to the ICU^– group triggered our interest to examine our cohort in relation to another, small publicly available dataset of critically ill patients.¹⁰ In that study, the authors performed a longitudinal sampling of critically ill patients with and without antibiotics administration to study abundance changes in the resistome profile of known pathogens. For the comparative analysis, we applied a microbial source tracking algorithm, namely FEAST.⁴³ When using our cohort as a sink we noticed how progressively the public ICU⁺ group becomes more taxonomically similar with our ICU⁺ group, as well as the compositional similarities between the two ICU^– groups (Figure 1(d)). From the methodological point of view, this analysis confirms that FEAST ⁴³ is a relatively sensitive method but also that despite the non-significant overall differences between ICU⁺ and ICU^– groups there is still a characteristic microbiome signature due to the antibiotic administration. Considering the relatively small differences between the ICU⁺ and ICU^– groups, we next investigated whether these subtle differences in microbial signatures could be integrated into an algorithm to correctly classify the patients. To this end, a random forest classifier integrating species and pathways was developed and achieved an area under the receiver operating characteristic (ROC) curve (AUC) of 0.934 (Figure 1(e)) using 20 features (Supplementary Figure 3). Enterococcus faecalis, a known pathogen⁴⁸ whose abundance has been associated to an increase in susceptibility to V. cholerae infection,¹⁵ and several amino acid pathways that can serve as precursors for the synthesis of short-chain fatty acids (SCFAs) were among the selected features.

Critically ill patients treated with antibiotics show an “infection-vulnerable” gut microbiome

We subsequently used HUMAnN2 ⁴⁹ for functional profiling of the microbial communities and identified 483 pathways in total. To further evaluate the community taxonomic and functional characteristics of the critically ill patients, we compared our cohort with a panel of publicly available datasets. For that, we retrieved 2,180 taxonomically and functionally annotated samples from 16 human diseases using the curatedMetagenomicData ⁵⁰ package in R and we applied FEAST,⁴³ treating our samples as sources and the disease panel as sinks. In these settings, FEAST ⁴³ revealed that there is a significantly lower similarity (Kruskal–Wallis, P < .0001) between the microbiome structure and function of individuals in the critically ill patients group with patients with cancer (colorectal cancer,^51–56 melanoma^57,58), inflammatory (ankylosing spondylitis,⁵⁹ atherosclerosis,⁶⁰ Behcet’s disease,⁶¹ hypertension,⁶² inflammatory bowel disease⁴⁶), and metabolic diseases (cirrhosis,⁶³ metabolic syndrome,⁶⁴ nonalcoholic fatty-liver disease,⁶⁵ obesity,^66,67 Type-1-Diabetes,^68,69 Type-2-Diabetes^70,71) compared to Healthy individuals (Figure 2(a)). Due to this unexpected association, we expanded the range of diseases as sinks using infectious diseases and specifically 66 samples from patients with acute diarrhea,⁷² C. difficile infection⁷³ and V. cholerae infection.⁷⁴ In contrast to what we observed with noninfectious diseases, the taxonomic composition of the critically ill patients is significantly more similar (Kruskal–Wallis, P < .0001) to infectious diseases than the gut composition of healthy individuals (Figure 2(b)). At the functional level, this difference is even more profound with the ICU⁺ group having the highest similarity (Kruskal–Wallis, P < .0001) with all three infectious diseases, suggesting that critically ill patients may be at a risk of severe infections by gut pathogens and that antibiotic treatment increases significantly that risk.

Figure 2. Comparative analysis of the microbiome of critically ill patients with other diseases

(a – b) FEAST ⁴³ estimation of microbial source contribution for each “sink”. Here, sinks are taxonomic and functional compositions from a panel of diseases. Source are taxonomic and functional compositions of samples from our Healthy, ICU^– and ICU⁺ groups. Box plots show the median (centerlines), first and third quartiles (box limits) and 1.5x interquartile range (whiskers), source contributions to sinks. (a) cancer, inflammatory and metabolic diseases and (b) infectious diseases. (a – b) Significant differences were determined using Kruskal-Wallis test (*P < .05; **P < .01, ***P < .001, ****P < .0001).

Immunomodulatory metabolites and their microbial producers are depleted from critically ill patients exposed to antibiotics

We subsequently focused on the impact of antibiotic administration in the ICU in relation to abundance changes in microbes with a potential role in immune regulation and host-immune homeostasis. Significant differences were observed at the phylum level among the three groups. The Healthy group had, on average, higher abundance of Bacteroidetes (17.9%, Kruskal–Wallis, P = .0063), Actinobacteria (9.9%, Kruskal–Wallis, P = .035), and Verrucomicrobia (4.7%, Kruskal–Wallis, P = .0098) (Supplementary Figure 4(a)). At the genus level, a striking difference in the abundance of Enterococcus was found, with the ICU⁺ group having the highest abundance compared to the other two groups (Kruskal–Wallis, P = 3.6e-6). In contrast, the relative abundance of Blautia was significantly higher in ICU^– (Kruskal–Wallis, P = .029) (Supplementary Figure 4(a)).

At the species level, we found 106 and 80 species significantly differentially abundant (Wilcoxon rank-sum test, P < .05) between the Healthy and ICU⁺ and Healthy and ICU^– groups, respectively. The species enriched in the Healthy group include known SCFAs and BAs producers, such as Ruminococcus bromii,⁷⁶ Faecalibacterium prausnitzii,⁷⁷ Eubacterium eligens,⁷⁸ Eubacterium hallii ⁷⁸ and Eubacterium rectale,^77,78 among others (Figure 3(a) and Supplementary Figure 4(b)). On the contrary, species enriched in ICU⁺ and ICU^– included known pathogens such as Klebsiella pneumoniae,⁷⁹ Klebsiella oxytoca,⁷⁹ E. faecalis,⁴⁸ and Enterococcus faecium ⁴⁸ (Figure 3(a) and Supplementary Figure 4(b)). Comparing ICU⁺ and ICU^– groups, we found 12 species significantly differentially abundant (8 enriched in ICU^– and 4 in the ICU⁺). R. gnavus, Clostridium symbiosum, and Veillonella parvula, known SCFAs and BAs producers,^80–83 were enriched in ICU^–, whereas from the species enriched in ICU⁺ (Figure 3(a)) Bifidobacterium animalis has been indicated as a SCFAs producer.⁸⁴ Since B. animalis was the only SCFA producer enriched in ICU⁺, we retrieved the genome-scale metabolic model of B. animalis from the AGORA repository ⁸⁵ and we simulated growth on an ICU media using flux balance analysis. We monitored the potential of B. animalis to produce SCFAs, including acetate, propionate, and butyrate; however, none of these compounds was predicted to be produced in our in silico simulations (data not shown).

Figure 3. Differential abundance analysis reveals changes in bacteria with important functional properties in critically ill

(a) Taxonomic trees visualized using R package metacoder.⁴⁷ Only taxa differentially abundant between Healthy, ICU⁺ or ICU^– (P < .05, Wilcoxon rank-sum test) are highlighted in the tree by color. Color of the taxa reflects the group with higher abundance. Bar plots show the relative abundances of significantly different (FDR<0.05, Wilcoxon rank-sum test) short-chain fatty acid producers (green circles), bile acid producers (red circles), disease-associated species (purple circles) or differentially abundant species identified in all pairwise comparisons (yellow circles).(b) Co-abundance networks of differentially abundant bacterial species (P < .05, Wilcoxon rank-sum test) between Healthy and ICU⁺ and between ICU⁺ and ICU^–. Only significant correlations (P < .05) with absolute correlation coefficient >0.4 were used for network construction. Nodes are colored based on their affiliated phyla. Node sizes reflect the mean abundance of the species. Edge colors reflect either negative correlation (blue) or positive correlation (red). Edge widths reflect the strength of the correlation.(c-d) Germ-free L1 larval stage C. elegans worms were populated with depicted bacterial strains in the anoxic chamber for indicated times, followed by transfer to normoxia and UV-killed OP50 E. coli diet. In (C) host survival was measured after 24 h of normoxic culture and in (D – left panel) the timely development was assessed after 48 h of normoxic culture; (D – right panel) depicts reproductive aging of nematodes following anoxic reconstitution with B. animalis.

We used the enriched species to reconstruct the species co-abundance network using all samples from our cohort. We observed a much more intense within group communication between the Healthy-enriched species (vs. ICU⁺) with 162 positive and 19 negative correlations (absolute correlation coefficient >0.4) compared to the species enriched in ICU⁺ (vs Healthy and/or ICU^–) that showed 43 positive and no negative correlations (Figure 3(b) and Supplementary Table 2). Alistipes putredinis, F. prausnitzii, and Bacteroides uniformis were interconnected healthy-enriched species showing negative correlations with Klebsiella and Staphylococcus species, common sources of serious infections in ICU.⁸⁶ From the species enriched in ICU⁺ compared to Healthy, E. faecium, Staphylococcus epidermidis, and Lactobacillus rhamnosus showed the highest number of correlations, suggesting that these may be important species to maintain the community structure. From the species enriched in ICU⁺ compared to ICU^–, Lactobacillus crispatus showed positive correlations with known pathogens such as Enterococcus and Klebsiella species, as well as negative correlations with Healthy-enriched SCFAs and BAs producers including F. prausnitzii, E. rectale, and A. shahii.

To examine in vivo host effects of the species enriched in ICU⁺ patients compared to ICU^–, we utilized nematode C. elegans as a model host. C. elegans recently emerged as a suitable simple model for discovering conserved host–microbiome interactions.⁸⁷ The experiments were conducted by reconstituting germ free L1 nematodes with the two, sufficiently annotated, bacterial strains enriched in ICU⁺ (B. animalis and L. crispatus). We found that both strains had a strong negative impact on the host homeostasis: L. crispatus exerted direct toxicity in the model host (Figure 3(c)), and B. animalis instigated a delay in the C. elegans development indicative of the physiological stress (Figure 3(d)). By including heat-killed bacteria as an additional control and performing OD normalization across conditions, we found that host toxicity of L. crispatus requires live bacteria, while B. animalis rather acts as a passive stressor (Supplementary Figure 5). Interestingly, the nematodes, which overcame the developmental hindrance by B. animalis, displayed a delay in reproductive aging (Figure 3(d)) consistent with the putative probiotic effect of B. animalis in the animals that were able to conquer the initial stress caused by this bacterium. To probe the potential probiotic effect of B. animalis at the mechanistic level, we tested its ability to induce nuclear translocation of DAF-16/FOXO transcription factor – a conserved mediator of stress resistance and longevity extension.^90–92 We found that live B. animalis indeed had the strongest capacity to induce DAF-16 activation among all conditions tested (Supplementary Figure 6). Our results are thus in line with the previously reported ability of B. animalis to cause disease in immunocompromised human patients, while it acts as a probiotic in healthy humans,⁹³ demonstrating the physiological relevance of our nematode findings. Collectively, our in vivo tests indicate that the rarefication of the microbiome, which is exacerbated by antibiotics exposure in ICU⁺ patients, facilitates the enrichment of microbes with potential of exerting direct detrimental effects on the host.

Functional shifts in the microbiome

Next, we compared MetaCyc pathway abundance to explore the gut microbiome functionality in the three groups. As shown in the ordination plot (Figure 4(a)), the differences between ICU⁺ compared to the other 2 groups at the functional level are driven by SCFA biosynthesis, AA biosynthesis, and fermentation to SCFA, whereas the abundance of pathways related to nucleotide and nucleoside degradation and vitamin biosynthesis are the main drivers in the comparison of Healthy vs ICU^–. Therefore, we then focused on AAs, SCFAs, and BAs due to their important role in the regulation of the immune system and their influence in diseases associated with dysbiosis.¹¹ In total, we identified 117 significantly differentially abundant pathways related to the metabolism of AAs, SCFAs, and BAs (Wilcoxon rank-sum test, P < .05). Even though the majority of these pathways had a higher relative abundance in the Healthy compared with ICU^– group, only a few pathways reached statistical significance. On the contrary, the differences in the relative abundances between the ICU⁺ with the Healthy group were striking (Figure 4(b)). The relative abundance of 98 of the 117 unique pathways was significantly lower in ICU⁺. Despite the small number of significantly differentially abundant species between ICU⁺ and ICU^–, there were more differences in the metabolic pathways between the two groups than between Healthy and ICU^–, suggesting that antibiotic treatment has a stronger effect in the metabolism of AAs, SCFAs, and BAs than critical illness. Several of these pathways were selected in the random forest model to classify the two groups of ICU patients (Figure 1(e) and Supplementary Figure 3).

Figure 4. Functional shifts contributed by the microbiome

(a) Multidimensional scaling (MDS) plot of the samples based on the pathway abundances. The top 15 pathways with the strongest significant correlation with the overall ordination (FDR<0.05 with function envfit from the R package vegan ⁷⁵) are highlighted with arrows where the length of the arrows reflects the strength of the association.(b) Numbered heatmaps showing the relative abundances of differentially abundant bacterial MetaCyc pathways. Additional heatmaps indicate significance (blue tiles, P < .05, Wilcoxon rank-sum test). Only pathways related to short-chain fatty acid, bile acid and amino acid metabolism, as identified by manual curation, are shown.

Metabolomics analysis

We then performed targeted metabolomic analysis and quantified the levels of 10 SCFAs and 27 BAs in 38 stool samples (see Methods). Consistently with the results from the MetaCyc pathway analysis, we observed a significant decrease in the abundance of SCFAs and BAs in the ICU⁺ group compared to the Healthy and ICU^– groups (Figure 5(a)). Of 10 identified SCFAs, 6 were significantly different between the groups (Figure 5(a)). Acetic acid, propionic acid, butyric acid, and valeric acid were significantly lower in the ICU⁺ group compared to both the Healthy and ICU^– groups (Wilcoxon rank-sum test, P < .05). The levels of these SCFAs were not found to differ significantly between the ICU^– and Healthy groups. A similar pattern was observed with the BAs; ketolithocholic acid, deoxycholic acid, glycolithocholic acid, hyodeoxycholic acid, isolithocholic acid, lithocholic acid, and ursodeoxycholic acid were all found significantly lower in abundance in ICU⁺ patients (Wilcoxon rank-sum test, P < .05) compared to both Healthy and ICU^– patients (Figure 5(a)). Among them, ursodeoxycholic acid is increasingly used in the clinical setting for a treatment of a variety of conditions.^94,95 Ketolithocholic acid, lithocholic acid, and ursodeoxycholic acid have been also found to provide resistance against C. difficile infections and to modulate the host inflammatory response during the infection.⁹⁶ The BA profiling of ICU⁺ patients may also explain the high similarity at the functional level revealed by FEAST with the C. difficile cohort (Figure 2(b)).

Figure 5. Associations between gut microbiome and short-chain fatty acids and bile acids

(a) Content of SCFAs and BAs (log_e) in three groups are plotted as bar plots. Wilcoxon rank-sum test, *P < .05, **P < .01, ***P < .001, ****P < .0001, ns: not significant.(b) Heatmap showing Spearman’s correlations between significantly different metabolites and previously identified significant species. The species selected were either having (i) significantly higher abundance in the Healthy group compared to ICU⁺ or (ii) significantly higher abundance in the ICU^– compared to the ICU⁺ group. Metabolites category is shown with row annotation and the origin of enriched species is shown with column annotation. Correlations with FDR<0.05 are marked with a circle. Heatmap in the lower panel shows the secretion of short-chain fatty acids obtained by flux variability analysis (FVA) ⁸⁸ of available bacterial metabolic models. Unavailable models of corresponding species are shown in gray. Bottom boxplots show the results from the Growth Rate Analysis (GRiD).⁸⁹ The red dashed line indicates GRiD ⁸⁹ value 1.

To identify which species were mainly responsible for the differences in the levels of the measured SCFAs and BAs in ICU⁺ compared to the other two groups, we performed growth rate analysis ⁸⁹ and Spearman’s correlation between the 13 metabolites and 38 species. The species selected were either having (i) significantly higher abundance in the Healthy group compared to ICU⁺ but not compared to the ICU^– group or (ii) significantly higher abundance in the ICU^– compared to the ICU⁺ group (Figure 5(b)). As shown in Figure 5b, A. putredinis, Lachnospiraceae bacterium 2 1 58FAA and Lachnospiraceae bacterium 1 1 57FAA showed positive correlations with SCFAs and/or BAs and were predicted through Flux Variability Analysis (FVA) ⁸⁸ to secrete SCFAs. Several Bacteroides species, including B. uniformis, were actively growing (GRiD ⁸⁹ > 1), they showed positive correlations with the measured SCFAs and they were predicted through FVA ⁸⁸ to secrete SCFAs supporting the identified correlation (Figure 5(b)). Since BA metabolism is not included in these metabolic models, we analyzed the BA biosynthesis potential of the species showing high correlation with BAs using differential analysis of enzymes involved in the BA biosynthesis pathway. The enzymes cbh and baiN were found in the genomes of Bacteroides xylanisolvens, L. bacterium 2 1 58FAA, R. obeum, and R. gnavus. The abundance levels of B. xylanisolvens and R. obeum cbh and the abundance levels of R. obeum baiN were found significantly higher in the Healthy group compared to the ICU⁺ (Wilcoxon rank-sum test, P < .05, data not shown). R. obeum has been previously shown to have an inhibitory effect on V. cholerae due to its capacity to degrade virulence-activating signals in the gut through the synthesis of bile salt hydrolases.¹⁵

Candida species and resistant genes flourish in critically ill patients

Since the microbiome composition and functional profile of critically ill patients appeared to be dysbiotic, we investigated next if this leads to systematic changes in the structure of the fungal community, with emphasis on Candida species and the antibiotic resistance gene (ARG) levels in the three groups. We built high-quality libraries for ITS2 sequencing of 74 available stool samples from the Healthy, ICU⁺ and ICU^– groups. We estimated the fungal relative abundance using the DADA2 pipeline.⁹⁷ The taxonomic profiling revealed that Ascomycota (85.63%) was the most abundant fungal phyla, followed by Basidiomycota (9.84%) and Muromycota (3.44%) (Figure 6(a)). There was no significant difference in the alpha diversity (Shannon, Simpson, Chao1) between any of the two groups and the type of beta-lactam administrated in the ICU⁺ group did not appear to influence this pattern (Wilcoxon rank-sum test, P > .05) (Supplementary Figure 7). However, when we calculated the Bray–Curtis distance of the three groups we observed a high community dissimilarity between ICU⁺ and Healthy (PERMANOVA, P = .001) and between ICU^– and Healthy (PERMANOVA, P = .001), whereas the ICU⁺ and ICU^– groups did not show significant differences (PERMANOVA, P = .703) (Figure 6(a), Supplementary Figure 10(a)). We subsequently investigated differences in fungal species relative abundance between the three groups. In total 19, 11, and 3 fungal species were found significantly different (Wilcoxon rank-sum test, P < .05) in the comparisons between Healthy vs ICU⁺, Healthy vs ICU^–, ICU⁺ vs ICU^–, respectively (Figure 6(b)). The Candida genus, which includes species that are opportunistic pathogens, was found to have the greatest number of differentially abundant species between the three groups. Candida albicans, Candida glabrata, Candida pseudolambica, and Candida tropicalis were all found to have the highest abundance in ICU⁺ patients (Figure 6(b)).

Figure 6. ICU patients have distinct mycobiome profiles

(a) Upper – PCoA plots of the gut mycobiome based on species-level Bray-Curtis distances. Color indicates the groups. Significant differences were determined using PERMANOVA (P < .05). Lower – the same PCoA plots as upper panel highlighting the top 3 most abundant fungal phyla (Ascomycota, Basidiomycota and Mortierellomycotina). Color scale indicates the relative abundance of each phylum. Groups are defined by shapes.(b) Heatmap showing the normalized abundances of significantly differentially abundant species in Healthy vs ICU⁺, Healthy vs ICU^–, ICU⁺ vs ICU^– comparisons (P < .05, Wilcoxon rank-sum test). Cumulative Sum Scaling (CSS) was used for species abundances normalization. The annotation bars on the left indicate the comparisons in which the species was found to be significant. The names of species are colored according to their affiliated genera. Differentially abundant species in ICU⁺ vs Healthy comparison after FDR correction (FDR<0.25) are marked with *.(c) Boxplots comparing the normalized abundances of specific Candida species (Wilcoxon rank-sum test).

To determine whether antibiotic treatment exerts selective pressure on the resistome as a whole, we analyzed the change in Pfams related to the resistome and mobilome, as well as the abundance differences of ARGs between the three groups. In the abundance comparison between the Healthy and the ICU⁺ groups, there were 71 statistically significant Pfams (Wilcoxon rank-sum test, P < .05) related to the resistome and/or mobilome, with 48 of them being more abundant in the ICU⁺ and 23 in the Healthy group (Supplementary Figure 8(a)). Interestingly, the differences between ICU^– and Healthy groups were also large with 30 and 16 Pfams being more abundant in the ICU^– and the Healthy group, respectively (Supplementary Figure 8(a)). In the comparison between ICU⁺ and ICU^– there were 19 significant Pfams with 14 being more abundant in the ICU⁺ and 5 in the ICU^– group (Supplementary Figure 8(a)).

Subsequently, we annotated the ARGs in the three groups using deepARG.⁹⁸ The overall ARG abundance profile of the ICU⁺ and ICU^– groups using the Bray–Curtis distances indicates that there is no substantial perturbation during antibiotic treatment (PERMANOVA, P = .261), but there are significant differences between ICU⁺ and Healthy (PERMANOVA, P = .001) and between ICU^– and Healthy (PERMANOVA, P = .001) (Supplementary Figure 8(b), Supplementary Figure 10(b)). Comparing the total accumulative relative ARGs abundance led to the same conclusion; critically ill patients have already a unique resistome profile compared to healthy individuals, which is only marginally disturbed by antibiotics administration (Wilcoxon rank-sum test, P = 5.97e-8, Supplementary Figure 8(c)). Considering the polypharmacology approaches often applied in ICU, our observation is in agreement to recent studies suggesting that human targeted drugs can significantly impact on the gut resistome profile.⁹⁹ Nevertheless, an abundance comparison of the individual ARGs revealed 24 ARGs, including 3 beta-lactams (SHV, PENA, AMPC), significantly higher in the ICU⁺ group and only 1 in the ICU^– group (dabestr, Confidence Interval (CI) = 95%, Supplementary Table 3). Among the ARGs that exhibit the highest abundance differences (dabestr, 95%) between the ICU⁺ and ICU^– groups were AAC(6ʹ)-I (aminoglycoside), ADEC (multidrug), ERMB (macrolide, lincosamide and streptogramin), VANS (glycopeptide), and MSRC (macrolide, lincosamide and streptogramin) (Supplementary Figure 8(d)).

Discussion

Sepsis, acute respiratory distress syndrome (ARDS), and multi-organ failure represent common conditions most frequently driven by an inappropriate host response to pathogens of the critically ill. These conditions are responsible for immense global mortality accompanied by a tremendous economic burden ¹. While breakthroughs of molecular medicine have revolutionized treatment in oncology and rheumatology, in critical illness research endeavors of decades have not resulted in any targeted therapies. In practice, intensive care can be considered mostly as supportive and antibiotics are a cornerstone of care for patients with sepsis, i.e. infection-driven organ dysfunction.⁷

The off-target effects of antibiotics on the microbiome are, however, also particularly obvious in the ICU population. Antibiotics often fail to resolve organ failure despite evidence of infection and, even more concerning, are frequently administered not to miss an occult infection. While the resulting increase of multi-resistant bacteria is an obvious problem, the negative impact on the “holobiont” in the ICU have largely been ignored. More to the point, considerations in intensive care are more dwelling around the early days of introduction of antimicrobial therapy when Paul Ehrlich propagated the concept of “therapia sterilisans magna” where only “parasitotropic“ effects in the absence of “organotropic“ effects of drugs were envisioned.¹⁰⁰ Thus, molecular therapies for these common and lethal diseases are desperately needed and depend on a better understanding of systems biology of the host metabolome–microbiome interplay.¹⁰¹ Up to now, the downstream consequences, such as host inflammation and cellular damage, and not the upstream sources, in particular the complex microbial ecosystems that reside in and on the human body, have been the priority of research. Nevertheless, two recent studies in mice ^102,103 indicate that Fecal Material Transplantation (FMT) and specific species in the gut microbiome could prevent sepsis opening up new clinical research avenues.

Here, we present evidence from a human study that the selective pressures to which critically ill patients are exposed (parenteral nutrition, polypharmacy, including e.g., proton pump inhibitors, shock states requiring invasive life support measures, such as catecholamine treatment) shape the microbiome of these patients in a unique way with highly distinct characteristics compared to healthy or other disease states, including metabolic, inflammatory, or malignant diseases. As we have shown, the microbiome structure and function of critically ill patients resembles signatures mainly observed in severe infections such as C. difficile and V. cholerae. This already “infection-vulnerable” microbiome structure in critical illness becomes severely dysbiotic after antibiotic treatment with an observed depletion of SCFAs, including propionate, butyrate, and acetate, and BAs. Similarly, in an elegant mice study Kim et al., ¹⁰² demonstrated using FMT from healthy littermates that high levels of butyrate (and potentially also propionate), provided mainly by Bacteroidetes, can rescue from lethal sepsis caused by a pathogenic mixture of K. oxytoca, E. faecalis, Serratia marcescens and C. albicans isolated from a septic patient. Three of these pathogenic species, K. oxytoca, E. faecalis and C. albicans, were found in significantly higher abundance in our critically ill patients compared to healthy individuals. However, the poor similarity in gut microbial taxonomic abundances between human and mice ¹⁰⁴ highlights the importance of investigating the relationship between critical illness and gut microbiota in human clinical samples. For example, while Kim et al., ¹⁰² hypothesized that critical illness itself may result in the depletion of gut butyrate. In our human study, this was not the case and only after antibiotics administration we observed a dramatic decrease in the SCFA levels. Nevertheless, we observed a significant change in taxonomy also in ICU^– compared to Healthy related to inflammation, such as a significant increase in the abundance of Staphylococcus and Enterococcus, which have been reported as key factors for the development of signs of systemic inflammation, nosocomial infection and complications in the ICU course.^48,86

Our study has several limitations. Critically ill patients represent a heterogeneous patient group characterized by comorbidities, past infections and age, leakage of alveolar, and intestinal barriers as well as impaired defense and repair mechanisms.^105,106 Low flow states up to overt shock, autonomic dysfunction, and lack of sufficient options for enteral feeding, e.g., due to dysphagia or impaired consciousness affect gut function, transit time, and defecation.¹⁰⁷ Many drugs that are applied routinely in addition to antibiotics, such as proton pump inhibitors or catecholamines can affect gut function substantially.⁹⁹ Nevertheless, while rodent studies allow to control many of the aforementioned confounders, only clinical studies allow to address the impact of antibiotics on the gut microbiome as it relates to clinically meaningful outcomes. Thus, we aimed at control of confounders through inclusion of patient cohorts requiring critical care but not receiving antibiotics.

Patients in ICUs represent a relatively small subgroup of hospitalized patients, but they reflect a specific at-risk population that accounts for approximately 25% of all hospital-acquired infections.¹⁰⁸ Infection as a typical complication of critical illness increases length of ICU stay, morbidity, mortality, and costs.¹⁰⁹ Moreover, the rise in multidrug-resistant strains prompts attention on nonantibiotic strategies in the prevention and treatment of nosocomial infections, such as pro- or synbiotics. However, mechanisms of action in the vulnerable population of critically ill patients requires further investigation into the mechanisms that shape the gut microbiome. From our study we suggest that the loss of a handful of species, that are highly connected with the production of SCFAs and BAs, during antibiotic administration in the ICU allows the expansion of pathogenic species, which exhibit potential to cause direct hindrance of host homeostasis. Despite the availability of more advanced antibiotics in ICUs death rates from sepsis following nosocomial infections keep increasing, indicating that these antibiotics do not increase survival but instead they produce a highly dysbiotic gut ecosystem that allows more aggressively resistant and lethal pathogens to thrive.¹ These changes are likely associated with or even to promote a state of “protracted critical illness”, a frequent observation in patients discharged after prolonged intensive care and characterized by persistent systemic infection. It seems warranted to design studies that aim to restore the gut microbiome or replace key metabolites, such as SCFAs or BAs, in this vulnerable patient population to restore homeostasis of the “metaorganism” after discharge from intensive care.

Materials and methods

Study design

A prospective observational study was undertaken from May 2018 until January 2019 at the Jena University Hospital. Adult critically ill patients either treated with systemic antimicrobial therapy (piperacillin/tazobactam or meropenem) for at least 2 days or without any systemic antimicrobial therapy within the last 7 days were eligible for this study. Patients with inflammatory bowel disease, major bowel resection, selective decontamination of the oral and digestive tract, oral vancomycin therapy, immunocompromised patients, history of chemotherapy during the last 6 months, or known travel history to areas of high antimicrobial resistance within the last 4 weeks were excluded from this study. The need for informed consent was waived by the IRB, since this investigation did not involve any intervention at the patient. The collected basic data were used only for calculating mean values per group and were not linked to the individual metagenomic profiles. Stool samples from healthy volunteers (>50 years old) with no antimicrobial therapy within the last 6 months served as a control group. Healthy volunteers interested in participating were invited for a consultation with medical doctors in the Jena University Hospital. Their health status and prior antibiotic use was self-reported. The basic demographic data for the healthy cohort is provided in Supplementary File 1. We collected up to 100 ml feces which were sampled immediately after natural defecation. Fecal specimens were transferred into two sterile containers, one that was mixed with liquid thioglycolate medium supplemented with catalase and 10% glycerol and one without buffer solution. Both containers were stored at –80°C. Each patient contributed with only one stool sample.

DNA extraction from stool samples

All stool samples were processed by Novogene (UK). DNA was extracted using the following protocol: Stool samples were thoroughly mixed with 900 μL of CTAB lysis buffer. All samples were incubated at 65°C for 60 min before being centrifuged at 12000 × g for 5 min at 4°C. Supernatants were transferred to fresh 2-mL microcentrifuge tubes and 900 μL of phenol:chloroform:isoamyl alcohol (25:24:1, pH = 6.7; Sigma-Aldrich) was added for each extraction. Samples were mixed thoroughly prior to being incubated at room temperature for 10 min. Phase separation occurred by centrifugation at 12,000 × g for 15 min at 4°C, and the upper aqueous phase was re-extracted with a further 900 μL of phenol:chloroform:isoamyl alcohol. Next, samples were centrifuged at 12,000 × g for 10 min at 4°C, and the upper aqueous phases were transferred to fresh 2-mL microcentrifuge tubes. The final extraction was performed with 900 μL of chloroform:isoamyl alcohol (24:1), and layer separation occurred by centrifugation at 12,000 × g for 15 min at 4°C. Precipitation of DNA was achieved by adding the upper phase from the last extraction step to 450 μL of isopropanol (Sigma-Aldrich) containing 50 μL of 7.5 M ammonium acetate (Fisher). Samples were incubated at – 20°C overnight, although shorter incubations (1 h) produced lower DNA yields. Samples were centrifuged at 7500 × g for 10 min at 4°C, and supernatants were discarded. Finally, DNA pellets were washed three times in 1 mL of 70% (v/v) ethanol (Fisher). The final pellet was air-dried and re-suspended in 200 μL of 75 mM TE buffer (pH = 8.0; Sigma-Aldrich).

Library preparation and sequencing for metagenomics

Sequencing library was generated based on Illumina technologies and following manufactures’ recommendations. Index codes were added to each sample. Briefly, the genomic DNA was randomly fragmented to a size of 350 bp, then DNA fragments were narrowly size selected with sample purification beads. The selected fragments were then end polished, A-tailed, and ligated with adapter. These fragments were filtered with beads again and amplified by PCR reaction. At last, the library was analyzed for size distribution and quantified using real-time PCR. The library was then to be sequenced on an Illumina platform Novaseq 6000 (Novogene) with paired-end reads of 150 bp.

Internal transcribed spacer sequencing

The concentration of genomic DNA was determined by Qubit, and the DNA quality was checked on the gel. 200 ng of DNA was used as input for PCR reaction with corresponding primer set specifically binding to different hypervariable regions. Each primer set had a unique barcode. PCR product was then run on the gel and DNA fragment with the proper amplification size was cut and purified. The purified PCR product was then used as template for library preparation. The PCR products were pooled together with equal amount and then end polished, A-tailed, and ligated with the adapter. These fragments were filtered with beads again. After PCR reaction (to make library fully double strand), the library was analyzed for size distribution and quantified using real-time PCR. The library was then to be sequenced on Hiseq2500.

Metabolomics analysis

We performed targeted metabolomics analysis for 38 of the 75 available samples. The remaining samples were destroyed during a prolonged stay in the customs during the COVID-19 pandemic.

Quantification of SCFAs: SCFAs were extracted by addition of 2 mg ultra-pure water pr. mg of sample. The samples were vortex mixed for 1–2 min until suspension is reached, and centrifuged at max speed for 10 min at 4°C. The supernatant was transferred to a spinX centrifuge filter and centrifuged for additional 5 min at 4°C. The filtrate was collected and stored at –20°C until analysis. Sample analysis was carried out by MS-Omics as follows. Samples were acidified using hydrochloric acid, and deuterium labeled internal standards where added. All samples were analyzed in a randomized order. Analysis was performed using a high polarity column (ZebronTM ZB-FFAP, GC Cap. Column 30 m x 0.25 mm x 0.25 μm) installed in a GC (7890B, Agilent) coupled with a quadrupole detector (5977B, Agilent). The system was controlled by ChemStation (Agilent). Raw data was converted to netCDF format using Chemstation (Agilent), before the data was imported and processed in Matlab R2014b (MathWorks, Inc.) using the PARADISe software described by Johnsen et al.¹¹⁰

Quantification of BAs: Bile acids were extracted by addition of 4 mg methanol pr. mg of sample. The samples were vortex mixed for 1–2 min until suspension is reached, and centrifuged at max speed for 10 min at 4°C. The supernatant was transferred to a spinX centrifuge filter and centrifuged for additional 5 min at 4°C. In a HPLC vial, 285 μL filtrate is mixed with 15 μL internal standard. The final volume of filtrate of CS10768 and CS10798 where below 285 μl. Therefore, 200 μl were combined with 85 μl Methanol and 15 μl internal standard. Sample analysis was carried out by MS-Omics as follows. The analysis was carried out using a Thermo Scientific Vanquish LC coupled to Thermo Q Exactive HF MS. An electrospray ionization interface was used as ionization source. Analysis was performed in negative ionization mode. The chromatographic separation of bile acids was carried out on a Waters Acquity HSS T3 1.8 μm 2.1 × 150 mm (Waters). The column was thermostated at 30°C. The mobile phases consisted of (A) ammonium acetate 10 mmol/l, and (B) methanol:acetonitrile (1:1, v/v). Bile acids were eluted by increasing B in A from 45% to 100% for 16 min. Flow rate was 0.3 min. Peak areas were extracted using Tracefinder 4.1 (Thermo Scientific). Identification of compounds were based on accurate mass and retention time of authentic standards.

Data processing

Quality control of sequence data

Quality control to remove low-quality reads was performed as described previously.¹¹¹ Briefly, all Illumina primer/adapter/linker sequences were removed. Subsequently, low-quality regions (consecutive regions with Phred quality <20) were trimmed. Finally, all reads were mapped to the human genome with BWA version 0.7.4 ¹¹² and reads with >95% identity and 90% coverage were removed as human DNA contamination.

Taxonomic profiling

Taxonomic annotation of the high-quality reads was performed using MetaPhlAn2⁴² version 2.7.7 with default settings, generating taxonomic relative abundances. Bacterial community profiles were constructed at phylum, genus and species level for further analyses.

Taxonomic annotation of fungal ITS was performed using the DADA2 pipeline⁹⁷ version 1.14 with default parameters including adapter removal, quality filtering and trimming, dereplication of identical reads, read-pair merging, ITS2 extraction and chimera removal. Remaining reads were binned as operational taxonomic units and aligned to the UNITE fungi database using RDP classifier.¹¹³ All samples were then normalized by cumulative sum scaling using R package metagenomeSeq.¹¹⁴

Functional annotation

The HUMAnN2 pipeline ⁴⁹ version 0.11.2 was used for functional annotation of the high-quality reads after the quality control. The quantified pathway and gene family abundances in the units of RPKs (read per kilobase) were then normalized to copies per million (CPM) units by the provided HUMAnN2 script, resulting in transcript-per-million-like (TPM) normalization. Gene families were then regrouped to Pfam domains for further analyses.

Abundance comparisons

Species, pathways, and Pfams were filtered by 10% prevalence across all samples and their relative abundances were used for statistical comparisons between the three groups. Differentially abundant features were identified by the Wilcoxon rank-sum test and were considered significantly differentially abundant if the P-value was less than 0.05.

Differentially abundant phyla and genera were identified by the Kruskal-Wallis test and were considered significantly differentially abundant if the P-value was less than 0.05.

Metagenomics sequences from HUMAnN2 profiles

Gene family abundances were clustered using mgs-canopy⁴⁶ version 1.0 software with standard parameters. Gene family clusters were considered metagenomic sequences (MGS) if they had at least 700 genes. Taxonomic annotation of MGS was done using species annotation information available for each gene family.

We calculated contributions of each species to an MGS. An MGS was annotated to the species with the largest contribution if: the gene contribution of that species was more than 50% and the second largest species was “unclassified” or contributed less than 10%.

Diversity analysis

Alpha diversity indices Shannon, Simpson, and Chao1 were calculated using the R packages vegan⁷⁵ and fossil¹¹⁵ based on relative species abundance. Wilcoxon rank-sum test was used to test for significant differences in alpha diversity. For estimating community dissimilarities, Bray–Curtis distances were calculated using the R package vegan⁷⁵ based on the relative species abundance. To test for significant differences in the microbial composition, permutational multivariate analysis of variance (PERMANOVA), as implemented in the function adonis from R package vegan,⁷⁵ was used to analyze beta-diversity.

Co-abundance networks

The relative abundance table for significantly different species was processed using SparCC ¹¹⁶ for co-abundance network inference. Species–species correlation coefficients were estimated as the average of 20 inference iterations and 100 permutations were used for the pseudo P-value calculation. For the visualization of the co-abundance network, only interactions with an absolute correlation coefficient >0.4 were used.

Metabolic modeling

To estimate the availability and composition of metabolites in ICU patients, the nutrition fed in ICU (https://www.fresenius-kabi.com/de/ernaehrung/fresubin-original) was considered. Metabolic composition of complex products such as fish oil was described by vmh diet designer (https://www.vmh.life/#nutrition/dietdesigner). The human genome-scale metabolic model Recon3D 3.01 ¹¹⁷ simulated based on the ICU specific diet (Supplementary File 2) was used to predict metabolites that can potentially be secreted by the host. Flux Variability Analysis (FVA) ⁸⁸ was used to determine feasible exchange reaction flux bounds that support metabolite secretion alongside optimal growth rate. Identified metabolites were assumed to be available for the bacterial species and strains to be consumed. Genome-scale metabolic models of the studied species and strains were collected from two different gut model repositories, AGORA 1.03 (https://www.vmh.life) ⁸⁵ and CarveMe.¹¹⁸ Taking into account the availability of ICU diet compounds and potential host-secreted metabolites in the gut, the maximum amount of SCFAs (acetic acid, propionic acid, and butyric acid) production by bacterial species and strains were predicted by applying FVA ⁸⁸ again alongside achieving maximum ATP yield as objective for the available bacterial metabolic models.

Abundance of ARGs

The metagenomic reads were analyzed using the deepARG pipeline,⁹⁸ which uses deep learning to identify and quantify ARGs. Reads were compared to the provided DeepARG-DB database using a prediction model to evaluate sequence similarities and predict antibiotic resistance. The pipeline was run in short sequence mode with a minimum probability cutoff of 0.8, an identity cutoff of 80%, an E-value cutoff of 1e-10 and a minimum coverage of 50%.

Testing for significant differences in ARG abundance was performed using R package dabestr ¹¹⁹ with a confidence interval of 95%.

Bacterial growth rate estimation

Bacterial growth rate was calculated using the growth rate index (GRiD) ⁸⁹ version 1.2. The algorithm first calculates the coverage of all contigs of a reference genome in the sample, sorts them from high to low, and reorders them to two groups, placing an ori-containing contig at start and a ter-containing contig at the mid-region of the genome. Next, it calculates coverage drops across a sliding 10 Kb window, with values representing the coverage ratio of the peak and trough of the curve. High values represent faster growth rates.

Random forest model

A Random Forest classifier was built to classify ICU patients into ICU^– and ICU⁺ based on bacterial taxonomic profiles and pathways. The model was implemented using R package caret ¹²⁰ with all bacterial species and pathway abundances as input features. The model was trained after centering and scaling the data and removing variables with near zero variance, using a tune length of 10 and fivefold cross-validation as resampling method, the rest of the parameters were left as default. Feature importance were calculated using function varImp from R package caret.¹²⁰ A random forest model was then built using only the 20 most important features. R package PRROC ¹²¹ was used for ROC calculation and plot.

Bacterial exposure and survival assay in C. elegans

Bacterial strains used were Escherichia coli OP50, Lactobacillus crispatus (DSM 20356), and Bifidobacterium animalis (DSM 20104). All strains except E. coli were acquired from the German Collection of Microorganisms and Cell Cultures and cultured following the supplier guidelines. Bacterial stocks were kept at –80°C until use. C. elegans strain used was N2 wild-type isolate obtained from the Caenorhabditis Genetics Center.

Bacterial stocks were inoculated in anoxic broth (MRS for L. crispatus and BSM for B. animalis) and incubated at 37°C for 72 h. 150 µL of each bacterial culture was spread onto medium sized NGM agar plates (5,5 cm diameter) and incubated for 24 h at 37°C in an anaerobic container (BD GasPak™ EZ container systems) prior to worm addition. E. coli was grown on NGM in normoxia. UV-killed OP50 was produced by exposing NGM plates to UV light (320 nm) in a Chemi-Doc XRS+ transilluminator (BioRad) for 10 min.

Age-synchronized germ-free worms were obtained by collecting eggs from gravid adults upon treatment with alkaline hypochlorite solution (composition per liter: 200 mL 1 M KOH, 250 mL bleach, and 500 mL ddH2O) and letting the eggs hatch overnight at 20°C in M9 buffer. Synchronized L1 larvae were seeded onto NGM plates containing respective bacterial strains (4 plates per strain). Approximately 150 worms were seeded per plate. Plates were incubated at 20°C in an anaerobic container for 5 h or 24 h. After incubation worms were washed with M9 buffer and transferred to UV-killed OP50 plates. UV-killed OP50 plates were incubated at 20°C in normoxic conditions for 24 h prior to survival assessment. Heat-killed bacteria was obtained by submerging bacterial cultures at 80°C for 60 min and OD600 was normalized to 0,2 before placing the cultures onto NGM plates. Live bacteria was also subjected to OD600 normalization in all tests, which involved heat-killed control conditions.

Survival rate was assessed by screening all the worms present in each plate. A worm was regarded as dead if it did not respond to gentle touch with a platinum wire. Survival was expressed as percentage of the total number worms. Each experiment was performed 3 times.

Developmental fitness and reproductive aging assays in C. elegans

Bifidobacterium animalis was grown for 72 h on anoxic BSM broth at 37°C. Afterward, 150 uL of bacterial culture was spread onto medium-sized NGM plates (4 plates) and incubated for 24 h at 37°C in an anaerobic container (BD GasPak™ EZ container systems) prior to worm addition. E. coli OP50 was grown on NGM at normoxic conditions. C. elegans N2 population was synchronized as described above and approximately 150 worms were seeded on either B. animalis or OP50 plates. Plates were incubated for 24 h at 20°C in the anaerobic container. Developmental assay was carried out as previously described.¹²² In brief, 30 worms per bacterium were put individually onto small UV-killed OP50 plates and incubated for 48 h at 20°C before developmental stage of each worm was visually assessed. Reproductive aging assay was carried out as described previously.¹²³ In brief, after incubation with B. animalis or OP50 the worms were washed with M9 and let to develop until L4 stage on UV-killed OP50 plates at 20°C, normoxia. At this moment 25 randomly picked worms (per condition) were transferred individually onto small-sized UV-killed OP50 plates. Every day the brood size of each worm was determined (sum of eggs and L1s) and parent worm was transferred to new plate until egg laying ceased. These experiments were performed 3 times.

DAF-16 nuclear translocation assay

Nematodes expressing DAF-16::GFP fusion protein were obtained from the Caenorhabditis Genetics Center (strain TJ356). Bacterial strains used were the same as described above. Anoxic broth was inoculated with an aliquot of actively growing bacterial culture and incubated at 37°C for 48 h. OD600 of all living cultures was normalized to 0,1. 150 µL of each bacterial culture was spread onto medium sized NGM agar plates (5,5 cm diameter) and incubated for 24 h at 37°C in an anaerobic container (BD GasPak™ EZ container systems) prior to worm addition. E. coli was grown on NGM in normoxia. Heat-killed bacteria was generated as described above and OD600 was normalized to 0,2 prior to seeding onto NGM plates.

Age-synchronized germ-free worms were obtained as described above and grown until L4 stage on UV-killed OP50. L4 worms were washed with M9 buffer and transferred to NGM plates containing bacteria. Approximately 150 worms were seeded per plate. Plates were incubated at 20°C in an anaerobic container (BD GasPak™ EZ container systems) on either live or heat-killed bacteria for 5 h. Control plates were incubated in normoxia on live bacteria only. After incubation, 30 worms were picked from each condition and transferred to empty NGM plates for imaging. Imaging was carried out using a ZEISS Axio Zoom.V16 microscope equipped with fluorescence light. Imaged worms were sorted into three categories (nuclear, intermediate, cytosolic) depending on the localization of the GFP tagged DAF-16 transcription factor. This experiment was performed three times.

Statistical analysis

To determine differential abundance of taxonomic, functional, and metabolic features between groups Wilcoxon two-tailed rank-sum test was applied when analyzing the differences between two groups, whereas Kruskal–Wallis test was used when more than two groups were compared. Correlation between microbial taxa and metabolites was assessed by Spearman’s correlation. The R package dabestr¹¹⁹ was used to test differential abundance of ARGs. Significant differences in source contributions to sinks using FEAST⁴³ were assessed using Wilcoxon two-tailed rank sum test. To assess differences in alpha diversity and beta dispersion, Wilcoxon two-tailed rank sum test was used, whereas PERMANOVA was used for beta diversity.

Ethics approval and consent to participate

The study was approved by the local ethics committee in Jena, Germany (5289-10/17). Informed consent was waived by the ethics committee because feces were only collected after natural defecation and no patient data were obtained.

Author contributions

GP and MB designed the study. MB was involved in sample collection and processing. GP designed and supervised the bioinformatics analysis and edited the manuscript. AMS and LZ carried out the data analysis and data integration. MM performed the metabolic modelling. BS and YN helped with the interpretation of results and revised the manuscript. RS carried out the quality control and pre-processing of the metagenomic data, and together with AB and CB helped with the interpretation of results. PAP, MM and ME designed and performed the animal studies. All authors read and approved the final manuscript.

Disclosure statement

The authors declare that they have no competing interests.

Availability of data and material

Raw sequence files were deposited into the NCBI Sequence Read Archive (SRA) under accession number PRJNA670323.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

GP, MB, AMS, AB, ME, and LZ would like to thank Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) CRC/Transregio 124 ‘Pathogenic fungi and their human host: Networks of interaction’, subproject B5 and INF and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Germany’s Excellence Strategy - EXC 2051 – Project-ID 390713860 for financial support. LZ and CB have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Project: 802736 MORPHEUS).