Integrated longitudinal multiomics study identifies immune programs associated with acute COVID-19 severity and mortality

Abstract

BACKGROUNDPatients hospitalized for COVID-19 exhibit diverse clinical outcomes, with outcomes for some individuals diverging over time even though their initial disease severity appears similar to that of other patients. A systematic evaluation of molecular and cellular profiles over the full disease course can link immune programs and their coordination with progression heterogeneity.METHODSWe performed deep immunophenotyping and conducted longitudinal multiomics modeling, integrating 10 assays for 1,152 Immunophenotyping Assessment in a COVID-19 Cohort (IMPACC) study participants and identifying several immune cascades that were significant drivers of differential clinical outcomes.RESULTSIncreasing disease severity was driven by a temporal pattern that began with the early upregulation of immunosuppressive metabolites and then elevated levels of inflammatory cytokines, signatures of coagulation, formation of neutrophil extracellular traps, and T cell functional dysregulation. A second immune cascade, predictive of 28-day mortality among critically ill patients, was characterized by reduced total plasma Igs and B cells and dysregulated IFN responsiveness. We demonstrated that the balance disruption between IFN-stimulated genes and IFN inhibitors is a crucial biomarker of COVID-19 mortality, potentially contributing to failure of viral clearance in patients with fatal illness.CONCLUSIONOur longitudinal multiomics profiling study revealed temporal coordination across diverse omics that potentially explain the disease progression, providing insights that can inform the targeted development of therapies for patients hospitalized with COVID-19, especially those who are critically ill.TRIAL REGISTRATIONClinicalTrials.gov NCT04378777.FUNDINGNIH (5R01AI135803-03, 5U19AI118608-04, 5U19AI128910-04, 4U19AI090023-11, 4U19AI118610-06, R01AI145835-01A1S1, 5U19AI062629-17, 5U19AI057229-17, 5U19AI125357-05, 5U19AI128913-03, 3U19AI077439-13, 5U54AI142766-03, 5R01AI104870-07, 3U19AI089992-09, 3U19AI128913-03, and 5T32DA018926-18); NIAID, NIH (3U19AI1289130, U19AI128913-04S1, and R01AI122220); and National Science Foundation (DMS2310836).

Keywords: Adaptive immunity; COVID-19; Immunology; Innate immunity.

Figure 1. Data overview and multiomics factor generation.

(A) Table with the number of samples used in the integrative analysis separated by assay (rows) and scheduled time of collection (columns). Cells are shaded to reflect the relative number of multiomics samples available. (B) Plot of clinical TG assignments for all IMPACC cohort participants (n = 1,164) and clinical descriptions for each TG. The x axis represents days from hospital admission, and the y axis represents the ordinal respiratory status score (1, 2 = discharged, 3–6 = hospitalized, 7 = fatal). The shaded region denotes the IQR of each TG. (C) Preprocessed data for different assays were split into training and test cohorts. (D) Dimensionality reduction was performed via MCIA on the training cohort assays to construct multiomics factors and loadings. (E) Baseline factor scores were used to train a classifier for predicting the TG with the model selected via cross-validation. The classifier was used to predict the TG for the testing cohort factor scores. This figure was created with BioRender.com.

Figure 2. Multiomics factor prediction results.

(A) Prediction schema (created with BioRender.com). (B) Box plots of Spearman correlations and AUROC values for the severity task [TG1, TG2/TG3, TG4/TG5 vs. Probability(TG4/5)] and mortality task [TG4, TG5 vs. Probability(TG5|TG4/TG5)] across bootstrapped iterations for each baseline model for the testing cohort. Significance was calculated by standard normal approximation of bootstrapped differences between models. (*P ≤ 0.05 and ***P ≤ 0.0001). (C) Dot plot of coefficients for the MCIA prediction model, with dot size and color representing magnitude and direction, respectively. (D) Alluvial plot showing the distribution of testing cohort individuals in each TG linked to their initial baseline respiratory status. The line color represents the predicted risk from the MCIA model. (E) Dot plot of P values from linear mixed-effects models with enrollment site as a random effect. Sex and age were further adjusted as fixed effects when associating baseline factor scores with various baseline clinical measurements, complications during hospital stay, and comorbidities in the training cohort, with dot size and color representing significance and direction, respectively. Values are only shown for adj. P ≤ 0.05.

Figure 3. The severity factor increased in severe COVID-19.

