Longitudinal metabolomics of human plasma reveals prognostic markers of COVID-19 disease severity.

There is an urgent need to identify which COVID-19 patients will develop life-threatening illness so that medical resources can be optimally allocated and rapid treatment can be administered early in the disease course, when clinical management is most effective. To aid in the prognostic classification of disease severity, we perform untargeted metabolomics on plasma from 339 patients, with samples collected at six longitudinal time points. Using the temporal metabolic profiles and machine learning, we build a predictive model of disease severity. We discover that a panel of metabolites measured at the time of study entry successfully determines disease severity. Through analysis of longitudinal samples, we confirm that most of these markers are directly related to disease progression and that their levels return to baseline upon disease recovery. Finally, we validate that these metabolites are also altered in a hamster model of COVID-19.

Introduction

Coronavirus disease 2019 (COVID-19), which is caused by infection with the novel coronavirus SARS-CoV-2, has led to a global health crisis. To date, more than 150 million cases of COVID-19 have been reported worldwide and resulted in more than 3.2 million deaths. The infection-fatality rate of SARS-CoV-2 can be reduced with the appropriate care (e.g., intensive care unit beds, oxygen, staff, extracorporeal-life support, and therapeutics). Unfortunately, hospital resources can quickly become depleted in situations when cases spike. Although vaccination efforts are underway worldwide, SARS-CoV-2 infections continue to increase rapidly in a number of countries, such as India and Brazil, where medical facilities are being overwhelmed.

Patients who develop critical illness from COVID-19 are best treated early in the disease course before the onset of severe symptoms., It is currently difficult to determine which subset of patients will develop life-threatening disease and, therefore, most benefit from receiving treatment when resources are limited. Early identification of these patients would allow optimal allocation of care. To that end, the objective of the current study was to identify metabolites in patient plasma that accurately predict life-threatening cases of COVID-19 before the onset of severe symptoms.

SARS-CoV-2 is an enveloped, single-stranded positive-sense RNA virus that gains entry into host cells through binding of the viral S protein to the angiotensin-converting enzyme 2 (ACE2) receptor., Multiple studies have established that patients infected with SARS-CoV-2 have metabolic dysregulation, possibly because of immune-triggered inflammation or other changes in host physiology.10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 To date, however, unique alterations in metabolites upon SARS-CoV-2 infection have not been validated in large patient cohorts that have been profiled longitudinally from early after infection through recovery. Longitudinal assessment of metabolic trends is necessary to assess which changes are indicative of disease course.

In this study, we performed untargeted metabolomics profiling on the polar and non-polar fractions of 700 human plasma samples collected from 339 patients as part of the WU-350 cohort recruited during the first phase of the pandemic in St. Louis, Missouri. Untargeted metabolomics allows for the unbiased profiling of the human metabolome and has been successful at discovering metabolite biomarkers associated with disease pathology. Using machine learning (ML), we built a model of COVID-19 disease severity based on the metabolic profiles of samples collected from patients at study enrollment. The model led us to identify a panel of unique metabolite biomarkers that were highly indicative of disease severity. We confirmed that most of these metabolites were directly related to SARS-CoV-2 infections through comparison with patient demographics, comorbidities, clinical measurements, and longitudinal samples taken from individuals over the course of disease progression. Lastly, we validated that the same biomarkers appeared in an established hamster model of SARS-CoV-2 infection.23, 24, 25

Results

Clinical cohort WU-350: demographics

The clinical cohort presented in this study consisted of 155 female and 184 male participants. Of the 339 patients, 272 were considered SARS-CoV-2 positive (COV+) whereas 67 were considered SARS-CoV-2 negative (COV-). The COV+ and COV- assignment of each individual was confirmed by nasopharyngeal swab polymerase chain reaction (PCR). Significant differences were observed in several demographic factors for the COV+ cohort compared with the COV- cohort (Table 1). The COV+ group contained significantly older study participants (p < 0.0001; Figure S1A). The COV+ group also had significantly more African American, male, and non-smoking individuals. There was no significant difference in the body mass index (BMI) between the two groups (Figure S1B).

Table 1

Demographics of all subjects

Parameter	COV-	COV+	p value
N	67	272
Gender (M/F)	28/39	156/116	0.0220
Age (years)	48 ± 16	60 ± 17	<0.0001
Race (African American/white/other)	32/35/0	202/66/4	<0.0001
Body mass index	30 ± 8	31 ± 9	0.3987
COVID-19-like symptomsa
Any number of COVID19-like symptoms	65	242	0.0437
Fever	29	128	0.5788
Chills	13	46	0.6300
Conjunctival congestion	1	1	0.2815
Nasal congestion	7	10	0.0229
Headache	18	23	<0.0001
Cough	35	147	0.7907
Sore throat	14	23	0.0034
Shortness of breath	44	168	0.5540
Nausea or vomiting	15	41	0.1487
Diarrhea	11	38	0.6098
Myalgia	16	66	0.9476
Fatigue	21	53	0.0353
Loss of taste or smellb	0	4	0.3180
Asymptomaticc	NA	21
Acute respiratory failured	9	100	0.0002
Acute renal failured	2	68	<0.0001
Comorbidities
Chronic kidney diseasee	4	56	0.0050
Diabetese	15	120	0.0011
Cancere	8	24	0.4345
COVID-19 drug treatmentf
Remdesivir	2	43	0.0056
Dexamethasone	13	63	0.5087
High/Low arterial pH	5	83	0.0001
Current smoker	18 (27%)	36 (13.2%)	0.0063
Hospital admissiong	26 (38.8%)	252 (92.6%)	<0.0001
ICU admissiong	10 (14.9%)	129 (47.4%)	<0.0001
Intubation and mechanical ventilation	4 (6.0%)	70 (25.7%)	0.0005
Deceased	6 (8.9%)	65 (23.9%)	0.0071
Deceased because of COVID-19h
30 day mortality	N/A	39 (76.5%)
60 day mortality	N/A	49 (96.1%)
90 day mortality	N/A	50 (98.0%)
overall mortality	N/A	51 (100%)

Table includes both training and test cohorts. Percentages are shown as the percentage of the group (COV- or COV+). Demographic data were last updated May 20, 2021. Data are presented as means ± SD, p values of numeric parameters were calculated by using a two-tailed Student’s t test with unequal variance, p values of categorical parameters were calculated by using a chi-square test. Abbreviations: M, male; F, female; N/A, not applicable.

Two COV- individuals were without symptoms but were exposed to a SARS-CoV-2-positive individual. Nine COV+ individuals had other symptoms (e.g., confusion, lethargy, an altered mental state, or breathing anomalies). Clinical metadata last updated October 16, 2020.

CDC guideline symptom was added to the symptom questionnaire late in the study; parameter is not available for most of the subjects.

A SARS-CoV-2 test was routinely administered at presentation to the hospital for reasons other than COVID-19 (e.g., accidents, pre-operation tests, regular checkups, cancer screening, injuries, or exposure to a SARS-CoV-2-positive individual).

During present hospitalization.

Recorded up to 1 year before the current admission or up to 1 year before the d0 sample for those who were not admitted to the hospital.

At any point during hospital stay.

Hospital and/or ICU admission of COV− group was for reasons other than COVID-19 (e.g., accidents, acute respiratory failure from bacterial pneumonia, intentional self-harm, possible heart failure, hypertension, trauma, and cancer).

Percentages shown as the percentage of total number of patients who died of COVID-19.

Of the 272 COV+ individuals, 252 were admitted to the hospital and 129 of those patients were admitted to the intensive care unit (ICU). As expected, the incidence of both factors (hospitalization and ICU admission) was significantly increased in the COV+ cohort. Treatment of severe COVID-19 cases often results in intubation and mechanical ventilation. In total, 70 of the COV+ patients required mechanical ventilation, whereas only four COV- individuals required mechanical ventilation. The therapeutics administered to treat COVID-19 in this cohort were dexamethasone and remdesivir. In total, 43 COV+ patients were treated with remdesivir and 63 with dexamethasone during their hospital stay (Table 1). The mortality rate in the COV+ group was 23.9%, which was significantly higher than the 8.9% mortality rate in the COV- group. A total of 65 COV- patients died, with 51 of the deaths being attributed to COVID-19 and 14 being attributed to other causes.

