Metagenomic assessment of gut microbial communities and risk of severe COVID-19

Background The gut microbiome is a critical modulator of host immunity and is linked to the immune response to respiratory viral infections. However, few studies have gone beyond describing broad compositional alterations in severe COVID-19, defined as acute respiratory or other organ failure. Methods We profiled 127 hospitalized patients with COVID-19 (n = 79 with severe COVID-19 and 48 with moderate) who collectively provided 241 stool samples from April 2020 to May 2021 to identify links between COVID-19 severity and gut microbial taxa, their biochemical pathways, and stool metabolites. Results Forty-eight species were associated with severe disease after accounting for antibiotic use, age, sex, and various comorbidities. These included significant in-hospital depletions of Fusicatenibacter saccharivorans and Roseburia hominis, each previously linked to post-acute COVID syndrome or “long COVID,” suggesting these microbes may serve as early biomarkers for the eventual development of long COVID. A random forest classifier achieved excellent performance when tasked with classifying whether stool was obtained from patients with severe vs. moderate COVID-19, a finding that was externally validated in an independent cohort. Dedicated network analyses demonstrated fragile microbial ecology in severe disease, characterized by fracturing of clusters and reduced negative selection. We also observed shifts in predicted stool metabolite pools, implicating perturbed bile acid metabolism in severe disease. Conclusions Here, we show that the gut microbiome differentiates individuals with a more severe disease course after infection with COVID-19 and offer several tractable and biologically plausible mechanisms through which gut microbial communities may influence COVID-19 disease course. Further studies are needed to expand upon these observations to better leverage the gut microbiome as a potential biomarker for disease severity and as a target for therapeutic intervention. Supplementary Information The online version contains supplementary material available at 10.1186/s13073-023-01202-6.

observations to better leverage the gut microbiome as a potential biomarker for disease severity and as a target for therapeutic intervention.
Most prior studies have largely focused on the presence, absence, or the differential abundance of specific microbes in COVID-19 [7, 9-16, 19, 20], and few have interrogated microbial network dynamics to identify which co-occurring or co-excluded species are foundational to maintaining microbial homeostasis. This represents a missed opportunity to identify potential bacterial targets to restore a more favorable, health-promoting gut configuration. Similarly, other studies have not considered how these shifts might influence gut metabolite pools. Finally, prior studies interested in exploring the gut microbiome in COVID- 19 have largely sought to characterize the differences in healthy controls compared to infected patients rather than those with moderate compared to severe disease [7, 10-12, 14, 16]. Establishing a predictive biomarker of disease severity may improve early identification of at-risk patient populations that require immediate intervention or those that are more likely to benefit from effective antiviral therapies [21].
It remains unclear what role the gut microbiome plays in regulating the severity of COVID-19 in hospitalized patients and what specific microbially-mediated mechanisms may underlie this relationship. To address these questions, we conducted a study of hospitalized patients with COVID-19 at a US tertiary medical center. Using metagenomic profiling of fecal samples collected from these patients, we demonstrate significant depletions of Fusicatenibacter saccharivorans and Roseburia hominis in severe COVID-19, reductions of which have previously been linked to post-acute COVID-19 syndrome (PASC) or long COVID [18,22]. Strikingly, we observed these declines during patients' index hospitalizations, suggesting the presence of an early microbial signal that may predict the development of a long-term complication. We further use network analysis to identify significant changes in microbial co-occurrence networks in severe COVID-19 and perform complementary predicted metabolite analyses to further link these changes to alterations in bile acid pool and short-chain fatty acid (SCFA) levels, offering biologically plausible mechanisms to explain the link between gut microbial communities and COVID-19 disease severity.