(A) The severity factor scores across clinical TGs at baseline (severity adj. P = 1.4 × 10^–30, mortality adj. P = 0.049). (B) Longitudinal trajectory of the severity factor for different clinical TGs (mortality slope adj. P = 7.1 × 10^–15). The shaded region denotes a 95% CI from a generalized additive mixed model of the fitted trajectory (thick black line), thin black lines show individual participant-fitted models, and light gray lines connect the participants’ time points. (C) Pathway enrichment of the severity factor. (D) Network of enriched pathways and selected high-contribution features. The full list of associated features is in Supplemental Table 5. (E) Heatmaps of differential expression tests for pathways in C that showed baseline separation between TG4 and TG5 with linear mixed-effects modeling. (F) Heatmap of differential expression tests of leading-edge metabolites from metabolism pathways in E. (G) Regression coefficients from linear mixed-effects modeling between the severity factor and normalized cell frequencies from whole blood (CyTOF) of both parent and child populations. Daggers indicate a significant association between the reduction of a child cell population frequency, which is significantly associated with the severity factor and severity factor apoptosis signaling in PGX (mortality/severity = baseline mortality/severity task, slope5|4 = TG5 vs. TG4 longitudinally; *P ≤ 0.05,**P ≤ 0.01,***P ≤ 0.001; joint = aggregated P value across omics).

Figure 4. Integrative multiomics network identifies upstream regulators and mediators of NET formation.

(A) Broad elevation of transcriptomics and proteomics features in NET formation and complement in the severity factor. Pathway connections are from the Kyoto Encyclopedia of Genes and Genomes (KEGG) NET formation pathway. (B) Top cytokines in the severity factor, when bound to their receptors, trigger downstream signaling pathways, including ERK and p38 signaling pathways, and are important in NET formation.

Figure 5. Multiomics mortality factor enriched for antibodies, IFN signaling, and cellular metabolic changes.

(A) Mortality factor scores across clinical TGs at baseline (severity adj. P = 0.14, mortality adj. P = 0.049). (B) Longitudinal trajectory of the mortality factor for different clinical TGs. The shaded region denotes a 95% CI of the fitted trajectory (thick black line), thin black lines show individual participant-fitted models, and light gray lines connect the participants’ time points. (C) Functional pathway enrichment of the mortality factor revealed downregulation of Igs, upregulation of the IFN response, cholesterol metabolism, and acetylated peptides. (D) Network of enriched pathways in C and top selected high-contribution features. The full list of associated features is given in Supplemental Table 7. (E) Spearman correlation test between the mortality factor and serum anti–spike IgG antibody using baseline samples; P values were calculated from a linear mixed-effects model controlling for TG, sex, and age. (F) Regression coefficients from linear mixed-effects modeling of the mortality factor with normalized cell frequencies from whole blood (CyTOF) of both parent and child populations. Daggers indicate a significant association between the reduction of a child cell population frequency, which is significantly associated with the mortality factor and severity factor apoptosis signaling in PGX. (G) Differential expression tests via mixed-effects modeling of leading-edge metabolites in highlighted metabolomic pathways. (H) Spearman correlation coefficient between the mortality factor and nasal SARS-CoV-2 quantitative PCR (qPCR) Ct using baseline samples; P values were calculated from a linear mixed-effects model controlling for TG, sex, and age (mortality/severity = baseline mortality/severity task, slope5|4 = TG5 vs. TG4 longitudinally; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001; joint = aggregated P value across omics).

Figure 6. Virus-centered integrative multiomics network of the mortality factor.

(A) Positive association of nasal SARS-CoV-2 viral load and inverse associations of total and SARS-CoV-2–specific antibodies were top features of the mortality factor. (B) JAK/STAT IFN signaling was positively associated with the mortality factor and viral load, and IFN signaling inhibitors were also positively associated with the mortality factor, potentially contributing to the dysregulation of IFN responsiveness and uncontrolled viral load in TG5 despite early ISG elevation. The elevation of apolipoproteins from the plasma proteomics may also have contributed to the heightened STAT activity. Longitudinal trajectories from generalized additive mixed-effects modeling of (C) hallmark IFN-α response genes in the NGX, (D) IFN inhibitors in the NGX, and (E) nasal SARS-CoV-2 viral load determined from RT-qPCR Ct.

Figure 7. Summary of highlighted host immune programs.

(A) The severity factor identified an integrated multiomics cascade associated with disease severity, characterized by dysregulated metabolisms, e.g., essential amino acid (aa) metabolism, elevated inflammatory soluble proteins and transcripts, an elevated signature of coagulation and NETosis, and reduced T cell circulation and signaling in patients with more severe disease. The dysregulated metabolisms potentially served as early modulators of this broadly dysregulated immune state. The links in this panel reflect hypotheses formulated on the basis of our findings. (B) The mortality factor revealed a virus-centered multiomics immune state as an early hallmark of mortality among the critically ill patients (ICU, ventilation, or mortality), including reduced Igs and B cell circulation, dysregulated IFN responsiveness, as suggested by elevated IFN inhibitor levels in both nasal and PBMC transcriptomics, along with persistently elevated viral loads in patients with fatal illness. This figure was created with BioRender.com.