A high-risk gut microbiota configuration associates with fatal hyperinflammatory immune and metabolic responses to SARS-CoV-2

ABSTRACT Protection against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection and associated clinical sequelae requires well-coordinated metabolic and immune responses that limit viral spread and promote recovery of damaged systems. However, the role of the gut microbiota in regulating these responses has not been thoroughly investigated. In order to identify mechanisms underpinning microbiota interactions with host immune and metabolic systems that influence coronavirus disease 2019 (COVID-19) outcomes, we performed a multi-omics analysis on hospitalized COVID-19 patients and compared those with the most severe outcome (i.e. death, n = 41) to those with severe non-fatal disease (n = 89), or mild/moderate disease (n = 42), that recovered. A distinct subset of 8 cytokines (e.g. TSLP) and 140 metabolites (e.g. quinolinate) in sera identified those with a fatal outcome to infection. In addition, elevated levels of multiple pathobionts and lower levels of protective or anti-inflammatory microbes were observed in the fecal microbiome of those with the poorest clinical outcomes. Weighted gene correlation network analysis (WGCNA) identified modules that associated severity-associated cytokines with tryptophan metabolism, coagulation-linked fibrinopeptides, and bile acids with multiple pathobionts, such as Enterococcus. In contrast, less severe clinical outcomes are associated with clusters of anti-inflammatory microbes such as Bifidobacterium or Ruminococcus, short chain fatty acids (SCFAs) and IL-17A. Our study uncovered distinct mechanistic modules that link host and microbiome processes with fatal outcomes to SARS-CoV-2 infection. These features may be useful to identify at risk individuals, but also highlight a role for the microbiome in modifying hyperinflammatory responses to SARS-CoV-2 and other infectious agents.


Introduction
Infection with SARS-CoV-2 leads to a wide variety of potential outcomes from asymptomatic responses to acute respiratory distress and death. 1,2 While certain demographic factors such as age, male gender, and comorbidities that include obesity, cardiometabolic diseases and diabetes are associated with an increased risk for more severe disease, the molecular mechanisms that underpin disease pathophysiology remain poorly understood. Indeed, we still do not know if severe outcomes are due to direct effects of viral replication within target cells, to a dysregulated host immune response to the virus, to preexisting deficits in mechanisms of host resilience to infection, or to a combination of these factors. [3][4][5] Initially, SARS-CoV-2 infects angiotensinconverting enzyme 2 (ACE-2) expressing epithelial cells of the upper respiratory tract. If the infection remains limited to the upper respiratory tract, then this is usually associated with a mild disease course and rapid recovery. If the virus is not eliminated and infection persists, then other types of ACE-2 expressing cells can become infected. 6 In addition, viralinduced metabolic reprogramming and exaggerated immune responses generate a wide range of inflammatory mediators that disrupt organ homeostasis, impact host metabolism, drive a hypercoagulation state, impair epithelial barrier function and destroy host cells and tissues. [7][8][9][10][11] However, even among those who develop this cytokine storm, many can still make a full recovery, suggesting that additional factors may modulate host susceptibility to the most severe outcomes associated with COVID-19. One of these resilience factors might include the microbiome. [12][13][14][15] Human mucosal surfaces and body cavities harbor diverse communities of commensal microbes that play essential roles in the regulation of host metabolic responses, epithelial barrier function, immune education and immune regulation. [16][17][18][19][20] These effects are partially induced by activation of host pattern recognition receptors to microbialderived danger signals, but increasingly the role of bacterial metabolites in shaping host immune function is being recognized. [21][22][23] Immunoregulatory bacterial metabolites can trigger host G proteincoupled receptors (GPCRs), aryl hydrocarbon receptors (AhRs), nuclear hormone receptors such as the farnesoid X receptor, or can directly modulate gene expression through epigenetic mechanisms. Importantly, many immunoregulatory bacterial metabolites are derived from dietary substrates (e.g. fiber), linking diet and lifestyle to protection from infection via microbial mechanisms.
In this study, our primary aim was to identify the immune-metabolic-microbial interactions and biomarkers that predict the most severe outcomes to SARS-CoV-2 infection in a well-characterized cohort of patients hospitalized with COVID-19. In addition, we wished to identify clusters of patient metadata features that might provide novel mechanistic insights into the disease pathophysiology. Lastly, we wished to extend our understanding of the molecular processes within the holobiont that mediate resilience to severe biological challenges, such as viral infection.