Of the 272 COV+ patients, 242 showed at least one COVID-19-related symptom mentioned by the Centers for Disease Control and Prevention (CDC), including fever, chills, conjunctival congestion, nasal congestion, headaches, cough, sore throat, shortness of breath, nausea or vomiting, diarrhea, myalgia, fatigue, and loss of taste or smell. Of the remaining COV+ cases, 21 showed no symptoms and were classified as asymptomatic. Nine subjects in the COV+ group displayed other symptoms, such as confusion, lethargy, an altered mental state, or breathing anomalies. Of the 67 COV- cases, 65 presented with at least one COVID-19-related symptom upon study entry, whereas two received a SARS-CoV-2 PCR test upon exposure to a SARS-CoV-2-positive individual. The frequency of COVID-19-related symptoms is shown in Table 1, and the distributions across the COV+ and COV- cohorts are depicted in Figure S1C. In both the COV- and COV+ groups, the number of COVID-19-related symptoms reported per individual was comparable. The breakdown of how many symptoms were experienced per individual in both the COV+ and COV- groups is shown in Figures S1D and S1E.

Next, we examined the distribution of comorbidities and self-reported medical history presented in the WU-350 cohort. Medical conditions were assigned based on ICD10 codes and the Elixhauser naming conventions. Chronic kidney disease (CKD) and diabetes (recorded up to 1 year before the current admission and up to 1 year before the day 0 (d0) sample for those who were not hospitalized) was significantly higher in the COV+ group compared with the COV- group (Table 1). The incidence of acute renal failure and acute respiratory failure during the present hospitalization was also significantly elevated in the COV+ group compared to the COV- group (Table 1). Furthermore, 31% of the COV+ individuals showed an abnormal arterial pH (acidosis or alkalosis) compared with 7.5% in the COV+ group. Of individuals in the COV+ group, 23.9% are now deceased. Of the COVID-19-related deaths, 76.5% occurred within 30 days, 96.1% within 60 days, and 98% within 90 days of study entry.

Study design

Blood was collected from study participants directly after enrollment in the WU-350 study (d0). Further longitudinal samples were collected 3 (d3), 7 (d7), 14 (d14), 28 (d28), and 84 (d84) days after the initial blood collection, when possible. The collection of longitudinal samples depended on survival of the study participants as well as the participants’ ability to complete study visits after being discharged from the hospital. A total of 703 human plasma samples from 341 patients were available for metabolomics profiling, including 324 d0 samples, 164 d3 samples, 111 d7 samples, 54 d14 samples, 31 d28 samples, and 19 d84 samples. All samples were divided into nine randomized sample batches and analyzed by liquid chromatography/mass spectrometry (LC/MS). An extract of the standard reference material (SRM) 1950 from NIST (National Institute of Standards and Technology, Metabolites in Frozen Human Plasma) was measured repeatedly as a quality control (QC), and blank samples were used to assess background signals. Polar and lipid metabolite fractions were extracted from each sample, and a global metabolomics profile was acquired in both positive and negative ion modes. Processing of the data led to the putative identification of 235 polar and 472 lipid metabolites based on accurate mass and MS/MS matching. Peak areas were extracted for those 707 metabolites to form the metabolic profile of each patient.

Given that the metabolic profiles were acquired over several months, the combined data showed strong batch effects as demonstrated by the principal component analysis (PCA) in Figure S2A. To remove the variance introduced by the individual batches, but not lose the differentiating biological variance within the research (WU-350) samples, we tested several normalization approaches (Figure S2B) and selected a combined batch correction (ComBat) approach that outperformed the other common normalization approaches tested (e.g., probabilistic quotient normalization, unit length, constant sum, quantile, etc.). After normalization, the metabolic profiles retained differences according to sample origin (WU-350, QC, and blank) as shown in Figure S2C but no longer clustered based on batch (Figure S2D).

The goal of this study was to find metabolic alterations that are predictive of disease severity in SARS-CoV-2-positive individuals. We used admission to the ICU during disease progression to classify patients as having severe or non-severe disease, as has been done previously., An ideal set of predictor metabolites would allow an individual presenting at the hospital and receiving a positive SARS-CoV-2 PCR-test result to be assessed for the likelihood of severe disease progression to best guide treatment at the earliest stage of hospitalization. Thus, we grouped the presented COV+ cohort into a non-severe (COV+ non-severe) group that did not require ICU admission and a severe group (COV+ severe) that did require ICU admission. For data interpretation purposes, three samples were excluded as the specimens originated from two patients who tested positive for other coronaviruses (NL63, HKu1). The final patient cohort consisted of 67 COV- cases, 143 COV+ non-severe cases, and 129 COV+ severe cases. Unsupervised analysis of the metabolic profiles for the 322 d0 samples available in our patient cohort demonstrated a clear trend in principal components space that separated COV+ severe, COV+ non-severe, and COV- patients (Figure 1A). Further, several significantly varying metabolites suggested that the metabolic profiles at d0 may indeed be predictive of disease severity. Hierarchical clustering analysis (HCA) of the 54 statistically significant metabolites (p < 0.05, Welch’s ANOVA) with an absolute fold change greater than two when compared with the COV- group revealed striking changes in multiple representatives of lipid classes, including lysophophatidylcholines (LPCs), phosphatidylcholines (PCs), and triglycerides (TGs). Further, several polar metabolites known to be related to COVID-19, including gluconate and dimethylguanosine, were also significantly altered (Figures 1B and S3).

Figure 1

Study design

(A) Principal component analysis based on all polar (n = 235) and lipid (n = 472) metabolites in SARS-CoV-2-negative individuals (COV-, n = 59, green), SARS-CoV-2-positive individuals with non-severe disease (COV+ non-severe, n = 140, orange), and SARS-CoV-2-positive individuals with severe disease (COV+ severe, n = 123, red). Data shown are only from the sample provided during study entry(d0).

(B) Hierarchical clustering analysis of metabolic profiles of COV-, COV+ non-severe, and COV+ severe patients at d0. Represented are 54 significantly changing polar and lipid metabolites (p < 0.05, Welch’s ANOVA, Benjamini-Hochberg correction). Each column is a sample, and each row is a metabolite. Cell color corresponds to the normalized metabolite intensity in each sample.

(C) Human cohort of 67 SARS-CoV-2-negative and 272 SARS-CoV-2-positive participants. The cohort was divided into a training cohort and a test cohort. The study design incorporated six blood draws for SARS-CoV-2-positive individuals on days 0 (d0), 3 (d3), 7 (d7), 14 (d14), 28 (d28), and 84 (d84) after WU-350 study entry.

Predictive model of COVID-19 disease severity

The global trends in the d0 metabolic profiles visible in the PCA and HCA visualizations prompted us to develop an ML model of disease severity that would predict ICU admission caused by SARS-CoV-2 infection. Patients presented to the hospital and enrolled in the study at different times during their disease course, hence the d0 sample is variable with respect to days- post symptom onset. The time from patient-reported symptom onset to d0 sample collection ranged from 0 to 107 days, with 4 days being the median. In addition, 45 of the COV+ patients were on mechanical ventilation and 98 COV+ patients were admitted to the ICU when the d0 sample was collected (Figure S4A). Although this diversity added additional variance to the d0 metabolic signatures, it also enabled a more clinically relevant model to be constructed because patients will present to the hospital at various points in the disease course. We note that COV+ severe and non-severe patients were not significantly different in time from symptom onset to the d0 sample collection (Figure S4B). The 707 metabolites that composed the metabolic profiles served as the predictors for our ML model. To assess predictive power, we split our dataset into two distinct groups of COV+ patients: a training set (163 patients) that we used to select, optimize, and train our ML model and an independent test set (100 patients) that was only used to evaluate the model’s performance (Figure 1C; Tables S1 and S2). Using our training set, we evaluated the efficacy of five ML algorithms with 20-fold cross-validation and found that a linear ElasticNet regression model was the most effective (Figure S4C). After training the model, we applied it to the patients in the test set and assessed performance by using the area under the receiver operating characteristic curve (AUC). On the test set, we saw strong predictive performance (AUC = 0.72) that outperformed a simple model that only uses BMI and age to predict disease severity (Figure 2A) and is significantly more predictive than a random model (Figure 2B, see permutation test in method details). As further validation, when the trained model was applied to the COV- patients (no COV- patients were in the training set), the mean scores output by the model were lower than those for the COV+ non-severe and the COV+ severe patients in the test set (Figure S4D). This indicates that the model can differentiate disease severity and distinguish COV+ and COV- patients. We wish to emphasize that PCR is the gold standard to diagnose SARS-CoV-2 infection. As such, we present this result only as confirmation that our model correctly predicts disease severity and not as a potential diagnostic for viral infection.