Study population
From April 2020 to May 2021, we prospectively enrolled 127 consecutive hospitalized patients aged ≥ 18 years with confirmed COVID-19 at the Massachusetts General Hospital to a longitudinal COVID-19 disease surveillance study. Patients were categorized as having severe COVID-19 if they required admission to the intensive care unit with acute respiratory failure (the need for oxygen supplementation ≥ 15 L per minute (LPM), non-invasive positive pressure ventilation, or mechanical ventilation) or other organ failures (such as shock requiring vasopressor initiation) [23]. Otherwise, they were categorized as having moderate COVID-19. Patients were screened daily for inclusion from among all admitted individuals for whom a designation of possible SARS-CoV-2 infection was flagged by hospital infection control. COVID-19 infection status was subsequently confirmed with at least one positive nasopharyngeal SARS-CoV-2 polymerase chain reaction (PCR) test. An optional biospecimen collection protocol was nested within this longitudinal study, which allowed collection of additional clinically relevant biospecimens, including stool samples. All consecutive eligible consenting patients were included.

Sample/data collection
Fresh stool was collected and refrigerated at 4℃ until aliquoting/freezing at − 80℃ (typically within 4 h of collection) from adult patients enrolled in the prospective biospecimen collection study (241 samples from 127 admitted patients). Participants were able to provide stool samples as frequently as once daily, as well as declining donation on any given day (while remaining in the study). Study coordinators blinded to case status abstracted data from the electronic health record using a double data entry approach with discrepancies adjudicated by re-abstraction or after discussions with supervising authors. We collected information on admission age (years), biological sex (male, female), race (White, Black, Asian, American Indian, Mixed, or Other), ethnicity (non-Hispanic or Hispanic), admission BMI (kg/m 2 ), comorbidities including history of cancer, pulmonary, or cardiac disease, hypertension, hyperlipidemia, and diabetes mellitus (each yes/ no), smoking history (active, former, never, unknown, and pack-years among smokers), and their composite admission Charlson Comorbidity Index, a validated score predictive of in-hospital mortality [24]. Information on hospital course, including admission Simplified Acute Physiology Score II (SAPS II) [25] and Sequential Organ Failure Assessment (SOFA) scores [26] were calculated from routine laboratory results and clinical assessments. The use of antibiotics, antivirals including remdesivir, hydroxychloroquine, corticosteroids, anti-IL-6 therapy, any form of oxygen support, highflow oxygen, bilevel positive airway pressure (BiPAP) ventilation, or mechanical ventilation (each yes/no) was collected. Mortality within 90 days of admission was ascertained in the post-study period.

Extraction protocols
Stool samples, reagent-only negative controls, and mock community positive controls (Zymo Research) were extracted using either the AllPrep PowerFecal DNA/ RNA 96 Kit (Qiagen) or the Maxwell HT 96 gDNA Blood Isolation System (Promega) [27]. SARS-CoV-2 viral load was quantified as per CDC guidelines [28] using the 2019-nCoV N1 primer and probe set [28], as well as human RNaseP as an internal control. Each RT-qPCR reaction contained TaqPath ™ 1-Step RT-qPCR Master Mix (Thermo Fisher), RNA template, the CDC N1 or RNaseP forward and reverse primers (IDT), probe, and RNase-free water to a total reaction volume of 10 μl. Viral copy numbers were quantified using N1 quantitative PCR (qPCR) standards (IDT) in tenfold dilutions to generate a standard curve. The assay was run in triplicate for each sample with three no-template control wells per 384 well plate.

Microbial sequencing
Samples were sequenced by two metagenomic sequencing facilities at the Broad Institute and Baylor College of Medicine according to their standard established platforms. DNA was prepared for sequencing using the Illumina Nextera XT DNA library preparation kit. All libraries were sequenced with a target of 3 GB output at 2 × 150 bp read length using the Illumina NovaSeq platform.