Systemic levels of immune mediators correlate with disease severity
While changes in circulating cytokine levels due to SARS-CoV-2 infection are already well described, the immune mediators that distinguish survivors from non-survivors in severely ill patients have not been clearly identified. To better understand the immune processes that might distinguish these patients, we measured the levels of 54 immune mediators in the earliest serum sample obtained following study enrollment after admission to the intensive care unit (ICU; severe  or the hospital ward (mild to moderate COVID-19) from 172 hospitalized patients with PCR-confirmed SARS-CoV-2 causing COVID-19. Patient demographic details are shown in Table 1. Those with mild/moderate COVID-19 (n = 42) were younger, more likely to be female, less frequently obese, required fewer medications and had fewer comorbidities compared to those with severe COVID-19 (n = 130). However, there were no differences in demographics, medication use or comorbidities in those severely ill patients that survived infection (n = 89), compared to those COVID-19 patients with a fatal outcome (n = 41). In contrast, principal component analysis of serum immune mediators demonstrated a clear separation between patients with different COVID-19 disease outcomes (Figure 1(a)). Compared to healthy volunteers (n = 29), levels of 36 circulating immune mediators were significantly differed (30 higher and 6 lower) in those hospitalized with COVID-19 ( Figure 1(b) and Figure S1). Of these mediators, levels of 28 were significantly different between patients with mild/moderate COVID-19 compared to patients with severe disease (Figure 1(b)). Within the severely ill group, the levels of eight circulating immune mediators (soluble intercellular adhesion molecule-1 (sICAM-1), monocyte chemoattractant protein-1 (MCP-1), interleukin (IL)-8, macrophage-derived chemokine (MDC), interferon gamma-induced protein-10 (IP-10), IL-15, IL-1 receptor antagonist (RA) and thymic stromal lymphopoietin (TSLP)) were significantly different between those that survived and those that died (Figure 1(b) and (c)).

Systemic metabolic responses associated with disease severity
In addition to measuring serum cytokines, we quantified and compared metabolite levels in the first serum sample obtained following study recruitment after admission to the ICU or hospital ward for COVID-19 patients with mild/moderate disease (n = 25), COVID-19 patients with severe disease that survived (n = 75) or COVID-19 patients with severe disease that succumbed to death (n = 39). Distinct differences in circulating metabolites were evident between each of the groups (Figure 2(a) and (b)). Metabolic processes were dramatically different in patients during acute SARS-CoV-2 infection, whereby levels of 377 metabolites were significantly different (adjusted p < .05) between healthy volunteers (n = 20) and those with mild/moderate COVID-19 ( Figure 2(b)). These differences were further exaggerated in COVID-19 patients with severe disease (583 metabolites, adjusted p < .05), in particular those with a fatal outcome (659 metabolites, adjusted p < .05), when compared to healthy volunteers. Within the severely ill patients, 140 metabolites distinguished those that survived versus those that died. The metabolites that contribute most to the differences between the groups included those involved in tryptophan metabolism, polyamine metabolism, histidine metabolism, lipid metabolism, bile acid metabolism and antioxidant responses such as the plasmalogens (Figure 2(c) and Figure S2). Random forest analysis suggested a good discriminatory power for distinguishing COVID-19 disease severity or fatality based solely on a selection of circulating metabolites (Figure 2c and Figure S3), underlining the robustness of these differences. Given the substantial and significant differences in metabolite levels, we examined in more detail the most significantly impacted pathways associated with COVID-19 severity (Figure 3(a)). Interestingly, levels of sulfonated bile acids were particularly disrupted with disease severity. Host tryptophan metabolism was associated with a heavy depletion of tryptophan, with enhanced generation of kynurenate, kynurenine and quinolinate, at the expense of serotonin synthesis in COVID-19 patients (Figure 3(a) and Figure S4(a)). In contrast, microbial tryptophan metabolites were present at lower levels in the serum of those with the worst outcome ( Figure 3(a) and Figure S4(b)). Changes in circulating microbial metabolites may be due in part to an impaired gut barrier (as indicated by increased serum SCFA levels and lower citrulline levels, Figure  S4(c) and (d)), or may reflect changes in the composition or metabolism of the gut microbiome. Overall, metabolites associated with microbial metabolism (as described by Bar et al.) 24 were significantly altered in those with severe disease and those with a fatal outcome ( Figure S4(e)).
Next, we performed a weighted co-expression network analysis restricted to the COVID-19 patients, to identify communities of co-abundant metabolites. Positive correlations between metabolites (Spearman, adjusted p < .0005) were used to build the network. The analysis identified six communities (c1-c6) of highly intercorrelated metabolites based on the Leiden algorithm [two iterations, ModularityVertexPartition, weighted network ( Figure 3(b)]. Primary and secondary bile acid metabolism are contained in c3, SCFA in c5, while tryptophan and histidine metabolism are in c1, c2 and c5 (Figure 3(c)). The central community (c1) with the most interconnected metabolites, central metabolites, and greatest influence on the global dynamics of the network includes mannose (Figure 3(d)), which is a known inflammatory biomarker and reported to be associated with COVID-19 severity. 25 Furthermore, the metabolites that are significantly different between COVID-19 severe patients with or without a fatal outcome are primarily found within community c1 (Figure 3(e)).

Differences in the gut microbiome associate with disease severity and death
To investigate the possible involvement of the gut microbiome in these immune and metabolic changes, we profiled the microbiome by sequencing 16S rRNA gene amplicons from the first fecal samples collected following study recruitment after admission to the ICU or hospital ward for COVID-19. From the 99 hospitalized COVID-19 patients with available stool samples for 16S amplicon sequencing, 32 had mild/moderate disease, 45 had severe disease and survived, while the remaining 22 patients had severe disease with a fatal outcome. Global measures of microbiome alpha diversity were not different between clinical groups, with no significant difference detected in Shannon indices as well as in the number of detected taxa at the level of Operational Taxonomic Units (OTUs), species or genus levels between the three disease outcome groups ( Figure S5). However, Envfit-based analysis of the Principal Coordinates revealed a significant difference in gut microbiome composition (beta diversity) between the three COVID-19 disease severity groups, irrespective of the distance measures used (Figure 4(a) and Figure S6). We next investigated these differences in microbiome profiles in an unsupervised manner, i.e. without utilizing the disease outcome information. Using an iterative enterotyping-based approach applied on the Principal coordinates (See Methods), the microbiomes could be optimally clustered into two configurations (MicrobiomeGroup1 and MicrobiomeGroup2), resolved clearly along the first Principal Coordinate (Figure 4(b) and (c)). Notably, there were significant differences in the proportions of the two distinct microbiome configurations in the clinical outcome groups (Chi-square test estimate = 11.23, p-value = 0.0036, Figure 4(d)). MicrobiomeGroup1 was over-represented in severe COVID-19 patients with a fatal outcome, while Microbiome Group2 was associated with those with mild/moderate symptoms. Strikingly, within the severe outcome group, individuals who were classified into the high-risk MicrobiomeGroup1 had significantly higher levels of cytokines associated with both fatality and severity (P = .02; Mann-Whitney Test), with higher (albeit not statistically significant) levels of cytokines associated only with disease severity (P = .12, Mann-Whitney Test) ( Figure S7(a)). A similar trend was observed for the association between fatality and severity associated cytokines and the microbiome configuration in COVID-19 patients with mild/ moderate symptoms, but this difference did not reach statistical significance ( Figure S7(b)). In addition, severity associated metabolite levels was higher in those patients with severe disease in the MicrobiomeGroup1 ( Figure S8).
We next investigated the genus-level composition differences across the two microbiome configurations by performing ordinary-least square (OLS)-based regression analysis to measure the association between abundance of microbial genera and the PCo1 axis values after adjusting for confounders (age, gender, BMI and hospital location). A number of genera showed significant associations with PCo1 with FDR ≤ 0.15 (Benjamini-Hochberg corrected), even after confounder adjustment. While genus-level groups such as Enterococcus and Lachnoclostridium were associated negatively with PCo1 (high relative abundance in the high risk MicrobiomeGroup1), others including Roseburia, Dorea, Faecalibacterium, Bifidobacterium, Fusicatenibacter and multiple Ruminococcus species showed the opposite trend ( Figure 4(e)). Relaxing the thresholds identified additional genera that showed nominally significant association with PCo1 (P ≤ .05). The high risk MicrobiomeGroup1 was characterized by higher levels of multiple pathobionts, as operationally defined in our previous work including Enterococcus, Eggerthella, Lachnoclostridium, Erysipelatoclostridium, Streptococcus, Flavonifractor and lower levels of multiple taxa known to be associated with antiinflammatory or protective immune responses (including Faecalibacterium, Agathobacter, Dorea, Coprococcus, Lachnospiraceae, Bifidobacterium) ( Figure 4(e)). 26,27 The relative abundance data for the genus-level microbiome features and the statistical comparisons between the two microbiome groups are included in supplementary Table S1 and Table S2, respectively. Many of the observed differences in the microbiome were significantly associated with changes in levels of circulating immune mediators ( Figure S9).
We did not perform 16S amplicon sequencing on the healthy controls in this study as we did not have fecal samples from them. However, we have compared the COVID-19 patient 16S microbiome composition data to that of the healthy control 16S sequence data from two cohorts recently published from our group. 28,29 Merging these two control cohorts gave us 531 genus-level microbiome profiles of fecal samples of control individuals of European ancestry with a similar age range to the COVID-19 patients in our current study. Using the median kendall distances of the genus-level profiles, the high-risk MicrobiomeGroup 1 had significantly higher distances from the healthy reference controls, compared to the low-risk MicrobiomeGroup 2 ( Figure S10(a)). This can also be observed across the PCo1 landscape ( Figure S10(b)). There were seven genera whose loss was associated with increase in distance from the healthy reference ( Figure S10 (c)). These were Anaerostipes, Ruminococcus_2, Faecalibacterium, Agathobacter, Ruminococcus_1, Bifidobacterium and Collinsella. The genera that are enriched with increase in distance from the healthy reference controls are Enterococcus, Coprobacter and Streptococcus.

Immune-metabolite-microbiome modules correlate with COVID-19 disease outcomes
Correlation network analysis is a powerful tool for revealing associations of diverse features within patient datasets. Feature-association networks were computed using the Weighted gene correlation network analysis (WGCNA) approach (see Methods) performed on 1,469 features (54 cytokines, 1,146 metabolites and 269 microbial genera) using signed Spearman correlations with a soft-power threshold of 7 ( Figure S11(a)) from the 70 hospitalized COVID-19 patients with complete data for all three data layers. A total of 14 modules (annotated as different colors) were identified, 5 of which had a significant association with disease outcome (Benjamini-Hochberg FDR ≤ 0.05) and 2 modules showed nominal associations (P ≤ .05 and FDR ≤ 0.1) ( Figure S11(b) and (c)). The module (annotated as 'turquoise') that showed significant positive association with disease severity and death contained most of the severity associated cytokines (as identified in Figure 1), metabolites ( Figure S3) and microbial genera identified above (Figure 4(e)), combined with kynurenine associated metabolism products and coagulation linked fibrinopeptides ( Figure 5 (a)). Two modules (annotated as 'brown' and 'tan') were nominally positively associated with a poor outcome. Of these, the brown module contained a triad of pathobionts linked to urobilinogen ( Figure S12(a)), while the tan module was enriched for sulfonated bile acids ( Figure S12(b)). In contrast, four modules, annotated as 'red', 'blue', 'black' and 'yellow', were significantly negatively associated with COVID-19 severity and death. The first module (red) contained the anti-inflammatory Ruminococcus_2 clade, linked with tryptophan, alanine and the SCFAs butyrate/isobutyrate and valerate ( Figure 5(b) and Figure S13). The second module (blue) that negatively associated with disease severity contains a cluster of beneficial microbial taxa (including Bifidobacterium), bilirubin degradation products, TARC and IL-17A ( Figure S14). The third module (black) exclusively contains metabolites, in particular fatty acid derivatives ( Figure S15), while the final significant module (yellow) contains Roseburia, Fusicatenibacter, Romboutsia linked with sphingomyelin and carnitine-derived products ( Figure S16).

Discussion
Despite the substantial literature published on SARS-CoV-2, the molecular mechanisms underpinning positive versus negative clinical outcomes remain poorly defined. In this study, we examined the differences in circulating inflammatory markers and metabolites in sera, and the composition of the gut microbiota, in a large group of hospitalized patients with COVID-19. We have identified several potential regulatory nodes whereby integrated immune, metabolic and microbiome processes contribute to susceptibility or resilience to SARS-CoV-2 infection associated damage.
Our identification of circulating inflammatory mediators that associate with COVID-19 disease severity such as CRP and IL-6 are consistent with previous reports and support the hypothesis that an overly aggressive immune response contributes to immunopathology and severity. 30,31 In addition to severity associated factors, we have identified a subset of eight cytokines that are further dysregulated in severe patients with a fatal outcome. Higher levels of IP-10 and IL-15 indicate greater activation of a T helper 1 (Th1)-associated innate anti-viral response, while a significant reduction in MDC levels may reflect the inhibitory effect of a Th1 environment on Th2 cytokines such as MDC. We were particularly interested in TSLP as this cytokine is an epithelial cell-derived alarmin, which is released by injured stromal cells to recruit and activate innate immune cells, and its blockade is currently being investigated in asthma clinical studies. [32][33][34] In combination with the chemokines MCP-1 and IL-8, and sICAM-1 (which modulates leukocyte adhesion and migration across endothelial cells), elevated TSLP levels indicate a greater amount of epithelial tissue damage and inflammatory cell recruitment to the damaged sites in patients who do not recover from SARS-CoV-2 infection. 35 As SARS-CoV-2 is a lytic virus, it is possible that viral replication in epithelial cells may directly drive TSLP levels in sera, although indirect effects on epithelial cells within the respiratory tract or gut might also induce TSLP release. Importantly, TSLP levels were previously shown to be elevated in patients with long COVID, suggesting that long-term impacts of SARS-CoV-2 on epithelial cells should be examined in more detail, potentially guiding future therapeutic interventions. 36 While elevated proinflammatory cytokine levels are a common feature also associated with severe responses to other respiratory viruses such as influenza A (H1N1), the pattern of polyfunctional cytokine responses described in the current study seems to be specifically associated with severe COVID-19, potentially due to the broad infective capacity of SARS-CoV-2 to invade tissues and organs outside the respiratory tract. 37,38 Significant metabolic reprogramming and compensatory responses are evident in COVID-19 patients with severe disease and particularly in those with a fatal outcome. Decreased serum levels of plasmalogens suggest a significant level of systemic oxidative stress as these sacrificial phospholipids are preferentially oxidized to protect more vulnerable membrane lipids such as polyunsaturated fatty acids. 39 Altered tryptophan metabolism was particularly interesting to observe as the profound shutdown in serotonin production coupled with accumulation of quinolinic acid indicated a shift from production of neuroprotective compounds to production of neurotoxic compounds, which might be clinically relevant. 40 An imbalance between host and microbial tryptophan metabolism was also evident as serum kynurenine levels increased, while products of bacterial tryptophan metabolism such as indoleacetic acid were significantly decreased in those with severe and fatal disease. 41 These are important AhR ligands that can contribute to immune regulatory responses, can drive an "exhaustion" phenotype in immune effector cells, and are important for maintenance of the gut epithelial barrier by induction of 43 Changes in host tryptophan metabolism have also been associated with a poor outcome to influenza infection, but changes in microbiota metabolism of tryptophan have not been described. 44 Other significantly different metabolites such as the polyamines putrescine and spermidine play important roles in protecting against inflammatory responses within the airways. 45 In addition, changes in secondary bile acid serum levels indicate significant disruption of microbial metabolism and/or changes in the gut barrier. Secondary bile acids significantly impact regulatory and effector immune responses, which may be relevant for the development of severe COVID-19. 46,47 Increased levels of sulfonated bile acids in serum also indicate significant disruption of bile acid metabolism in severely ill COVID-19 patients as sulfonation is an important detoxification mechanism that prevents reabsorption of bile acids from the gut and promotes their elimination in feces. 48 We identified a high-risk gut microbiome configuration associated with an inflamed host phenotype and increased risk of the worst disease outcomes.
Several pathobionts including Enterococcus were enriched in severe disease, while well described immune regulatory microbes such as Bifidobacterium and Ruminococcus were enriched in those who survived. 49,50 Similar microbiome configurations have been described in other settings such as increasing age, whereby a decrease of the core protective microbiome accompanied by an increase of pathobionts was observed. 51 In addition, acquisition of this subset of disease-associated taxa has been shown to shift the metabolic state to a disease-like state. 27 These changes in the microbiome may have happened gradually over time and could potentially make the host less resilient to SARS-CoV-2 infection. It is possible that this highrisk gut microbiome configuration may not only associate with a poor outcome to SARS-CoV-2 infection but might also be important for appropriate responses to other respiratory infections, such as influenza.
The hyper-inflammatory state observed in COVID-19 patients with a fatal outcome implies a failure in the negative feedback mechanisms that should restrain the devastating overproduction of inflammatory cytokines and soluble mediators, which lead to multiorgan failure. Our integrated analysis of microbiome features, cytokines and metabolites suggests that important microbial-derived immunoregulatory processes that contribute to negative feedback mechanisms may be lacking in those with the most severe outcomes to SARS-CoV -2 infection. Alternatively, increased levels of proinflammatory pathobionts may drive excessive proinflammatory responses that cannot be contained by the regular feedback mechanisms. While further studies will be required to determine causal interactions, this study supports the hypothesis that successful responses to infectious agents such as SARS-CoV-2 involve the gut microbiome mediated by effects on metabolism and host inflammatory processes. Microbiome profiling may assist in the early identification of patients at high risk of severe symptoms, while targeting the microbiome via appropriately selected probiotics and/or prebiotics may enable the immune system to respond to infectious challenges in a robust, effective and well controlled manner.

Study cohort
We performed an investigator-initiated, prospective multicentre cohort study of adult (≥18 years) patients who were admitted with Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV -2) to four different hospitals in Switzerland and Ireland. Infection was confirmed by SARS-CoV-2 polymerase-chain reaction (PCR) from an upper or lower respiratory specimen. Exclusion criteria included COVID-19 diagnosis after discharge from the ICU. Recruitment started in August 2020 and in total we recruited 172 hospitalized patients from St. Gallen, Switzerland (n = 37), Geneva, Switzerland (n = 50), Ticino, Switzerland (n = 77) and Cork, Ireland (n = 8). All patients or patient representatives signed a patient informed consent. The study was approved by local ethics committees (EKOS 20/058 for the three Swiss sites and The Clinical Research Ethics Committee of the Cork Teaching Hospitals for Cork University Hospital). Patients were enrolled typically within 24-48 h after admission to the intensive care unit (ICU) or a hospital ward. Baseline characteristics, underlying comorbidities and medication use at the time of sampling were collected and are summarized in Table 1. All medical procedures and treatments were left at the discretion of the treating physicians but documented in the database such as complications during ICU stay and outcomes until hospital discharge. Patients were categorized to have mild disease when there were no radiographic indications of pneumonia and moderate disease if pneumonia with fever and respiratory tract symptoms were present. Severe disease was defined as a respiratory rate ≥30 breaths per minute, oxygen saturation ≤93% when breathing ambient air or PaO2/FiO2 ≤ 300 mm Hg, or anyone that required mechanical ventilation. Only those that died during their hospital stay were recorded as a SARS-CoV -2-related death in this study. Serum and fecal samples were collected as soon as possible following enrollment into the study and immediately stored frozen at −80°C at the clinical site. Serum from healthy volunteers was obtained prior to the start of the COVID-19 pandemic and healthy volunteers were prescreened for any acute or chronic immune, metabolic or infectious disorder that would exclude them from being used as a healthy control.

Cytokine analysis
We examined the levels of 54 cytokines and growth factors (using MSD multiplex kits according to manufacturer's instructions) in the serum of 172 hospitalized COVID-19 patients. Serum from patients was typically obtained within 24 hours after study enrollment. Sera obtained prior to the pandemic from 29 healthy volunteers were ana-

Metabolomics
Untargeted metabolomics on patient sera was performed by Metabolon TM using the HD4 platform. Briefly, all methods utilized a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization (HESI-II) source and Orbitrap mass analyzer operated at 35,000 mass resolution. The sample extract was dried and then reconstituted in solvents compatible to each of the four methods. One aliquot was analyzed using acidic positive ion conditions, chromatographically optimized for more hydrophilic compounds. Another aliquot was also analyzed using acidic positive ion conditions; however, it was chromatographically optimized for more hydrophobic compounds. Another aliquot was analyzed using basic negative ion optimized conditions using a separate dedicated C18 column. The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1 × 150 mm, 1.7 µm) using a gradient consisting of water and acetonitrile with 10 mM Ammonium Formate, pH 10.8. The MS analysis alternated between MS and data-dependent MSn scans using dynamic exclusion. The scan range varied slightly between the methods but covered 70-1000 m/z.

16S sequencing
Fecal samples were obtained as soon as possible following hospitalization. Total community DNA was extracted from fecal samples by a combined Repeat Bead Beating -Qiagen DNA extraction method, and the V3 dash V4 region of the 16S gene was amplified and sequenced as previously described. 52 The uniquely barcoded amplicons were sequenced on an Illumina MiSeq platform (Illumina, California, USA) utilizing 2 × 300 bp chemistry.

Bioinformatic analysis
From the Log2 transformed metabolomics data obtained from Metabolon, any metabolite with no variance among samples was removed. Pairwise differential abundance analysis was performed between conditions using R package LIMMA.
Benjamini-Hochberg correction (BH) was applied for each comparison. R packages Boruta was applied for feature and tree number selection before random forest analysis. Random forest classifiers were built with the most important features, 1000 trees, mtry of 1 and 10-fold cross-validation using R packages caret and randomForest. They were evaluated using confusion matrices/roc curves. For association analysis, significant positive correlations (Spearman, FDR <0.0005) were extracted and used to build the network using python igraph (https://igraph.org/python/). The strength of the connections and relevance of the network were evaluated by plotting distribution of correlation coefficients and comparison of the network to a random network with similar dimensions. Community detection was performed using the Leiden algorithm from the python module leidenalg (https://leidenalg.readthedocs.io/en/stable/ index.html). For each community large enough (N > 30), metabolite set enrichment analysis (MSEA) was performed. For metabolite set enrichment analysis (MSEA), all MetabolonTM terms were extracted with their corresponding metabolites as reference. Python 3 gseapy package was used to perform a hypergeometric test between list of significant metabolites and reference. Importance plots, dot plots, bar plots, pca plots were produced with R package ggplot2. Heatmaps were designed with the R package ComplexeHeatmap. Networks were represented using Cytoscape 3.6.1 and metabolites of interest highlighted.
For the microbiome analysis, the raw Illumina reads obtained for each sample were quality-filtered using the trimmomatic program, using the default parameters. 53 The quality filtered reads were then taxonomically classified using both DADA2 54 (for read-level genus classification and identification of amplicon sequence variants or ASVs within each sample) and Spingo 55 (for species level classification). Amplicon Sequence Variants obtained using DADA2 for all the samples were then further merged by performing into Operational Taxonomic Units (OTUs) using the denovo-sequence-based clustering using the qiime package. 56 Subsequent downstream analyses of the taxonomic profiles (at all three levels, namely genus, species and OTU) as well as integrated analysis of taxonomic profiles with cytokine profiles and the metabolome were performed using various modules/packages of the R programming interface (v 4.0.3; R Core Team 2020). Estimates of alpha diversity were computed using the diversity function of the vegan package of R. Principal Coordinate Analyses (PCoA) were computed using the ade4 package. The envfit function of the vegan package was used to perform the envfitbased analysis using the three top Principal Coordinates. The analysis fits environmental or host metadata factors (in this case the Covid outcome) onto the ordination space (in this case, the PCoA) and attempts to identify significant associations. Enterotyping of the gut microbiome profiles was performed as described in a previous study from our group. 57 Two group comparison of microbiome abundances were performed using the Mann-Whitney tests (using the wilcox.test function of R stats package). For more than two-group comparisons, pairwise comparisons within groups were computed using Mann-Whitney tests. The p-values were corrected using Benjamini-Hochberg FDR correction (p.adjust function of the stats package). For associating microbiome features with variations in Principal Coordinates (PCo) (specifically the Principal Coordinate 1), ordinary Least Squares (OLS) regressions were performed after adjusting for confounders (age, gender, BMI and hospital location) using the glm function of the stats package (as has been performed in previous studies). 27,58 Specifically, for each feature, in this case the Microbiome Genus, two models were computed: Model 1: glm(Microbiome Genus ~ Hospital Location (Randomly assigned numeric codes from 1 to 4) + Age + Gender (encoded as 1 or 2) + BMI) Model 2: glm(Microbiome Genus ~ Hospital Location (Randomly assigned numeric codes from 1 to 4) + Age + Gender (encoded as 1 or 2) + BMI + Principal Coordinate) The adjusted R-squares and AIC values were then obtained to judge the goodness-of-fit of both the models. The directionality of the association of the Principal Coordinate with a Microbiome genus in Model 2 was assigned based on the sign of the corresponding coefficient in Model 2. The significance of improvement (P-values) of Model 2 with respect to Model 1 was judged using log-likelihood tests. All P-values were subsequently adjusted (across features) using the Benjamini-Hochberg approach to compute the FDR or Q-values.
Correlation analysis of associations among features in three data layers (genus-level microbiome, metabolome and cytokine profiles) was performed using the Weighted Gene-Coexpression Network Analysis (WGCNA). 59 While originally devised for computing gene coexpression networks, WGCNA is now being used in studies to integrate data from multiple OMICs layers. 60,61 In this study, the WGCNA was performed using an optimal soft-power threshold of 7 for scale-free topology. Using hierarchical clustering and topology overlap measures (TOM), we identified that the features from the three data layers could be optimally grouped into 14 modules, which were then investigated for association with disease symptoms using OLS models. The association networks within each module were then computed using the ReBoot approach as implemented in the ccrepe workflow. 62