Multi-omic longitudinal study reveals immune correlates of clinical course among hospitalized COVID-19 patients

Abstract

The IMPACC cohort, composed of >1,000 hospitalized COVID-19 participants, contains five illness trajectory groups (TGs) during acute infection (first 28 days), ranging from milder (TG1-3) to more severe disease course (TG4) and death (TG5). Here, we report deep immunophenotyping, profiling of >15,000 longitudinal blood and nasal samples from 540 participants of the IMPACC cohort, using 14 distinct assays. These unbiased analyses identify cellular and molecular signatures present within 72 h of hospital admission that distinguish moderate from severe and fatal COVID-19 disease. Importantly, cellular and molecular states also distinguish participants with more severe disease that recover or stabilize within 28 days from those that progress to fatal outcomes (TG4 vs. TG5). Furthermore, our longitudinal design reveals that these biologic states display distinct temporal patterns associated with clinical outcomes. Characterizing host immune responses in relation to heterogeneity in disease course may inform clinical prognosis and opportunities for intervention.

Keywords: COVID-19; SARS-CoV-2; immunophenotyping; longitudinal modeling; multi-omics; systems immunology.

Graphical abstract

Figure 1

Overview of IMPACC cohort, sample collection, and immunophenotyping (A) Clinical trajectory group assignment of IMPACC cohort participants (N = 1,164).^,, (B) The total number of collected tissues (whole blood, serum, plasma, and nasal swab samples) for all acute visits (up to day 28 post-admission, including escalation visits). A total of 15,193 samples were profiled from 540 participants across 20 hospital recruitment sites from 15 biomedical centers. (C) The total number of participants profiled by 14 different immunophenotyping assays over the course of the scheduled acute visits (visits 1–6). (D) Data analysis involved a rigorous data quality and confounder analysis, dimensionality reduction to combine features into modules, and association of module levels at visit 1 and their longitudinal pattern with the clinical trajectory group. Expression levels of modules at visit 1 are depicted as boxplots, while longitudinal patterns are shown as line graphs.

Figure 2

SARS-CoV-2 viral loads and antibody responses were associated with clinical trajectory group (A) Viral sequencing identified 60 PANGO lineages across the cohort. (B) The clinical trajectory group was not associated with any of the 9 lineages detected. (C) Viral loads (SARS-CoV-2 N1 gene Ct values) measured from samples collected at visit 1 (significantly higher in participants with more severe disease [adj.p = 0.037]). For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range. (D) Viral loads (SARS-CoV-2 N1 gene Ct values) from samples collected during the acute visits (shape: adj.p = 0.001, average: adj.p = 1.68e−5). (E) Anti-RBD IgG area under the curve (AUC) values measured from samples collected at visit 1 (lower in TG5 [adj.p = 0.68]). For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range. (F) Anti-RBD IgG AUC values from samples collected during the acute visits (shape: adj.p = 0.07, average: adj.p = 0.3). (G) Seroreactivity (log10 summed RPK across SARS-CoV-2 regions) across samples collected from the acute visits were measured longitudinally in two distinct regions (highlighted in gray within spike and N annotations): decreased seroreactivity in the NTD (shape: adj.p = 6.78e−6, average: adj.p = 0.058) and decreased overall seroreactivity in the LINK domain of the nucleoprotein (shape: adj.p = 0.023). (C and D) Because lower Ct values indicate higher viral loads, the y axis is reversed. (C and E) Shown are median values (horizontal lines), IQRs (boxes), and 1.5 IQRs (whiskers), as well as all individual points.

Figure 3