Sequence bioinformatics
Taxonomic and functional profiles from both locally recruited patients and publicly available sequences from our external validation cohort [19] were generated using the bioBakery 3 shotgun metagenome workflow 3.0.0, the details of which have previously been described [29]. Briefly, human reads were filtered using KneadData 0.10.0 and taxonomic profiles generated using MetaPhlAn 3.0.0 [30]. Functional profiling was conducted using HUMAnN 3.0.0 [30], resulting in gene family abundance tables assembled into higher order MetaCyc pathways [31].
Given the tight coupling and relatively conserved nature of gut taxonomic and metabolite profiles [32], we used the MelonnPan-predict 0.99.023 workflow [33] to interrogate the functional relationship between COVID-19 severity and microbial community metabolism. In brief, MelonnPan uses an elastic net model to conservatively predict putative metabolite levels based on stool UniRef90 gene family abundance.

Statistical analysis
To compare patient characteristics between study groups, we used standard statistical tests, including chisquared (χ 2 ) tests or Fisher's exact testing for categorical variables, the Student's t-test for normally distributed, non-categorical variables, and nonparametric Wilcoxon rank sum tests for all others. Differences with two-tailed p-value ≤ 0.05 were considered significant.ɑ-diversity was calculated using the Shannon index with the "diversity" function from the R package vegan [34]. Principal coordinates analyses (PCoA) were performed using specieslevel Bray-Curtis dissimilarity metrics with the "vegdist" function in the vegan package.
After filtering out features with no variance and low (< 10%) prevalence, we performed differential abundance testing of species-level taxonomy, MetaCyc pathways, and predicted stool metabolites using linear mixedeffects models to account for a nested data structure from repeated sampling of non-independent samples: Machine learning model building and evaluation were conducted using the SIAMCAT v. 1 to remove biomarkers with low overall abundance and z-transformed. A nested cross-validation procedure was applied to calculate prediction accuracy by splitting data into training and testing sets for twice-repeated, fivefold cross-validation. To account for longitudinal sampling [35], data splits were stratified by participant ID, ensuring samples from the same individual were used in the same fold. For each split, a random forest (RF) regressor was trained and subsequently used to predict COVID-19 disease severity. To evaluate model performance, we used the lambda parameter to maximize the area under the receiver operator characteristic curve (AUROC) with a 95% confidence interval (CI) for cross-validation error. We used the make.predictions function of SIAMCAT to assess model performance on our external validation dataset.
To assess whether ecological dynamics may help explain observed differences in taxonomy, we performed dedicated microbial network analyses. To account for our longitudinal data structure and the non-independence of longitudinal samples from the same individual, we restricted this analysis to each participant's first collected stool (all other analyses used the entire dataset). Network construction was conducted using the "netConstruct" function in NetCoMi v.1.0.2 [36], normalized using a modified centered-log ratio and limited the resulting network to microbes with an absolute Pearson correlation ≥ 0.4 (approximately equal to the 95th percentile of correlation matrix distribution). Network hubs were identified as those in the top quintile of degree, betweenness, and closeness centrality in each network (moderate vs. severe COVID-19, respectively). Finally, comparison of moderate and severe networks was performed using the "netCompare" function with 10,000 permutations.

Participant characteristics and overall gut community structure
We enrolled 127 hospitalized COVID-19 patients. 79 (62.2%) had severe disease and 48 (37.8%) had moderate disease. Collectively, they provided 241 stool samples (Fig. 1a, Additional file 1: Fig. S1). While BMI was higher in the severe group, there were no statistically significant differences observed between severity groups based on age, sex, race, ethnicity, various comorbidities, and smoking history (Additional file 1: Table S1). Patients with severe COVID-19 had a higher mean body mass index (BMI) as well as Simplified Acute Physiology Score II (SAPS II) [25] and Sequential Organ Failure Assessment (SOFA) scores [26], each a validated clinical assessment tool to risk stratify hospitalized patients' risk of mortality [37,38]. Severe COVID-19 patients more frequently received antibiotics, antivirals, and ICU therapies. Patients with severe COVID-19 had higher 90-day mortality compared to those with moderate disease (22.8% vs. 4.2%, p-value = 0.01).
Gut microbial diversity was significantly reduced in severe COVID-19 after adjusting for factors such as recent antibiotic use (Fig. 1b, p-value < 0.0001). We found that COVID-19 disease severity explained a statistically significant proportion of variance in Bray-Curtis distances (4.04%, FDR p-value = 0.01), while other demographic factors and details related to hospital course had either a modest and/or non-statistically significant effect on overall community structure-this finding was not fully explained by characteristic trade-offs along the Bacteroidetes/Firmicutes axes of variation [39] or prior antibiotic usage (Fig. 1c, d). No major batch effects attributable to sequencing center were observed, and thus, subsequent analyses were conducted on pooled samples (multivariable PERMANOVA R 2 for batch = 1.2%, p-value = 0.12, Additional file 1: Fig. S2).