Figure 2

Predicting SARS-CoV-2 severity by machine learning

(A) Receiver operating characteristic (ROC) curves of prediction model on training set (green) and test set (blue). Random performance is shown in gray. ROC of BMI and age as predictors for severe COVID-19 (red) results in nearly random performance.

(B) Permutation test results from permuting training set labels and training the model on the permuted data. With every permutation, the area under the ROC curve (AUC) for the test set was computed. The histogram shows the distribution of these AUC values for 1,000 random permutations. The model performance on the test set when trained on the non-permuted training data results in an empirical p value of 0.002, as shown in blue.

(C) Variable importance in reduced ElasticNet prediction model for disease severity of SARS-CoV-2 infection in humans. Negative values are predictive of non-severe disease and positive values are predictive of severe disease. Variable importance is after the model is trained on the complete dataset.

(D) Profile plot of the normalized intensity (d0) of the prediction model metabolites grouped into COV- (control, n = 59), COV+ non-severe (n = 140), and COV+ severe (n = 123). The color of the lines reflects the mean intensity of each metabolite in the COV- patients.

(E) Boxplots showing predictor metabolite intensities (d0) in the COV-, COV+ severe, and COV+ non-severe groups. Box limits represent the quartiles of each sample group. Whiskers are drawn to 1.5× of the inter-quartile range.

We next sought to interpret which metabolites were most salient to the model’s predictions. First, we computed the variable importance of the model when trained on the complete dataset, which found 92 unique metabolites that contributed to the model’s predictions. Among that group of 92 compounds were metabolites that have been previously implicated in SARS-CoV-2 infection, such as kynurenate, nicotinamide, creatinine, LPCs, PCs, and others.,,, The mean intensity of each metabolite in the COV-, COV+ non-severe, and COV+ severe groups can be seen in Figure S5. Next, we aimed to assess the robustness of the metabolites selected by the ML model. We used bootstrap resampling of our training dataset to construct confidence intervals for the variable importance of each of the 707 metabolites we profiled. The analysis led to the identification of 22 predictor metabolites that significantly contributed to the model’s fit. The structural identities of these metabolites were rigorously confirmed (see method details). Strikingly, manyof the predictor metabolites were LPCs. Using the reduced predictor set, we re-trained and re-optimized our ElasticNet model on the training set and assessed the predictive power on our test set. Using the reduced predictor set resulted in nearly an identical AUC to that of the full set of metabolites (AUC = 0.70) and still performed better than a random model or a model that used only BMI and age as predictors (Figures S6A and S6B). To ensure that our model is indeed capable of predicting future disease severity, we additionally evaluated our model on the portion of the test-set patients who were not admitted to the ICU on or before when the d0 sample was collected. We again found that our model’s performance remained high (AUC = 0.71) and continued to outperform a model based only on BMI and age (Figures S6C and S6D). The variable importance of the predictor metabolites when trained on the entire dataset is shown in Figure 2C. The mean intensity of the metabolites in the COV-, COV+ non-severe, and COV+ severe groups is shown in Figures 2D and S7. In addition, the intensities of these metabolites among test-set patients and patients not admitted to ICU on d0 are provided in Figure S8. To exclude the possibility that the metabolite alterations were due to the application of therapeutics, we additionally removed all d0 samples when dexamethasone or remdesivir were administered on the same day as sample collection. Using that subset, all predictor metabolites were still significantly altered. All LPCs and PCs that contributed to the model, as well as serine, presented a significant downward trend in signal abundance with disease severity. Conversely, the other polar metabolites (kynurenate and 1-methyladenosine), phosphatidylethanolamines (PEs), and ceramides exhibited a significant upward trend in signal intensity (Figure 2E).

Demographics, laboratory values, comorbidities, and COVID-19 severity

After evaluating the efficacy of the ML model, we wished to deduce the relationship of our metabolite predictors to COVID-19 disease severity. We examined whether these metabolites were reflective of an underlying condition/risk factor for severe disease or if they were related to the disease progression of COVID-19. We addressed the former by asking whether any of the predictor metabolites correlated with demographic factors, laboratory values, current severity predictors, or individual patient comorbidities available for the patient cohort. A comparison of the COV+ non-severe and severe groups identified several significantly different parameters (Table S3). The COV+ severe group was significantly older than the non-severe group (Figure 3A), but there was no significant difference in BMI (Figure 3B). CO2 levels were not significantly altered between groups (Figure 3C), with values mostly being in the normal range. In contrast, there were significantly increased levels of the inflammatory marker C-reactive protein (CRP; Figure 3D). D-dimer, absolute neutrophil count, and neutrophil percentage were also increased (Figures 3E–3G). Absolute lymphocyte count and lymphocyte percentage were decreased (Figures 3H and 3I). These data indicate more severe inflammation in the COV+ severe group compared with the non-severe group and are consistent with reports from previous studies.,, Neutrophil recruitment has also been shown to be dysregulated in severe cases of COVID-19.38, 39, 40, 41 Lymphopenia, abnormally low levels of blood lymphocytes, has been found to correlate with disease severity in COVID-19 and even to be predictive of disease severity.42, 43, 44, 45, 46

Figure 3

COV+ patient parameters

(A–I) Demographics, comorbidities, and laboratory values of SARS-CoV-2-positive individuals grouped by disease severity (non-severe, severe) for age (A), BMI (B), CO2 (C), C-reactive protein (D), D-dimer (E), absolute neutrophil levels (F), neutrophil percentage (G), absolute lymphocyte levels (H), and lymphocyte percentage (I). Statistical significance (p < 0.05) was assessed by using a two-tailed Student’s t test with unequal variance for data shown in (A-I).

(J) Proportion of COV+ severe and non-severe patients with specific comorbidities and laboratory test results.

(K) Pearson correlation of listed demographic/laboratory results/comorbidities with abundances of the predictor metabolites.

(L) Pearson correlation of interleukin levels with abundances of predictor metabolites. Cells in the heatmaps (K and L) shaded white do not have a statistically significant Pearson correlation (p < 0.05, Benjamini-Hochberg correction).

(M) AUC values for patient comorbidities, demographics, laboratory values, the metabolite model, and a model that uses both the metabolites and laboratory values (absolute lymphocyte, absolute neutrophil, CRP, and D-dimer) when predicting disease severity on the test-set patients with complete information for all patient parameters.

Given that specific comorbidities increase the risk of having a severe case of COVID-19,,, we asked which co-morbidities were enriched in the COV+ severe group compared to the COV+ non-severe group (Figure 3J; Table S3). Diabetes, cancer, and chronic kidney diseases were recorded up to 1 year before hospital admission or up to 1 year before the d0 samples for non-hospitalized individuals. Although diabetes was significantly more prevalent in the COV+ group, the number of individuals with cancer and/or chronic kidney disease was not significantly different between the groups. Because acute respiratory failure and/or acute renal failure is a feature of severe COVID-19, the proportion of individuals suffering from either or both is significantly greater in the COV+ severe patients. Further, laboratory tests showed an increased proportion of individuals having an abnormal arterial pH (acidosis or alkalosis) in the COV+ severe group compared to the COV+ non-severe patients. The presented laboratory results (arterial pH, neutrophils, lymphocytes, D-dimer, CRP, and CO2) were obtained within 48 h of the first sample after study entry (d0).