Association of serum proximity extension assay (Olink) and plasma proteomics modules with clinical trajectory groups (A–F) Analysis of serum Olink data identified significant associations in the expression levels of (A–C) Olink.mod3 and (D–F) Olink.mod2 among clinical trajectory groups. ImmuneXpresso, a text-mining tool linking cytokines/chemokines to cells, was used to annotate (A) Olink.mod3 (activator of cytotoxic NKs) and (D) Olink.mod2 (pro-inflammatory cytokines). (A and D) Significant enrichments (i.e., Fisher’s exact test p ≤ 0.05) are presented in the network. Blue arrows correspond to negative correlation/repression, while red arrows indicate positive correlation/production/activation. (B and C) Levels of Olink.mod3 (B) at visit 1 and (C) over time. (E and F) Levels of Olink.mod2 (E) at visit 1 and (F) over time. (G–O) Analysis of targeted and global mass spectrometry-based plasma proteomics data identified significant associations of (G–I) Targeted.Prot.mod1, (J–L) Targeted.Prot.mod3, and (M–O) Global.prot.mod4 with the clinical trajectory group. (G) MSigDB hallmark pathway analysis of the 58 proteins of Targeted.Prot.mod1 identified an association with coagulation. (H and I) Levels of Targeted.Prot.mod1 (H) at visit 1 and (I) over time. (J) MSigDB hallmark pathway analysis of the 26 proteins of Targeted.Prot.mod3 identified an association with coagulation and complement hallmark gene sets. (K and L) Levels of Targeted.Prot.mod3 at (K) visit 1 and (L) over time. Analysis of global mass spectrometry-based plasma proteomics data identified significant associations of Global.prot.mod4 with the clinical trajectory group. (M) MSigDB hallmark pathway analysis of the 54 proteins of Global.prot.mod4 identified an association with apical junctions, myogenesis, and epithelial mesenchymal transition. (N and O) Levels of Global.prot.mod4 at (N) visit 1 and (O) over time. (B, E, H, K, and N) For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range. (B, C, E, F, H, I, K, L, N, and O) Each point is a sample from an individual participant. Light gray lines connect samples from the same participant. Thick black lines correspond to a smooth spline fit for all participants in each trajectory group.

Figure 4

Association of plasma metabolomics modules with clinical trajectory groups (A–F) Analysis of plasma metabolomics data identified significant levels of (A–C) globalmet.mod6 and (D–F) globalmet.mod8 among clinical trajectory groups. (A–C) Levels of globalmet.mod6, comprised of mostly branched amino acid and urea cycle metabolites, (A and B) at visit 1 and (C) over time. (D–F) Levels of globalmet.mod8, which is comprised of phospholipid metabolites, were associated with severity at (D and E) visit 1 (adj.p = 7.33e−5) and (F) longitudinally. (B and E) For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range.

Figure 5

Association of cell subset frequencies with clinical trajectory groups (A) Visit 1 analysis identified the frequency of hematopoietic progenitor cells (HPCs) among non-granulocytes as different among clinical trajectory groups (adj.p = 6.34e−3), with higher average expression in the more severe groups. (B) The frequencies of CD14⁺CD16⁻, CD14⁺CD16⁺, and CD14^dimCD16⁺ monocyte subsets among parental monocytes at visit 1. (A and B) For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range. (C–E) Longitudinal analysis of (C) neutrophil subset frequencies (CD16hi, average adj.p = 9.76e−4, shape adj.p = 6.74e−3; CD16low, average adj.p = 0.0109, shape adj.p = 0.0310), (D) T cell frequencies (average adj.p = 6.01e−7, shape adj.p = 0.0123), and (E) cytotoxic NK cell frequencies among non-granulocytes.

Figure 6

Association of PBMC transcriptomic and nasal transcriptomic modules with clinical trajectory groups (A–I) Analysis of PBMC transcriptomic data identified 21 modules with significant differences in expression levels between clinical trajectory groups at visit 1, including (A–C) PBMC.mod2, (D–F) PBMC.mod14, and (G–I) PBMC.mod8. (A, D, and G) These modules were interpreted using the top 10 enriched terms by MSigDB Hallmark, Reactome, and KEGG pathway databases ranked by p value after filtering for significant pathways with p <0.05. (B, E, and H) Module expression over trajectory groups at visit 1. (C, F, and I) Module expression by trajectory group over time. (J–O) Analysis of nasal transcriptomic data identified 7 modules with significant differences in expression levels among clinical trajectory groups, including (J–L) module 3 (NasalRNAseq.mod3) and (M–O) module 4 (NasalRNAseq.mod.4). Networks of protein-protein interactions among genes in (J) module 3 and (M) module 4 were retrieved from STRINGdb. Size of a node denotes degree, and edge thickness denotes strength of interaction as provided by STRINGdb. (B, E, H, K, and N) For each boxplot, the vertical line indicates the median, the box indicates the interquartile range, and the whiskers indicate 1.5 times the interquartile range.

Figure 7

Markers of disease severity overlapping across assays Overlapping pathways associated with more moderate or more severe trajectory groups (A) at the time of hospitalization (visit 1; left) or (B) during the longitudinal follow up during the acute phase of the disease (right). For each overlapping pathway (row), the assays contributing to its identification as a marker of COVID-19 disease severity (column) are indicated. The color of each cell reflects whether the pathway is associated with moderate (blue) or severe (red) disease or both (purple). Pathways were manually separated into groups of biologically related processes based on their names.