Differential abundance testing
Using multivariable linear mixed-effects modeling accounting for SARS-CoV-2 stool viral load, (which has previously been linked to increased COVID-19-related mortality [40]), age, sex, antibiotic use, race/ethnicity, and other relevant clinical metadata (Methods), we observed statistically significant differences in 48 species-level taxa between severe and moderate COVID-19 (FDRcorrected p-value < 0.05, Fig. 2a and Additional file 1: Table S2). All but two of these taxa (Candida albicans & Enterococcus faecalis) were relatively depleted in severe disease ( Fig. 2a, b), a trend concordant with the observed decrease in species richness and evenness. While not directly comparable, the highest absolute β-coefficients from our multivariable modeling for antibiotic use was 3, while 27 of 48 significant taxonomic associations demonstrated coefficients >|3|, suggesting a consistently stronger link between COVID-19 severity and alterations in relative microbial abundance than antimicrobial therapy (Additional file 1: Table S2). We identified significant depletions of Fusicatenibacter saccharivorans and Roseburia hominis (Fig. 2b), consistent with prior work showing the relative contraction of each in patients with post-acute COVID-19 syndrome (PASC), also known as "long COVID" [18,22]. The abundance of F. saccharivorans and R. hominis were not significantly associated with clinical factors included in our model, including age, sex, and BMI (all FDR p-values > 0.05), though there was a trend towards increased E. faecalis with greater time between hospital admission and sample collection (FDR p-value 0.052; Additional file 1: Table S2).
Eight taxa were positively associated with stool SARS-CoV-2 viral load, including Methanobrevibacter smithii and Bilophila wadsworthia, as well as several Alistipes spp (Additional file 1: Table S2). Interestingly, an expansion of R. hominis was associated with increased stool viral load (Additional file 1: Table S2). Corresponding to community-wide depletions in microbial diversity, biochemical pathways encoded by gut bacteria were also significantly altered in severe COVID-19, including reductions in amino acid biosynthesis (e.g., glutamine synthesis), isoprenoid biosynthesis, and short-chain fatty acid production (SCFA) pathways, including glycerol degradation, acetyl-CoA fermentation, and methanogenesis from acetate (Additional file 1: Table S3 and Additional file 1: Fig. S3).
To ensure the robustness of our findings, we performed a sensitivity analysis in which we performed our multivariable differential abundance testing on stool collected within 30 days of admission (the median length of stay). We showed that other than an anticipated loss of power from decreased sample size, our findings were not materially altered. Of the 48 differentially abundant species in our primary analysis, 32 remained significant with this more stringent criteria. When similarly restricting our analysis to samples preceding the use of antibiotics (if any), 35 of the 48 differentially abundant species remained statistically significant (Additional file 1: Table S4).