Considering the number of significant associations in the patient parameters between COV+ severe and non-severe patients, we wanted to test whether any of our predictor metabolites significantly correlated with the clinical data after we controlled for disease severity (method details). To that end, we computed the Pearson correlation (for continuous parameters) or the point biserial correlation (for binary parameters) between each predictor metabolite and patient parameter (Figure 3K). The analysis revealed that many of the patient comorbidities were not significantly associated with the predictor metabolites. Smoking, BMI, and cancer were not significantly correlated with any of the predictor metabolites. Age and CKD had weak correlations with 1-methyladenosine (r = 0.32, 0.33), kynurenate (r = 0.20, 0.36), and serine (r = −0.19, −0.15). Kynurenate has been described previously as being associated with CKD. The most striking results are the significant positive correlations of the PEs (r = [0.20, 0.23], [0.26, 0.37]), ceramides (r = [0.14, 0.21], [0.26, 0.43]), and 1-methyladosine (r = 0.26, 0.20) with inflammatory markers (neutrophils and CRP) as well as the incidence of acute renal failure and diabetes. The PEs (r = [−0.07, −0.05]), ceramides (r = [−0.14, −0.12]), and 1-methyladenosine (r = −0.21) were also negatively correlated with lymphocyte levels, suggesting that these metabolite changes may be associated with inflammation. The LPCs, PCs, and serine showed the opposite trend as the PEs and ceramides. These molecules had positive associations with lymphocyte levels (r = [0.15, 0.34]) and negative associations with CRP (r = [−0.37, −0.23]) and neutrophil levels (r = [−0.28, −0.08]), suggesting that they too are associated with inflammation. We further correlated severity-associated interleukin levels that have been reported for a portion of our cohort in a prior publication (Figure 3L). Interleukin 6 (IL-6) levels, which have previously been shown to be a marker for disease severity,53, 54, 55 were inversely correlated with LPCs, PCs, and serine after controlling for severity (r = [−0.43, −0.26]). Similarly, 1-methyladenosine (r = 0.47), kynurenate (r = 0.44), Cer-NS d18:1_16:0 (r = 0.32), and PE 16:0_18:2 (r = 0.22) had significant positive correlations with IL-6, giving further evidence that many of the metabolite alterations are related to inflammation.

We then sought to assess the predictive power of our ML model relative to the predictive power of the patient comorbidities and laboratory values collected within 48 h of the d0 sample (CRP, D-dimer, neutrophil, and lymphocyte levels). For each patient parameter, we computed the AUC when predicting disease severity for each patient in the test set without missing values for any of those parameters (n = 68; Figure 3M). For all evaluated comorbidities, the model achieved a higher AUC. Conversely, the laboratory values (D-dimer, CRP, neutrophil, and lymphocyte levels) were all more predictive of disease severity than the metabolite model. However, when a new model was trained that used these laboratory values in combination with the biomarker metabolites, the performance of the combined model (AUC = 0.77) outperformed both the metabolite model (AUC = 0.66) and each laboratory value (AUC = [0.72, 0.73]), suggesting that the metabolite model is capturing complementary information to the laboratory values. In total, these results suggest that our predictor metabolites are indeed relevant to the pathogenesis of SARS-CoV-2 infection and not merely markers of other risk factors.

Longitudinal progression of predictor metabolites

To give further confidence that our predictor metabolites are associated with COVID-19 pathogenesis, we aimed to determine how the levels of these metabolites changed over the course of the disease progression. First, we considered the portion of the COV+ severe cohort that survived SARS-CoV-2 infection. We sought to determine the temporal behavior of their metabolic profiles as patients reached peak disease severity and after recovery. Accordingly, we analyzed the longitudinal metabolite abundances from individuals who had severe disease but survived and were discharged from the hospital. We compared their initial d0 plasma sample to the sample taken closest to the day of ICU admission, when the disease had progressed to peak severity. We also compared their initial d0 plasma sample to the last sample provided by the patient at or after hospital discharge. For several LPCs and one PC, a V-shaped trend was observed (Figure 4A). After the initial sample (d0), the levels of these metabolites dropped further as the disease worsened, but then began to return to d0 levels during recovery. The reverse trend was observed for Cer-NS d18:1_16:0. Its levels significantly increased until the patients were admitted to the ICU, but its levels subsequently dropped to below the initial d0 values in the final sample obtained.

Figure 4

Course of disease progression

(A) Prediction model metabolites that significantly vary in intensity as a function of disease progression for SARS-CoV-2-positive patients surviving severe disease (COV+ severe). d0 denotes the first sample after hospital admission, ICU denotes the sample collected closest to ICU admission, and the last sample is the final sample collected for the patient. Only patients for whom these time points were distinct samples were plotted. Statistical significance was assessed by a repeated-measures one-way ANOVA with a Benjamini-Hochberg correction.

(B) Principal component analysis showing the trajectory of the mean metabolic profile of the predictor metabolites in COV+ non-severe patients, surviving COV+ severe patients, and deceased COV+ severe patients. Deceased and surviving status was based on 90-day mortality. Patient samples were binned based on days post symptom onset (0–3, 4–6, 7–11, 12–18, 19–40, and 41–100). Samples taken when dexamethasone or remdesivir were administered were excluded. Asymptomatic patients were excluded. The mean d0 metabolic profile of COV- patients is shown in green.

These pronounced longitudinal trends in surviving COV+ severe patients raised the question of how the trajectory of disease progression (as marked by our predictor metabolites) differed among COV+ non-severe patients, surviving COV+ severe patients, and deceased COV+ severe patients. We also sought to compare the end points in these groups to the COV- d0 patients. We constructed representative metabolite profiles for the groups by using the levels of the predictor metabolites at different time windows relative to symptom onset (0–3, 4–6, 7–11, 12–18, 19–40, and 41–100 days post symptom onset) and performed a principal component analysis that enabled the trajectory of each group to be visualizedin two dimensions (Figure 4B). Samples from asymptomatic patients or patients administered the COVID-19-related therapies dexamethasone or remdesivir on the day of sample collection were excluded from this analysis. Notably, the analysis revealed three distinct trajectories with starting points that trended with disease severity. The groups then followed a common trajectory toward the COV- d0 sample. However, the COV+ non-severe patients recovered much quicker than the COV+ severe patients, reaching a similar metabolic profile to the d0 COV- patients in less than 40 days post symptom onset. The deceased COV+ severe patients had a similar metabolic profile during the 41–100 day window when compared to surviving COV+ severe patients. These results point toward distinct metabolic programs associated with the time course of infection and disease outcome. When examining the individual metabolite levels within the four groups at each time window, the deceased COV+ severe patients showed the same direction of dysregulation across the predictor metabolites as the surviving COV+ severe patients, but the magnitude of the perturbation was increased. In the 41–100 day window, less than 20% of the predictor metabolites were significantly altered, suggesting that metabolite levels return to control levels over time (Figure S9).

Next, we sought to compare the longitudinal progression of the predictor metabolites between the surviving COV+ patients. In the COV+ severe group, the LPC levels increased over the disease course to levels that were comparable to the COV- group (Figure 5A). In the COV+ non-severe group, the LPC levels recovered faster, reaching control levels in less than 40 days post symptom onset (Figure 5B). In total, 15 of the predictor metabolites showed a significant change (p < 0.05, Welch’s ANOVA) across the longitudinal time points in the COV+ severe group (Figure 5C). Of the LPCs, 11 were significantly increased over time. PC 38:6, PC 20:4_20:4, and PC 18:2_22:6 were also significantly increased with time. After initially being increased compared to the d0 sample of the COV- group, kynurenate showed a significant decreasing trend. Because of the lower sample numbers, the COV+ non-severe group had only six metabolites (five LPCs and one PC) that showed a significant trend (p < 0.05, Welch’s ANOVA; Figure 5D).

Figure 5

Longitudinal trends in COV+ patients

(A) Profile plot of the mean predictor metabolite intensities relative to d0 COV- samples (n = 59, gray) in symptomatic SARS-CoV-2-positive individuals with severe COVID-19 (n = 123, COV+ severe).