Accurate classification of COVID-19 severity using a microbiome-based random forest learner
Given our findings of both community-wide and featurelevel alterations linked to severe COVID-19, we next Study overview and overall community structure. a Study enrollment of hospitalized patients with confirmed COVID-19 with weekly stool sampling until the time of discharge or death, whichever occurred first. b Marked reduction in species richness and evenness in severe COVID-19 (inverse Simpson ɑ-diversity metric, p-value < 0.0001 from multivariable linear modeling adjusting for age, sex, prior antibiotic use, race, ethnicity, body mass index, Charlson Comorbidity Index, use of remdesivir or corticosteroids, days since admission, SARS-CoV-2 stool viral load, sequencing depth, and a participant-level random effect). Boxes represent median and interquartile range, while whiskers represent 95%ile. c Community-level disturbances in severe vs. moderate COVID-19 as depicted by joint ordination and principal coordinates analysis (PCoA), not fully explained by characteristic trade-offs in Bacteroidetes/Firmicutes or prior antibiotic use. d Ominibus testing of Bray-Curtis distances demonstrates that COVID-19 severity had a modest and statistically significant impact on the overall community structure. Other demographic information, covariates, and hospital course information were not significantly associated (FDR p-value > 0.05) used a machine learner to predict whether metagenomic features could serve as inputs to classify samples derived from patients with severe vs. moderate COVID-19. To assess whether non-microbial metadata (i.e., participant characteristics) should be jointly considered with microbial taxa in training our classifier, we generated an entropy heatmap to quantify the unique row-wise information with respect to column-wise data (in which noninformative variables would have a value of 0). As all the covariates used in our prior linear modeling (Methods) contributed unique information to label/disease severity prediction (Additional file 1: Fig. S4), each was included in our machine learning workflow.
Using both differentially abundant microbial features and clinical characteristics as our input with five-fold twice-repeated cross-validation (Fig. 3a), our random forest regressor achieved an area under the receiver operating characteristic (AUROC) of 0.925 when tasked with predicting whether stool was obtained from patients with severe vs. moderate COVID-19 (Fig. 3b). Our findings were only modestly attenuated when modeled without clinical metadata (AUROC 0.922) and stool SARS-CoV-2 viral load (AUROC 0.923), respectively. To robustly assess this result, we trained our model using only the top 20 differentially abundant microbial features, which only modestly degraded task performance (AUROC 0.898). Finally, though we ensured samples from the same individual were confined to a single cross-fold, to minimize the possibility of overfitting data from personalized gut microbial communities, we trained and assessed our model using only the first stool sample from each participant, which again performed with excellent accuracy (AUROC 0.871), further supporting the role of metagenomic profiling as a diagnostic biomarker for disease severity.    Table S2. b Highlighted box and scatter plots of taxa abundance by COVID-19 severity. For visualization purposes, technical/true 0s were imputed with a given taxa's minimum non-zero value. Boxes represent median and interquartile ranges, while whiskers represent 95%ile A comprehensive literature review identified one metagenomic cohort with publicly available information on COVID-19 disease severity [19]. After uniform preprocessing of raw sequences (Methods), we tested our model on this external dataset of 38 patients with mild/ moderate vs. severe/critical COVID-19. Despite heterogeneity in case definition, collection methods, country of origin, and the lack of additional clinical metadata beyond disease severity, testing our taxa-only classifier achieved an AUROC of 0.741 on this external dataset, independently verifying a strong association between COVID-19 disease severity and alterations in gut microbial communities (Fig. 3c).