(B) Profile plot of the mean predictor metabolite intensities in symptomatic SARS-CoV-2-positive individuals with non-severe disease (n = 140, COV+ non-severe).

(C and D) Heatmaps showing relative mean intensity of predictor metabolites in longitudinal profiles of symptomatic COV+ severe patients (C) or COV+ non-severe patients (D). The mean COV- d0 profiles are included as the control for reference. Patient samples were grouped based on symptom-onset windows. Patient samples where dexamethasone or remdesivir were administered were excluded from the analysis. Asymptomatic patients were excluded. ∗p < 0.05. Statistical significance was assessed by using a one-way Welch’s ANOVA with a Benjamini-Hochberg correction.

Syrian hamster model confirms metabolite changes in COVID-19

Lastly, we aimed to validate that the predictor metabolites we found in the human samples were directly related to COVID-19 and not to some unidentified correlate that could have confounded our analysis. As such, we used an established hamster model of SARS-CoV-2 infection., Syrian hamsters have been found to be susceptible to SARS-CoV-2 infection, with the virus mainly replicating in the upper and lower respiratory tract of intranasally challenged animals for approximately 6 days post infection. Diseased animals exhibit a loss of body weight and pathological lung inflammation. We obtained plasma samples from golden Syrian hamsters that were intranasally inoculated with SARS-CoV-2 according to Imai et al., influenza virus according to Iwatsuki-Horimoto et al., or nasally treated with saline solution as a mock infection. Viral infection of the lung and body-weight loss were confirmed in all SARS-CoV-2- and influenza-infected hamsters and were consistent with previous reports., All SARS-CoV-2-infected hamsters experienced greater body-weight loss (>5.1% of their initial body weight) compared to influenza-infected hamsters, and virus was detected in both the lungs and nasal cavity in both models (see method details; Figure S10).

After 2, 4, 6, and 14 days (d2, d4, d6, and d14) post infection, plasma was harvested from the SARS-CoV-2- and influenza-infected hamsters (Figure 6A) because those time points correlate with the well-established kinetics of pathology and infection in the hamster model of SARS-CoV-2.,, For the mock group, plasma was harvested on days 4 and 14 relative to the infection timeline. All plasma samples were subjected to the same LC/MS workflow as described above. Differences in the sample matrix prevented a direct comparison of metabolite intensities between rodent and human plasma, but we are able to make comparisons between the different experimental groups. We compared the predictor metabolite levels of samples from the three groups (SARS-CoV-2, influenza, and mock) harvested on d4 as a representative time point for early disease. Figure 6B shows the 11 significantly changing predictor metabolites (p < 0.05, Welch’s ANOVA) across the three groups. Of the LPCs, eight were significantly altered in animals infected with SARS-CoV-2. Although PC 20:4_20:4 was not detected in the hamster plasma, possibly because of matrix effects, the two detected PCs were significantly dysregulated, along with PE 16:0_20:4. All 11 of the significantly altered metabolites showed the same direction of dysregulation as what was observed in the d0 human samples. For all significantly varying metabolites except PE 16:2_20:4, infection with the influenza virus showed a similar trend of dysregulation relative to control samples but a very different magnitude compared to SARS-CoV-2 infection. This is consistent with milder disease in influenza-infected hamsters compared to hamsters infected with SARS-CoV-2 and suggests that the depression of LPCs and PCs may not be as specific to SARS-CoV-2 infection as changes in PEs.

Figure 6

Syrian hamster model confirms SARS-CoV-2-dependent metabolite changes

(A) Experimental design of the Syrian hamster model. Hamsters (n = 3–6 per group) were infected through intranasal inoculation of SARS-CoV-2 (1e5 plaque-forming units [PFUs]), influenza virus (1e5 PFU), or nasally treated with a saline solution (mock) on day 0 (d0). After 2, 4, 6, and 14 days (d2, d4, d6, and d14) post infection, animals were sacrificed and exsanguinated.

(B) Comparing metabolite intensities among hamsters infected with influenza (n = 6), SARS-CoV-2 (n = 5), and mock (n = 6) on d4 shows many of the predictor metabolites are significantly altered in the hamster model (p < 0.05, Welch’s ANOVA). Box limits represent the quartiles of each sample group. Whiskers are drawn to 1.5× of the inter-quartile range.

(C and D) Metabolite changes during disease progression in SARS-CoV-2 (C) and influenza (D). Infected animals show a faster recovery during influenza. All groups are n = 6 with the exception of SARS-CoV-2 hamsters at d2 (n = 3) and d4 (n = 5). ∗p < 0.10, ∗∗p < 0.05. Statistical significance was assessed by using a two-tailed Student’s t test with unequal variance between d2 and d4 samples. All values were corrected with the Benjamini-Hochberg procedure.

Finally, we wanted to determine whether the predictor metabolites in hamsters showed longitudinal changes over the time course of infection. Indeed, when compared to the d4 control (mock infection) samples, there was a similar trend in hamsters infected with SARS-CoV-2 (Figures 6C and S11A) as in human COV+ samples. There was minimal time-dependent variation in the mock-infection samples, and none of the predictor metabolite levels were significantly different between d4 and d14 in those animals (Figure S11B). In SARS-CoV-2-infected hamsters, plasma LPC levels decreased on d4 compared to d2 and slowly recovered toward the control levels on d14. By comparison, for the influenza-infected animals, metabolite levels approached those of the control group more rapidly (Figures 6D and S11C). The similar patterns of alteration and recovery in LPC and PC levels across both hamsters and humans suggest that these lipid species are directly related to SARS-CoV-2 infection.

Discussion

The current study sought to predict COVID-19 disease severity based on the metabolic profiles of human plasma samples obtained early in the disease course, before the onset of critical illness. We applied untargeted metabolomics to profile a patient cohort of 339 individuals, which amounted to 700 study samples in addition to QC and method blanks. From the experiments, we putatively identified 235 polar metabolites and 472 lipid species. Using these compounds, we applied ML techniques to build a predictive model that accurately classified a patient’s disease severity from their d0 metabolic profile obtained at the time of initial study entry. Currently, risk assessment is primarily based on BMI and age. Even though we see a significant difference between the ages of patients in the non-severe and severe disease groups (p < 0.0001), the results of our study show that risk assessment based on a small panel of predictor metabolites is more reliable than age and BMI, thereby providing a better metric for resource allocation.

Although our linear ElasticNet model is relatively simple compared to other popular ML models that have been applied to metabolomics datasets previously (including artificial neural networks, support vector machines, or ensemble-based approaches, such as random forest), linear models can be easily interpreted and provide robust performance even when modeling innately non-linear biological systems., Notably, a previous study successfully used an ElasticNet model to predict disease severity from a multi-omic dataset. Other studies have used non-linear ML models, such as random forest, to predict disease severity from metabolomics, lipidomics, and/or proteomics profiles.,, While those studies found higher AUC scores than that of our model, they used considerably smaller patient cohorts than what our model was trained and evaluated on. Further, AUC scores are not directly comparable, unless they are calculated from the same dataset. When we tested non-linear models (support vector machines and random forest), we found worse cross-validated performance relative to ElasticNet (see Figure S4C).

Interpretation of our model led us to identify metabolites that predicted disease severity. Using a reduced predictor set, we were able to retrain our model and found similarly strong predictive ability. Of the predictor metabolites, most were LPCs and PCs that were more decreased in patients with severe COVID-19 when compared to patients with non-severe disease. The observation that LPC and PC levels are altered in patients with COVID-19 has been corroborated by multiple studies from around the world. Early studies of COVID-19 patients from China and Italy showed a decrease in PC levels and an increase in select LPCs,,, but recent studies from South America and Canada have shown that both circulating PC and LPC levels decrease in patients with COVID-19., Not only do our results support that PCs and LPCs are directly related to COVID-19 disease severity, they also build upon previous findings by showing that these metabolites are decreased early in the disease course and are, therefore, predictive of ICU admission. Additionally, upon disease recovery, we find that LPCs are restored to the levels of control individuals. We also confirmed that PCs and LPCs are decreased in a hamster model of SARS-CoV-2 infection and that the levels of those metabolites in animal plasma similarly recover over time.

The cellular mechanisms underlying the decreased LPC and PC levels in the plasma of patients with severe COVID-19 may be connected to uncontrolled inflammatory responses or viral pathogenesis. Multiple studies have shown that individuals with severe COVID-19 display a hyper-inflammatory state characterized by dysregulated production of cytokines, such as IL-1, IL-6, IL-10, and others.61, 62, 63 Fatty acids in plasma have also been correlated with sustained IL-6 release and renal dysfunction in a smaller cohort of patients with COVID-19, suggesting that altered inflammatory responses may affect levels of circulating metabolites. Consistent with prior studies, in our cohort, we found that PC and LPC levels in patients with COVID-19 correlated with inflammatory markers (e.g., IL-6, CRP, lymphocyte, and neutrophil levels). In addition, reductions in circulating LPC levels have also been found in patients with pneumonia-induced inflammation, infection of the lung, and sepsis in which a cytokine storm occurs. Alternatively, metabolic alterations may be directly related to the pathogenesis of SARS-CoV-2. Reprogramming of cellular metabolism is critical to sustain viral replication during coronavirus infection., SARS-CoV-2-infected monocytes upregulate genes involved in lipid biosynthesis and lipid remodeling, leading to the formation of lipid droplets. Additionally, a CRISPR-based screen for host factors required for the replication of coronaviruses identified genes involved in the sterol regulatory-element binding-protein (SREBP) pathway and the gene encoding the low-density lipoprotein (LDL) receptor. The SREBP pathway regulates cholesterol metabolism, and the LDL receptor influences circulating lipid levels. As such, the alterations in plasma metabolites found in our study may reflect systemic alterations in metabolism induced by sustained viral replication, especially in severe disease.

In summary, our large sample size that included longitudinal measurements of patient plasma, collection of patient metadata (laboratory values, comorbidities, and demographics), and use of an animal model allowed us to discover a small panel of metabolite predictors for COVID-19 severity and rigorously investigate the relationship of those metabolites to the pathology associated with SARS-CoV-2 infection. The model we present here to assess the risk of a severe case of COVID-19 does not require intensive computation or untargeted metabolomics, making it immediately applicable to most hospitals around the world. The metabolites presented could, for example, be quantitatively measured by triple quadrupole mass spectrometers, which are widely available in clinical laboratories.

Limitations of the study

Although the size of our patient cohort is a major strength of this study, sample collection started early in the pandemic (March 2020), when the appropriate treatment of SARS-CoV-2 infections was not yet established. Over the time of study enrollment, treatment plans changed (e.g., fewer individuals were treated by mechanical ventilation in the second half of the study). Additionally, although the goal was to collect longitudinal samples 3, 7, 14, 28, and 84 days after initial sample collection, we were unable to obtain samples from all patients at every time point. Subjects were less likely to provide plasma after being discharged from the hospital, and samples could not be collected after patients were deceased. Moreover, even though the d0 sample would ideally have been taken at or preceding hospital admission, this was not always possible. Some patients enrolled in the study several days after being admitted, or they were enrolled only after being transferred from another hospital. These practical challenges affected the statistical power of our temporal comparisons. Furthermore, it is important to note that the patients studied here were all from hospitals in the St. Louis region and were mostly collected before the emergence of variants, such as B.1.617.2 (delta). Further longitudinal studies of patients from other geographic regions who have been infected with variants are thus needed to confirm the broad applicability of our findings.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to the lead contact, Dr. Gary J. Patti (gjpattij@wustl.edu).

Materials availability

This study did not generate new unique reagents.

Resources used in this study are available to the scientific community upon request from Dr. Gary J. Patti (gjpattij@wustl.edu).

Experimental model and subject details

Human subjects and samples

Over the period of March to August of 2020, blood specimens of 341 individuals who presented at Barnes Jewish Hospital or Christian Hospital located in Saint Louis, Missouri, USA were collected. Inclusion criteria were a physician-ordered SARS-CoV-2 nasopharyngeal swab polymerase chain reaction (PCR) test with a positive or negative outcome, availability of gender and age information, and an age greater than 18. Informed consent was obtained from all study participants. The clinical cohort consisted of 155 female and 184 male participants. Out of the 339 patients, 272 were considered SARS-CoV-2 positive (COV+) and 67 were considered SARS-CoV-2 negative (COV-). Samples were collected at the time of enrollment (d0) to the study, and 3, 7, 14, 28, or 84 days post study entry. Clinically relevant medical information (e.g., patient-reported symptoms, date of symptom onset, age, race, and BMI) was collected at the time of enrollment from the subject or from the medical record. Data were retrieved for the current study on October 16th 2020.

Animal studies

All animal studies were performed at Mount Sinai School of Medicine. Outbred female LVG golden Syrian hamsters (6-8 weeks of age) were obtained from Charles River Laboratories (Kingston, NY). Experiments were performed as previously published., In short, the hamsters were anesthetized by intraperitoneal injection of a mixture of ketamine and xylazine prior to intranasal inoculation with 0.1 mL of 1e5 plaque-forming units (PFU) of SARS-CoV-2 (WA-1) or H1N1 influenza A virus (A/California/04/2009). On day 2, 4, 6, and 14 after infection, 3-6 anesthetized hamsters per infection group were euthanized by exsanguination followed by intracardiac injection of veterinary euthanasia solution (SleepAway; Fort Dodge). Plasma samples were exposed to germicidal UV-C light. A disease score was calculated based on virology and weight-loss data.

Study approval