Systems approaches to interrogate microbial assemblages
To explore the possible biological mechanisms underlying our observations, we next sought to compare microbial co-occurrence networks in moderate vs. severe COVID-19 disease (Methods). We hypothesized that the community-wide and feature-level alterations observed in moderate vs. severe COVID-19 would change microbial network topology. First, we evaluated global microbial network properties. The adjusted Rand Index (ARI) is a measure of similarities in clustering, quantifying the likelihood that pairs of microbial species would be assigned to the same cluster in both networks. An ARI value of 0 indicates random clustering across comparator groups, a value of 1 indicates identical clustering, and a value of -1 indicates perfect disagreement [41,42]. When comparing moderate to severe COVID-19, the ARI was 0.199 (p-value < 0.001), a modest but statistically significant finding indicating somewhat similar clustering of microbial species between networks. Jaccard's index (JI) evaluates differences among central nodes between our two severity-specific networks, where a value of 0 indicates completely different sets of central nodes and a value of 1 indicates identical central nodes [43]. While there were no statistically significant differences in overall centrality measures when comparing moderate to severe cases, there were alterations in proportion of positive edges network-wide (92.9% vs 100%, p-value < 0.001), indicating a loss of moderate negative correlations in severe COVID-19. For example, C. albicans, which was relatively more abundant in severe compared to moderate COVID-19, has 0 vs. 3 negative edges in each disease state, respectively, raising the possibility that the loss of negative selective pressure can promote the growth of certain microbial species in severe COVID-19. We identified 16 taxa as network hubs, i.e., species with high putative importance given their centrality to the surrounding microbial networks (Fig. 4, Additional file 1: Fig. S5, and Additional file 1: Table S5). Five species were identified as hubs in both moderate and severe disease (Blautia wexlerae, Eubacterium hallii, Gordonibacter pamelaeae, Odoribacter splanchnicus, and Alistipes shahii), while 11 were unique to one network or the other (Fig. 4, Additional file 1: Table S5 and Additional file 1: Fig. S5). Critically, 9 of these 16 identified hubs, including Blautia wexlerae and Eubacterium hallii, were shown to be differentially abundant by disease severity (Fisher's exact p-value = 0.03, Additional file 1: Table S2), and the relative abundance of two hubs, Eubacterium rectale and Alistipes putredenis, were associated with stool viral load. We further observed that highly connected clusters in moderate disease become fragmented in severe COVID-19, as evidenced by an increase in singletons (χ 2 p-value < 0.001). We also observed a decrease in the number of hub taxa and dynamic taxa-level cluster reassignment (Fig. 4). Notably, all but one of the hubs shown to be differentially abundant by disease severity belonged to the same cluster, suggesting that significant loss of these central taxa in severe disease may contribute to the observed network instability.

Predicted stool metabolites linked to disease severity
We next sought to evaluate whether changes in microbial communities affected capacity for local metabolite production. Using a validated computational workflow to generate putative metabolic profiles from stool metagenomes [33] (Methods), we found 57 of 80 well-predicted known stool metabolites to be differentially perturbed based on COVID-19 disease severity (all FDR-corrected p-value < 0.05; Fig. 5a and Additional file 1: Table S6). We identified disrupted bile acid metabolism in severe COVID-19, with relative enrichment of primary bile acids (chenodeoxycholate, cholate, and ketodeoxycholate) alongside depletion of secondary bile acids (lithocholate, lithocholic acid, and deoxycholic acid) (Fig. 5b). Similar to our microbial pathway analysis which revealed reductions in MetaCyc pathways related to SCFA production, predicted levels of butyrate, isobutyrate, and propionate were also reduced in severe COVID-19 (Additional file 1: Table S6). Furthermore, we confirmed prior data showing relative enrichment of bilirubin [44], creatine and polyamines (e.g., acetyl-spermidine [45]), and pantothenic acid [46] in severe COVID-19, as well as a relative depletion of deoxyinosine [46] (Additional file 1: Table S6).

Discussion
In a comparatively large US hospital-based cohort of diverse patients admitted with confirmed COVID-19 during the initial year of the pandemic, we found community-and species-level alterations linked to disease severity. Using a random forest machine learner, these microbial features could accurately classify patients based on disease severity, indicating that specific gut microbial configurations may be linked to a more severe disease course, a finding we validated in a separate independent cohort. Network analyses identified significant disruptions to gut ecologic topology in severe COVID-19. Differential abundance testing of microbial pathways and separate predicted stool metabolite-based analyses suggest that these disruptions may change the balance of bile acids and SCFAs in the gut, identifying novel treatment opportunities that may ameliorate the severity of COVID-19. We also found significant depletions of two microbes previously associated with long COVID, suggesting early gut microbial disturbances may precede the development of a long-term complication.
Determining who will require a higher level of care remains one of the most challenging questions facing clinicians caring for patients with COVID-19. Our machine learning algorithm demonstrated excellent discrimination between moderate and severe COVID-19 using only gut microbial features. Notably, the inclusion of clinical data did not significantly improve the classification accuracy of our model. Prior work has incorporated such information from the initial presentation [47], multicytokine panels [48], and previously validated illness severity scores [49] to forecast whether a given patient will suffer from a more severe COVID-19 course. However, based on their performance characteristics, these approaches appear to be less accurate than our microbiome-centered approach.
Prior investigations have observed similar community-and taxa-level alterations in microbial composition  Fig. 4 Comparative microbial assemblages in moderate vs. severe COVID-19. We assembled discrete microbial networks for moderate vs. severe disease to demonstrate significant ecological heterogeneity characterized by fractured clustering and taxa-level reassignment in severe disease. Species are represented by circles (nodes) and species-species correlations were weighted by strength of correlation (edges drawn if absolute Pearson's ⍴ > 0.4). Node size indicates normalized relative abundance, and node colors indicate cluster membership. Cluster colors are retained across networks if two or more taxa are shared. Edge color reflects the direction of correlation, with red edges indicating a negative, and green edges indicating a positive correlation, respectively. Node hubs have been labeled, while clusters are referred to by their nominate node, or the taxa with the highest edge count in a given cluster by network. Node positions are fixed between networks (full node maps and labels found in Fig S4 and Table S5  in COVID-19. In the earliest phase of the pandemic, a study from Hong Kong (n = 36) also demonstrated relative reductions in the group Eubacterium among the gut metagenomes of COVID-19-infected patients compared to referent populations, and like our work, found widespread depletion of typical gut colonizers such as Faecalibacterium and Roseburia spp. in severe COVID-19 [9]. In an expanded population of 100 patients, the same group reaffirmed a reduction in diversity and a loss of health-associated gut commensals in severe COVID-19 [13]. Finally, a study of 30 SARS-CoV-2 infected patients in mainland China using 16S rRNA-based sequencing similarly demonstrated a change in gut community structure with reductions in ɑ-diversity compared to referent counterparts [14]. Notably, they also achieved success in classifying stool samples from patients with COVID-19 compared to those from healthy controls or those infected with influenza, indicating the relatively distinct gut ecology of COVID-19. However, their classification tasks were conducted in a smaller population using supervised feature selection (i.e., the top results from their linear discriminant analysis) of genus-level taxa, and arguably, the role of a gut microbial biomarker in discriminating COVID-19 from non-infected individuals is uncertain now that SARS-CoV-2 testing is more widely available [50]. Our work offers insights beyond these broad characterizations of the gut microbiome in COVID-19. Among the eight taxa that were positively associated with stool SARS-CoV-2 viral load, several contribute to pro-inflammatory sulfur metabolism, such as Methanobrevibacter smithii and Bilophila wadsworthia [51][52][53]. Our finding of enriched R. hominis with increased stool viral load despite a corresponding decrease among patients with severe COVID-19 may suggest an interaction between stool SARS-CoV-2 viral load, R. hominis, and severe COVID-19. It is appreciated that gut microbial ecology influences the host immune response to viral respiratory infections [3][4][5][6]. Our identification of Blautia wexlerae and Eubacterium hallii as network hubs depleted in severe COVID-19 (both Lachnospiraceae implicated in other immune-mediated diseases [54]) suggests these bacteria may engage in important roles in the regulation of immunity to SARS-CoV-2. Predicted alteration of secondary bile acid metabolism in severe disease provides another mechanism by which changes in gut microbial communities may influence the immune response to SARS-CoV-2. Bile acids regulate mucosal and systemic immunity in several ways [55]. Prior work has suggested that secondary bile acids are the primary ligand for TGR5 [56] through which they may suppress pro-inflammatory signaling [55,57], resulting in impaired immunity to viral infections [58,59]. The predicted shift in bile acid pools may also result in increased regulation of bile acid-sensitive transcription factors, as increased primary bile acids will preferentially activate farsenoid X Receptor, while depletions in secondary bile acids will reduce activation of vitamin D receptor (VDR) [60,61] and pregnane X receptor (PXR) [62]. Decreased VDR/PXR signaling during active infection is associated with increased systemic inflammation and increased morbidity and mortality [63,64], possibly contributing to the clinical milieu observed in severe COVID-19. This is a particularly noteworthy hypothesis given emerging epidemiologic data on the link between diet [65], vitamin D status [66], and COVID-19 disease risk and severity, as well as early work linking depletion of secondary bile acids to COVID-19-related mortality [67].
Our study has several key strengths. First, we assembled a large representative cohort of patients at a U.S.-based tertiary care center for whom we collected relevant clinical metadata to complement serial stool sampling. Second, our computational workflow allowed us to not only link community-level changes in gut microbial ecology but species-resolved signatures of severe COVID-19, which we were able to validate in an external cohort of patients. Third, complementary MetaCyc pathway and predicted metabolite analyses further link these changes to alterations in bile acid pool and SCFA levels. Taken together, these observations serve as proof of principle that using NGS to interrogate gut microbial ecology may generate tractable hypotheses to be explored in follow-up investigations. Finally, our results fit well in the context of independent works from other groups-ending credence to our findings-and using a machine learning classifier, we demonstrate excellent accuracy in discriminating samples from moderate vs. severe COVID-19. These findings hint at the possibility that modulating gut microbial communities may be a viable disease prevention or therapeutic strategy in COVID-19.
We acknowledge several limitations. We were not positioned to assess whether findings differed on the basis of SARS-CoV-2 strain or variants. Our study enrolled patients from April 2020 to May 2021 during which genomic surveillance infrastructure in the USA was not equipped to comprehensively explore this question. Prior to the Delta variant wave beginning in June 2021, the majority of COVID-19 cases were either Alpha or other less consequential variants of interest [68]. As our study enrolled hospitalized patients with moderate COVID-19 to minimize differences between those hospitalized with severe disease, we are not positioned to explore what differences-if any-may exist between each group and their non-SARS-CoV-2 infected counterparts and the degree to which between-group differences are attributable solely to critical illness or prolonged hospitalization. However, we either adjusted for participant factors that differed between groups in our multivariable modeling (e.g., BMI, comorbid disease, hospital length of stay, and antiviral therapy) or were limited by the fact that distinguishing features such as advanced oxygen delivery and other ICU-level interventions were, by definition, established markers of severe COVID-19. Several sensitivity analyses including restricting our cohort to those unexposed to antibiotics, using just the first stool sample provided, or those provided prior to the median length of stay were each consistent with our main findings. Since most patients with severe disease were admitted to the ICU shortly after presentation, we were unable to prospectively collect a substantial number of pre-ICU samples in these patients, limiting our ability to classify or predict the development of severe COVID-19. Given the observational nature of our study, we cannot exclude the possibility of residual confounding. However, we adjusted for multiple potential confounders. All enrolled patients were hospitalized, which may minimize study heterogeneity at the expense of overall generalizability. We also assessed the gut microbiome at the earliest feasible time point on admission. This resulted in variation in the timing of collection, which limits our ability to infer causality. Absolute microbial abundance measurements could not be obtained. Finally, our collection protocol did not allow for the measurement of stool metabolites to validate our computational approach, without which it may be more accurate to consider these results as suggestive of altered capacity for metabolite class production rather than actual differences in quantifiable metabolite pools. Relatedly, diet and other unmeasured determinants of stool metabolite production are unlikely to be stable in a hospitalized population. Despite these limitations, our findings are intended to be hypothesis-generating to inform the continuum of research that may logically follow.

Conclusions
Leveraging the gut microbiome as a potential biomarker for disease severity and modulating this fragile ecology to improve COVID-19 outcomes each hold significant appeal in the fight to end this pandemic. Multidisciplinary approaches will be needed to confirm our early findings. Prospective validation of a non-invasive indicator predictive of disease severity could readily identify and target at-risk individuals for more aggressive therapy.