Portions of the human study relevant to Barnes Jewish Hospital, Christian Hospital, and Washington University were reviewed and approved by the Washington University in Saint Louis Institutional Review Board (WU-350 study approval #202003085, and plasma metabolomics study approval #202004204). All animal studies were approved by the Institutional Care and Use Committee at Mount Sinai School of Medicine, following the humane care and use guidelines set by the institution.

Method details

Metabolomics sample preparation

Participant plasma, which had been stored at −80°C upon collection, was thawed on ice. A 50 μL aliquot was transferred onto a solid-phase-extraction (SPE)-system CAPTIVA-EMR Lipid 96-wellplate (Agilent Technologies) before addition of 250 μL of acetonitrile containing 1% formic acid (v/v) and 10 μM internal standard (consisting of uniformly 13C and 15N labeled amino acids from Cambridge Isotope Laboratories, Inc). The samples were mixed for 1 min at 360 rpm on an orbital shaker at room temperature prior to a 10 min incubation period at 4°C. Afterward, 200 μL 80% acetonitrile in water (v/v) was added to the samples. The samples were mixed on an orbital shaker (360 rpm) for an additional 10 min at room temperature. The samples were then eluted into a 96-deepwell collection plate by centrifugation (10 min, 57 x g, 4°C followed by 2 min, 1000 x g, 4°C). Polar eluates were stored at −80°C until the day of LC/MS analysis.

The SPE-plates were then washed twice with 500 μL 80% acetonitrile in water (v/v). Lipids still bound to the SPE-material were then released into a second elution plate, in two elution steps applying 2x 500 μL 1:1 methyl tert-butyl ether:methanol (v/v) onto the SPE cartridge and centrifuging for 2 min at 1000 g and 4°C. The combined eluates were dried under a stream of nitrogen (Biotage SPE Dry Evaporation System) at room temperature and reconstituted with 100 μL 1:1 2-propanol:methanol (v/v) prior to LC/MS analysis.

Hamster plasma samples were diluted 1:4 with methanol (v/v), vortexed for 30 s, and incubated at −20°C for 2 hours. Samples were centrifuged for 10 minutes at 13,500 x g at 4°C and supernatant was transferred to a new centrifuge tube, concentrated, and stored at −80°C until reconstitution as described above.

LC/MS analysis of polar metabolites

A 2 μLaliquot of polar metabolite extract was subjected to LC/MS analysis by using an Agilent 1290 Infinity II liquid-chromatography (LC) system coupled to an Agilent 6540 Quadrupole-Time-of-Flight (QTOF) mass spectrometer with a dual Agilent Jet Stream electrospray ionization source. Polar metabolites were separated on a SeQuant® ZIC®-pHILIC column (100 × 2.1 mm, 5 μm, polymer, Merck-Millipore) including a ZIC®-pHILIC guard column (2.1 mm x 20 mm, 5 μm). The column compartment temperature was maintained at 40°C and the flow rate was set to 250 μL·min-1. The mobile phases consisted of A: 95% water, 5% acetonitrile, 20 mM ammonium bicarbonate, 0.1% ammonium hydroxide solution (25% ammonia in water), 2.5 μM medronic acid, and B: 95% acetonitrile, 5% water, 2.5 μM medronic acid. The following linear gradient was applied: 0 to 1 min, 90% B; 12 min, 35% B; 12.5 to 14.5 min, 25% B; 15 min, 90% B followed by a re-equilibration phase of 4 min at 400 μL·min−1 and 2 min at 250 μL·min−1. Metabolites were detected in positive and negative ion mode with the following source parameters: gas temperature 200°C, drying gas flow 10 L·min−1, nebulizer pressure 44 psi, sheath gas temperature 300°C, sheath gas flow 12 L·min−1, VCap 3000 V, nozzle voltage 2000 V, Fragmentor 100 V, Skimmer 65 V, Oct 1 RF Vpp 750 V, and m/z range 50-1700. Data were acquired under continuous reference mass correction at m/z 121.0509 and 922.0890 for positive ion mode and m/z 119.0363 and 966.0007 for negative ion mode. Samples were randomized prior to analysis. In addition, a quality control sample was injected after every 12th sample to monitor signal stability of the instrument.

LC/MS analysis of lipid metabolites

A 2 μL aliquot of lipid metabolite extract was subjected to LC/MS analysis by using an Agilent 1290 Infinity II LC-system coupled to an Agilent 6545 QTOF mass spectrometer with a dual Agilent Jet Stream electrospray ionization source. Lipids were separated on an Acquity UPLC® HSS T3 column (2.1 × 150 mm, 1.8 μm) including an Acquity UPLC® HSS T3 VanGuard Pre-Column (2.1 × 5mm, 1.8 μm) at a temperature of 60°C and a flow rate of 250 μL·min−1. The mobile phases consisted of A: 60% acetonitrile, 40% water, 0.1% formic acid, 10 mM ammonium formate, 2.5 μM medronic acid, and B: 90% 2-propanol, 10% acetonitrile, 0.1% formic acid, 10 mM ammonium formate (dissolved in 1 mL water). The following linear gradient was used: 0-2 min, 30% B; 17 min, 75% B; 20 min, 85% B; 23-26 min, 100% B; 26 min, 30% B followed by a re-equilibration phase of 5 min.

Lipids were detected in positive and negative ion mode with the following source parameters: gas temperature 250°C, drying gas flow 11 L·min−1, nebulizer pressure 35 psi, sheath gas temperature 300°C, sheath gas flow 12 L·min−1, VCap 3000 V, nozzle voltage 500 V, Fragmentor 160 V, Skimmer 65 V, Oct 1 RF Vpp 750 V, and m/z range 50-1700. Data were acquired under continuous reference mass correction at m/z 121.0509 and 922.0890 in positive ion mode and m/z 119.0363 and 966.0007 in negative ion mode. Samples were randomized before analysis. In addition, a quality-control (QC) sample was injected after every 12th sample to monitor signal stability of the instrument.

Data preprocessing and normalization

Polar metabolite identifications were supported by matching the retention time, accurate mass, and MS/MS fragmentation data to our in-house retention time and MS/MS library created from authentic reference standards (Mass Spectrometry Metabolite Library supplied by IROA Technologies, Millipore Sigma, St. Louis, MO, USA) and online MS/MS libraries (Human Metabolome Database (HMDB, https://hmdb.ca ), Mass Bank of North America (MoNA, https://mona.fiehnlab.ucdavis.edu/ ), and mzCloud (https://www.mzcloud.org). Lipid iterative MS/MS data were annotated with the Agilent Lipid Annotator software. All data files were then analyzed in Skyline (Version 20.1.0.155) to obtain peak areas. m/z values of the metabolite and lipid target lists obtained from the metabolite identification workflow, which had at least an MS/MS match to an online library, were extracted under consideration of retention times.

Due to the risk of handling plasma samples from SARS-CoV-2-positive patients and not knowing how many batches of samples we would receive, we refrained from preparing a pooled sample and instead used the NIST SRM 1950 plasma reference material as a QC sample in each batch. The QC sample was injected after every 12th sample. After peak area extraction, batch effects were observed in the research samples (see Figure S2A). The research samples and QC data were used to test typical batch normalization methods (see Figure S2B), including constant sum, unit length, scale, percentile shift, minimum-maximum, probabilistic quotient normalization, quantile, and ComBat correction used in metabolomics.,72, 73, 74, 75 In Figure S2B, the variance remaining in the research samples normalized to the variance in the QC samples is shown for each method. The higher this ratio, the more variance that remains in the research samples and the more batch-derived variance in the QC samples is reduced. ComBat correction outperformed the other batch-correction approaches tested by using this metric. After correction, samples are well clustered according to sample type (WU-350, QC, blank) as shown in Figure S2C. In addition, within the research samples, there is no clustering by batch (see Figure S2D).

Machine learning

Samples were split into two distinct cohorts for training and testing the ML model. d0 COV+ patient samples within batches 1-6 made up the training set and d0 COV+ patient samples from batches 7 through 9 made up the test set. Training and tests sets were treated independently except for batch normalization, which was carried out for all patients (including samples collected after d0 and COV- samples) together. Demographics of both training and tests sets are available in Tables S1 and S2.

Model selection was based on 20-fold cross validation of the training set. Five different ML models: logistic regression, ElasticNet linear regression, partial least-squares discriminant analysis (PLSDA), support vector machine (SVM), and random forest were selected for consideration based on interpretability and previous studies.,,, Hyperparameters of all models and feature selection strategies were optimized by using 20-fold cross validation and a grid search. Two separate feature selection strategies were tested: a correlation-based approach and a statistic-based approach. In the correlation-based approach, the Pearson correlation was computed between each metabolite’s intensity and disease severity. Then, the top % of metabolites sorted by absolute correlation was taken as the predictors for the ML model. In the statistic-based approach, a Student’s t test was performed to assess the statistical significance of the differences in each metabolite’s intensity between COV+ severe and COV+ non-severe patients. Absolute fold-change and p value cutoffs were used to select metabolites. Performance was assessed with the area under the receiver operating characteristic curve (AUC). After optimization, ElasticNet regression achieved the highest AUC on the cross validated training dataset. The ElasticNet model is given below in Equation 1, where is the matrix of metabolic profiles (# of samples by # of metabolites), is the bias term, is the sample labels (0 = COV+ non-severe, 1 = COV+ severe), is the weight of each metabolite to the model prediction, is the weight of the regularization, and is the mixing parameter between the and norm regularization.After optimization, the correlation-based feature selection was used by taking the top 33% most correlated metabolites and model hyperparameters were set to and . In the reduced predictor model, no feature selection was performed and model hyperparameters were set to and .

The variable importance of each metabolite in the ElasticNet model is easily computed from the optimized weights, . To normalize for the different abundances of the metabolites, each weight was normalized by the median abundance of the metabolite across all samples. The more positive the variable importance, the more predictive that metabolite is to severe disease. The more negative the variable importance, the more predictive the metabolite is to non-severe disease. To find the metabolites that significantly contribute to the model fit, the training dataset was resampled with replacement 10,000 times. At each iteration, the ElasticNet model was trained and the variable importance was calculated. After the iterations were complete, the 95% confidence interval of the variable importance was calculated for each metabolite by using the 2.5 and 97.5 percentiles. If this interval included zero, the metabolite did not significantly contribute to the model fit.

All ML analyses were carried out by using Python (v3.7), with extensive use of the packages SciPy (v1.4.1) and Scikit-learn (v0.23.1).

Controlling for disease severity

To control for disease severity in the correlation analysis of the predictor metabolites, interleukin levels, and patient parameters (laboratory results, comorbidities, demographics), the mean value within each patient group (COV+ severe, COV+ non-severe, and COV-) was subtracted from each patient’s metabolite levels, interleukin levels, and patient parameters. This prevents weak correlations from occurring as a result of bulk changes associated with disease severity. After correction, the Pearson correlations between each corrected metabolite/interleukin measurement and patient parameter were calculated.

Confirming metabolite identities of predictor metabolites

The identities of the predictor metabolites were rigorously confirmed with authentic standards. For the polar compounds, authentic standards were purchased to not only match MS/MS spectra but also retention times. For lipids, one or two standards per lipid class were matched to an authentic standard to compare MS/MS spectra and retention times. PCs were identified based on m/z and the two characteristic fragments 184.0733 and 86.0964 in positive ionization mode. For PCs where no peaks for the acyl-chains were observed, only the sum composition can be given. PEs were identified based on the neutral loss of phosphorylethanolamine (141.0191) in positive mode. The fatty acyl composition could be derived from the spectra, but no differentiation of regioisomers was possible, as was the case for ceramides. To denote regiospecificity, metabolites whose regioisomers could be differentiated have their acyl-chains separated with a “/” while those that could not have a “_.” Cer-NS d18:1_16:0 was matched to its authentic standard. MS/MS data for Cer-NS d18:2_16:0 were matched to MS/MS library spectra. Cer-NS d18:3_16:0 eluted slightly before Cer-NS d18:1_16:0, as expected due to having one less double bond. LPCs were identified based on MS/MS spectral matches. Standards were available for LPC 14:0/0:0 and LPC 18:1/0:0. Their retention times were used as a reference for the other LPCs. The two regioisomers of LPCs (sn1 and sn2) were separated by liquid chromatography, with the sn1 isomer eluting later. They are also distinguished by their MS/MS spectra. 1-acyl-LPC (sn1) shows two main fragments (m/z 184.0733 and 104.1070), whereas the 2-acyl-LPC (sn2) has a more pronounced 184.0733 fragment. The 104.1070 fragment (choline) has been previously reported as being more abundant from sodiated LPCs. We note that sn2 LPCs can be converted to sn1 during sample preparation, and our sample preparation was not dedicated to preserve isomers.,

Acquiring MS/MS data

MS/MS spectra for polar metabolites and lipids were acquired by using an iterative approach in the MassHunter Acquisition Software (Version 10.1.48, Agilent Technologies) on an Agilent 6540 and 6545 QTOF, respectively. The same source settings as for MS1 data acquisition were used. MS/MS spectra were acquired at a scan rate of 3 spectra/s with different intensity thresholds and collision energies of 10, 20, and 40 V to increase identification rates.

To improve matching to Orbitrap spectral databases, MS/MS data for polar metabolites were acquired on an Orbitrap ID-X Tribrid mass spectrometer (Thermo Scientific). A Vanquish Horizon UHPLC system, with the same chromatographic conditions as described in the Method details, was interfaced with the mass spectrometer via electrospray ionization in both positive and negative ion mode with a spray voltage of 3.5 and 2.8 kV, respectively. The RF lens value was 35%. Data were acquired in data dependent acquisition (DDA) mode by using the built-in deep scan option (AcquireX) with a mass range of 67-900 m/z. MS/MS scans were acquired at 15K resolution on a NIST SRM 1950 plasma sample and from 4 individual samples (d0, d3, d7, and d14) in both positive and negative ion mode with different collision energies in the range of 20 NCE to 50 NCE for HCD and 30 NCE for CID to maximize identifications.

Quantification and statistical analysis

Statistical analysis

All statistical analyses were performed by using the SciPy (v1.4.1) and statsmodels (v0.11.1) Python packages and with the Mass Profiler Professional Software (Agilent Technologies, v15.5). All p values were corrected for multiple hypothesis testing by using the Benjamini-Hochberg procedure. All statistical tests used are described in the figure legends.

Permutation test

To assess the significance of the model fit and compare the predictive power to what is expected from random chance, we performed a permutation test. After the feature selection and model hyperparameters were optimized, the training dataset labels were permuted, and the model was retrained on the permuted data. Then, the performance of the model was assessed on the non-permuted test set and the AUC was computed. This process was repeated 1,000 times. The empirical p value was computed by calculating the percentage of the 1,000 permutations that achieved an AUC higher than that of the model’s performance when trained on non-permuted data.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Bacterial and Viral vectors
H1N1 influenza A virus (A/California/04/2009)	NIH	N/A
SARS-CoV-2 (WA-1)	BEI Resources	Cat# NR-52281
Biological samples
Human plasma samples	Washington University in St. Louis, Tissue Procurement Core (TPC), https://siteman.wustl.edu/research/shared-resources-cores/tissue-procurement-core-tpc/	WU-350
Reference human plasma	National Institue of Standards and Technolgoy	SRM1950
Chemicals, peptides, and recombinant proteins
Acetonitrile LC/MS OPTIMA	Fisher Scientific	Cat#A955 4, CAS 75-05-8
Euthanasia Solution SleepAway	Fort Dodge Animal Health	NAC No.: 10031902
Formic Acid for LC-MS LiChropur	Sigma-Aldrich	Cat#5438040100, CAS 64-18-6
Ketamine	Vedco, Inc., St Joseph, MO	NDC 50989-996-06
Xylazine	Akorn, Inc., Lake Forest, IL	N/A
Mass Spectrometry Metabolite Library	IROA Technologies	Cat#MSMLS
Medronic Acid	Sigma-Aldrich	Cat#64255, CAS 1984-15-2
Methanol LC/MS OPTIMA	Fisher Scientific	Cat#A456 4, CAS 67-56-1
Methyl tert-butyl ether HPLC Plus	Fisher Scientific	Cat#650560, CAS 1634-04-4
OPTIMA LC/MS 2-propanol (IPA)	Fisher Scientific	Cat#A461 4, CAS 67-63-0
Water LC/MS OPTIMA	Fisher Scientific	Cat#W64, CAS 7732-18-5
Deposited data
Raw and analyzed data WU-350	NMDR, https://www.metabolomicsworkbench.org	NMRD:ST001849, https://doi.org/10.21228/M80981
Metadata WU-350	NMDR, https://www.metabolomicsworkbench.org	NMRD:ST001849, https://doi.org/10.21228/M80981
Raw and analyzed data Syrian Hamster Model	NMDR, https://www.metabolomicsworkbench.org	NMRD:ST001853, https://doi.org/10.21228/M80981
Metadata Syrian Hamster Model	NMDR, https://www.metabolomicsworkbench.org	NMRD:ST001853, https://doi.org/10.21228/M80981
Experimental models: Organisms/strains
LVG golden Syrian hamsters	Charles River Laboratories (Kingston, NY)	LVG(SYR)
Software and algorithms
Biorender	BioRender.com	RRID:SCR_018316
GraphPad Prism v3.8	GraphPad	RRID:SCR_002798
Lipid Annotator v1 B1.0.54.0	Agilent Technologies	N/A
Mass Profiler Professional Software v15.5	Agilent Technologies	N/A
MassHunter Workstation LC/MS Data Acquisition v10.1 B10.1.48	Agilent Technologies	N/A
Python v3.7	python	RRID:SCR_008394
Scikit-learn v0.23.1	N/A	https://scikit-learn.org
SciPy v1.4.1	N/A	https://SciPy.org
Skyline 64-bit v20.2.0.343	MacCoss Lab Software	RRID:SCR_014080
statsmodels v0.11.1	N/A	https://statsmodels.org
Thermo Scientific Xcalibur v4.2.28.14, AcquireX workflow	Thermo Fisher Scientific	RRID:SCR_014080
Other
ACQUITY UPLC HSS T3 Column, 100Å, 1.8 μm, 2.1 mm X 150 mm	Waters	Cat#186003540
VanGuard Pre-Column HSS T3, 1.7 μm, 2.1mm x 5mm	Waters	Cat#186003976
SeQuant® ZIC®-pHILIC HPLC Column, 5μm polymer 100 × 2.1 mm	Merck Millipore	Cat#150462
SeQuant® ZIC®-pHILIC Guard Kit, 5μm polymer 20 × 2.1 mm	Merck Millipore	Cat#150438