Longitudinal analysis reveals that delayed bystander CD8+ T cell activation and early immune pathology distinguish severe COVID-19 from mild disease.

The kinetics of the immune changes in COVID-19 across severity groups have not been rigorously assessed. Using immunophenotyping, RNA sequencing, and serum cytokine analysis, we analyzed serial samples from 207 SARS-CoV2-infected individuals with a range of disease severities over 12 weeks from symptom onset. An early robust bystander CD8+ T cell immune response, without systemic inflammation, characterized asymptomatic or mild disease. Hospitalized individuals had delayed bystander responses and systemic inflammation that was already evident near symptom onset, indicating that immunopathology may be inevitable in some individuals. Viral load did not correlate with this early pathological response but did correlate with subsequent disease severity. Immune recovery is complex, with profound persistent cellular abnormalities in severe disease correlating with altered inflammatory responses, with signatures associated with increased oxidative phosphorylation replacing those driven by cytokines tumor necrosis factor (TNF) and interleukin (IL)-6. These late immunometabolic and immune defects may have clinical implications. Crown

Introduction

The immune pathology associated with COVID-19 is complex (Wang et al., 2020; Zhou et al., 2020). Most infected individuals mount a successful anti-viral response, resulting in few, if any, symptoms. In a minority of patients, there is evidence that ongoing cytokine production develops, associated with persistent systemic inflammation, end-organ damage, and often death (Del Valle et al., 2020; Lucas et al., 2020). The relationship between the initial immune response to SARS-CoV-2, viral clearance, and development of the ongoing inflammatory disease that drives severe COVID-19 is not clearly established, nor have the kinetics of the immune changes seen in COVID-19 been fully assessed as disease progresses. Defective immune recovery might drive ongoing disease and perhaps contribute to secondary immunodeficiency with an increased risk of subsequent infection.

Severe COVID-19 is associated with profound abnormalities in circulating immune cell subsets. There is a decrease in many peripheral blood subsets of both CD4+ and CD8+ T cells but an increase in activated and differentiated effector cells (Arunachalam et al., 2020; Hadjadj et al., 2020; Kuri-Cervantes et al., 2020; Laing et al., 2020; Mann et al., 2020; Mathew et al., 2020; Su et al., 2020). Cells expressing programmed cell death protein-1 (PD1) and other inhibitory molecules are increased, though whether these reflect genuine T cell exhaustion or changes accompanying T cell activation has not been firmly established. There is, nonetheless, evidence of functional impairment in both CD8+ and CD4+ T cells in a number of studies (Chen and John Wherry, 2020). Data on T helper (Th) cell subsets are variable, but there is evidence of increased Th17 cells and markedly reduced T follicular helper (Tfh) cells (Chen and John Wherry, 2020; Mathew et al., 2020; Su et al., 2020). There have been conflicting reports regarding B cell immunity, but increased circulating plasmablasts (Arunachalam et al., 2020; Hadjadj et al., 2020; Laing et al., 2020; Mathew et al., 2020) and reduced germinal center responses (Su et al., 2020) are consistently observed in severe COVID-19. Innate T cell subsets, including gamma-delta (gd) T cells and mucosal-associated invariant T (MAIT) cells, are also reduced, as are non-classical monocytes and both plasmacytoid and myeloid dendritic cells (pDCs and mDCs) (Kuri-Cervantes et al., 2020; Laing et al., 2020)

By analyzing transcriptome, serum cytokine, and immunophenotyping data of longitudinal samples from COVID-19 patients with a range of disease severities, for up to 3 months from symptom onset, we were able to address two important questions regarding the immune response to SARS-CoV-2. (1) How does the very early immune response in patients who cleared the virus and recovered from disease with few or no symptoms compare with those who progressed to severe inflammatory disease? This provided insight into the features of the initial immune response that correlate with severe inflammatory outcomes and whether systemic inflammation is an early or later development in those who progress to severe disease. (2) How rapidly do the profound immune defects that accompany severe COVID-19 recover, and do the kinetics of recovery relate to ongoing inflammation and clinical status?

We recruited 207 patients with COVID-19, ranging from asymptomatic healthcare workers in whom SARS-CoV-2 was detected on routine screening through to patients requiring assisted ventilation, and compared their results to 45 healthy controls. We performed detailed immune phenotyping at multiple time points up to 90 days from symptom onset, reporting absolute cell counts rather than proportions. All data have been made available at https://www.covid19cellatlas.org/patient/citiid/. We found that the immune response in patients with progressive COVID-19 was delayed compared to those with mild disease and was inflammatory in nature from the outset. Early immune cellular changes predicted severe disease course, and the variable recovery of these cells over 3 months is associated with marked changes in the nature of the systemic inflammation seen in severe COVID-19.

Results

Patient cohorts

SARS-CoV-2 PCR-positive subjects were recruited for this study between 31st March and 20th July 2020. Those without symptoms, or with mild symptoms, were recruited from routine screening of healthcare workers (HCWs) (Rivett et al., 2020). COVID-19 patients were recruited at or soon after admission to Addenbrooke’s or Royal Papworth hospitals.

Study participants were divided into five categories of clinical severity, used throughout this paper unless otherwise stated (Figures 1A; Data S1). These were: (A) asymptomatic HCWs; (B) HCWs who either were still working with mild symptoms insufficient to meet the criteria for self-isolation (Rivett et al., 2020) or were symptomatic and self-isolating; (C) patients who presented to hospital but never required oxygen supplementation; (D) patients who were admitted to hospital and whose maximal respiratory support was supplemental oxygen; and (E) patients who at some point required assisted ventilation. Three patients who died without admission to intensive care were also included in this severe group. In addition, 45 healthy controls (HCs), distributed across a range of age and gender, were recruited.

Figure 1

Cohort characteristics and changes in inflammatory markers over time

(A) Study participant and sample numbers split by severity categories and 12-day time bins post screening (group A) or symptom onset (groups B–E).

(B and C) Distribution of participant age (B) and gender (C) across severity categories.

(D) Boxplots showing CRP (mg/L) and SARS-CoV-2 PCR cycle threshold in 12-day time bins. Gray band, interquartile range of HCs or the SARS-CoV-2 swab cycle negative threshold (CT > 38).

(E) Heatmap showing log2 fold change in median CRP and serum cytokine and complement measures between COVID-19 cases and HCs, 12-day time bins.

Wilcoxon test FDR adjusted p value: ∗p < 0.05, ∗∗p < 0.005, ∗∗∗p < 0.0005.

See also Figure S1 and Table S1.

We analyzed the immune phenotype of 605 blood samples from 246 participants out to 90 days from the onset of symptoms (Figure 1A; Tables S1 and S2). As the clinical severity category increased, patients were more likely to be older and to be male (Figures 1B and 1C), as expected (Wang et al., 2020). A high-sensitivity C-reactive protein (CRP) assay demonstrated that classifying disease severity on the basis of maximal respiratory support is reflected in the CRP (Figure 1D). Patient courses are measured in time since symptom onset for groups B through to E. As they are asymptomatic, those in group A are measured from the date of their first positive swab; they are therefore likely to have been sampled later post-infection than patients in groups B–E, and are therefore not directly comparable to them in terms of time course. CRP (and, later, other variables) is compared to the interquartile range of 45 healthy controls. Nasopharyngeal swabs were assessed for SARS-CoV-2, allowing inclusion in this study, and were repeated in some patients. Initial viral titers, reflected by low PCR cycle threshold (CT) values, were higher in group E. With only occasional exceptions, patients in all severity groups had cleared the virus by 24 days after symptom onset (Figures 1D and S1A). Of the six patients with positive swabs after 30 days, four were overtly immunosuppressed (three solid organ transplants with recent induction/rejection treatment, one myeloma on B cell depletion therapy) and one was a peritoneal dialysis patient admitted with peritonitis.

We will outline the major datasets collected in this study before integrating them to study early and recovering disease.

Changes in cytokines and complement components with time and disease severity

Asymptomatic HCWs in group A had no evidence of cytokine or complement dysregulation, while those with mild symptoms (group B) showed an early, transient increase in C3c and the terminal complement complex (TCC), but not in CRP or cytokine concentrations (Figures 1E and S1B). Once patients developed symptoms severe enough to attend hospital (group C or above), a different picture was apparent. Both interleukin (IL)-6 and tumor necrosis factor (TNF)-α were raised in serum, along with other cytokines, as were all complement components measured. These abnormalities were maximal at the first bleed and largely persisted in group E. Abnormal IL-6 and TNF-α persisted in groups C and D despite clinical improvement (all had been discharged by days 49–60). Interferon-gamma (IFN-γ) was briefly raised in only a subset of patients. Increased C3c was prominent in group B, while C3a became the dominant complement component elevated in more severe disease (groups C–E).

IL-6, TNF-α, IL-10, and IL-1β rose in those with more severe disease (groups C–E), but, in contrast, there was no increase in inflammatory cytokines in groups A and B, pointing to marked differences in very early inflammatory responses between resolving and progressive disease. In addition, the persistence of cytokine abnormalities even beyond 60 days from symptom onset could have implications for resolution of clinical disease.

Both onset and recovery of immune cellular abnormalities vary with disease severity

Using standardized flow cytometry panels, we explored the size of 64 cell populations over time across the five clinical strata. Trucount analysis enabled calculation of absolute cell numbers. Cellular changes with time were assessed (examples in Figure 2A), and outcomes for 30 cell types have been shown relative to the median for HCs (Figure 2B). CyTOF, which uses whole blood rather than peripheral blood mononuclear cells (PBMCs), was also used in a subset of patients to allow quantification of granulocytes and non-classical and intermediate monocytes (Figure S2; STAR Methods). See Data S2 for details of all cell types, including data beyond 48 days of symptom onset and time as a continuous variable.

Figure 2

Profound immune abnormalities in moderate/severe COVID-19

(A) Boxplots showing absolute counts (cells/μL) for two representative cell populations, by severity groups and 12-day time bins post screening (group A) or symptom onset (groups B–E). Gray band; interquartile range of HCs.

(B) Heatmap showing the log2 fold change in median absolute cell count (left) or proportion of major subset (right) between COVID-19 cases (samples, n = 362) and HCs (n = 45), 12-day time bins.

Wilcoxon test FDR adjusted p value: ∗p < 0.05, ∗∗p < 0.005, ∗∗∗p < 0.0005. Population hierarchy is shown to the left. Population proportions are calculated within parent populations listed to the right. PB or PBMC, peripheral blood mononuclear cells (flow cytometry); WB, whole blood (CyTOF).

See also Figure S3 and Data S2.

Few changes in the immune phenotypes were seen in patients with asymptomatic (A) and very mild (B) disease, but once symptoms warranted presentation to hospital, the picture changed. Widespread immune abnormalities in groups C–E were most marked at the first bleed (Figure 2B), even when this was within 0–2 days of symptom onset (Data S2). Almost all CD4+ T cells subsets were reduced, as were many CD8+ T cell subsets and both naive and memory B cells; in contrast, plasmablast numbers rose in all groups. Innate lymphoid subsets, including MAIT cells, various gd T cell subsets, and natural killer (NK) cells, were also reduced, as were mDCs, and both non-classical and intermediate monocytes. These changes were correlated with, and were predictive of, severity, as discussed below.

We also calculated leucocyte number as a proportion of “parent” populations, either of total PBMCs (Figure S3A) or of major lymphoid compartments (Figures 2B and S3A), to allow comparison with most published literature (and all single-cell studies). Considering proportion underestimates COVID-19-associated pathology, for example, missing the early severity-correlated reduction in lymphocyte subsets, and the late persistence of low numbers of T cell subsets in severe disease. In agreement with this, analysis of CyTOF and CITE-seq data, which also examine cell proportions, did not identify most alterations in cell populations between severity groups that were observed by flow cytometry (Figures S3B and S3C).

Blood transcriptomic inflammation-related signatures vary with severity and time

Whole-blood RNA was isolated, and transcriptomes were generated by RNA sequencing (RNA-seq) and analyzed (in two time bins: 0 to 24 days and 25 to 48 days) using PLIER, which performs matrix factorization to identify interpretable latent factors. The contribution to each latent factor by immune cell subsets was then calculated (Figures S4A and S4B). These RNA-seq-derived latent factors were broadly aligned with the pattern observed in the cell count data (Figure 2B). An exception was the pronounced neutrophil signature seen early across groups C to E and persisting at days 25–48 in group E, demonstrating more pronounced neutrophil dysregulation across severity categories than suggested by increasing neutrophil number alone. An erythrocyte gene expression-driven latent factor was also prominent in group E at late time points, may be associated with heme metabolism, and is discussed later.

We then used unbiased weighted gene correlation network analysis (WGCNA) in the whole-blood transcriptome data to identify modules of co-regulated genes, where each could be summarized as an “eigengene.” Prominent modules correlated with both disease severity and time (Figures 3A, 3B, S4C, and S4D; Table S3). The module enriched for TNF-α/IL-6 genes correlates well with the cytokine concentrations determined in Figure 1, rising early in groups C–E and then largely resolving by 25–48 days. A neutrophil activation module was also prominent early across groups C–E, as was one associated with glycolysis. Thus, there was clear transcriptional evidence of activation of broad inflammatory pathways at early time points, and these largely recovered in most patient groups (with the exception of group E, in which many patients had persistent disease). In contrast, an interferon-related module was upregulated in groups B–E at days 0 to 24 from symptom onset, before declining (Figures 3A). The relative contributions of type I, II, and III interferons to this interferon-stimulated gene (ISG)-associated module cannot be easily distinguished (Banchereau et al., 2017), but a more detailed analysis of its kinetics showed that, while expression peaked at different eigengene values in each severity group, it then declined in all of them by around 30 days (Figure 3C), coincident with viral clearance and occurring irrespective of clinical and inflammatory state (Figure 3D).

Figure 3

Whole-blood transcriptomic signatures over time (n = 183)

(A) Heatmap derived from WGCNA, correlating whole-blood co-expression gene modules (colored blocks, y axis) with COVID-19 severity groups (x axis) in 24-day time bins post screening (group A) or symptom onset (groups B–E). Boxes are colored by strength of correlation. For details of annotation (by Enrichr) and gene content of all modules, see Figure S4 and Table S3. Boxplots show eigenvalues within key transcriptomic modules, according to disease severity and time.

(B) Mixed-effects model with quadratic time trend showing the longitudinal expression of eigengenes over time by severity. Gray band, interquartile range of HCs. Nominal and adjusted p values for the time × severity group interaction term shown.

(C and D) Mixed-effects model showing longitudinal expression of eigengene capturing interferon-stimulated genes (ISG) (C) and equivalent mixed-model showing changes in SARS-CoV-2 PCR cycle threshold (viral load) by time and severity (D). y axis inverted in (D).

(E) GSEA enrichment for Hallmark genesets against HC in COVID-19 cases split by severity in 24-day time bins post screening (group A) or symptom onset (groups B–E).

FDR adjusted p value is shown by circle diameter, with color representing normalized enrichment score of the associated gene set. See also Figure S4 and Table S3.

Finally, a supervised gene set enrichment analysis (GSEA) performed using Hallmark gene signatures (Figure 3E) (Liberzon et al., 2015) showed results largely consistent with the unbiased approaches just described, including demonstrating a late upregulation of genes associated with reactive oxygen species (ROS) and oxidative phosphorylation (OXPHOS) discussed below in the context of immune recovery.

Immune phenotype at presentation correlates with severity and may predict outcome

To determine whether the immune phenotype at presentation correlated with, or indeed could predict, subsequent disease course, we first performed a principal-component analysis using cell numbers across 24 primary immune cell populations from the 84 participants with blood draws taken between 0 and 10 days after the development of symptoms. Those in groups A and B clustered with HCs and away from those in groups C–E (Figure S5A). Hierarchical clustering of absolute cell counts identified two clusters (Figure 4A), one almost entirely comprised of HCWs from groups A and B (cluster 2) and the other containing all patients who progressed to ventilation and/or death and most who required supplementary oxygen support (cluster 1). Clustering using RNA-seq data obtained from 1–10 days after onset largely recapitulated that observed using cell number and was driven by ISG-, TNF-α-, and IL-6-associated gene pathways (Figure S5B). The severe cluster 1 was associated with increased age, CRP, TNF-α, and IL-6 (Figure 4B; Data S3A). Early differences between cell types drove this clustering, irrespective of their different subsequent trajectories (Data S3B). We used a sparse partial least-squares discriminant analysis (sPLS-DA) to determine which cell subsets were most informative for cluster prediction: clusters 1 and 2 could be discriminated with a minimum classification error rate of 0.07 ± 0.02 (93% accuracy) based on 13 key cell populations selected by the model (Figures 4C and S5C–S5E). The area under the receiver operator characteristic (AUROC) curve for patient cluster classification based on these 13 cell types was 0.98 (98% chance of accurate cluster prediction) (Figure S5F). CRP alone was inferior to the cell types in classifying patients and did not improve their performance when added (Figure S5G). These cell types were often the most profoundly and persistently affected by severe COVID-19 (particularly MAIT, gd T, Tfh, and CD4+ Temra cells).

Figure 4

Multivariate analysis of immune-cell populations in early disease correlates with clinical outcome

(A) Unsupervised clustering of absolute cell counts across 24 cell populations (normalized to the median of HCs) for COVID-19 samples taken ≤10 days from screening (group A) or symptom onset (groups B–E). Cases group into two clusters (cluster 1, orange, n = 46; cluster 2, purple, n = 38) by Euclidean distance and Ward D hierarchical clustering.

(B) Boxplots comparing age and inflammatory characteristics of individuals in clusters 1 and 2 at the time of sampling.

(C) Thirteen cell types selected by sPLS-DA as most informative in predictive models discriminating patients in clusters 1 and 2. Bars indicate loading coefficient weights of selected features (ranked from most to least informative in cluster prediction, from bottom to top).

(D) Addition of age, CRP, serum cytokine, and complement measures to unsupervised clustering of cellular data in (A) results in tighter grouping of COVID-19 patients by severity (cluster 1, orange, n = 17. cluster 2, purple, n = 38).

(E) AUROC curve showing sensitivity and specificity of severity group prediction (derived from clustering in D), based on absolute counts of 24 key cell types, CRP, or serum measures alone compared to all available measures.

(F) Kaplan-Meier plot of escalation-free survival in individuals within severity cluster 1 or cluster 2 split by hospitalization status. Escalation defined as a step up in respiratory support or death.

p value for the chi-square test of the difference between cluster 1 (n = 17) and hospitalized patients in cluster 2 (n = 13) is shown; numbers denote non-escalated patients in each group from days 0 to 30 post symptom onset.

See also Figure S5 and Data S3.

The ability to cluster patients and separate those with mild or no symptoms (groups A and B) from those presenting to hospital (groups C–E) underlines the profound association of immune subset abnormalities with disease severity. It does not, however, predict outcome in a clinically useful way, as such predictions are only of practical value in those who present for medical attention. We therefore incorporated additional data related to inflammation, including cytokines, age, CRP, and complement components; re-clustering with this combined dataset led to a smaller severity-associated group (Figure 4D). Input from both immune cell subsets and inflammation-related markers was required for optimal classification (Figure 4E). When considering only the subset of patients recruited at hospital presentation (groups C–E) and bled within 10 days of symptom onset, this clustering predicted disease progression, as defined by the subsequent need for increased respiratory support, or death, after blood sampling (Figure 4F). This analysis requires validation with larger patient numbers but is comparable to similar observations made by others (e.g., Laing et al., 2020; Mathew et al., 2020) and indicates that a combination of immune phenotype data and inflammatory markers could provide potentially clinically useful prediction of disease progression.

Early immune responses correlate with COVID-19 severity

HCWs in groups A and B did not progress to severe COVID-19 disease and also clustered apart from those with more severe disease when using either early immune cell counts or transcriptome data (Figures 4A and S5B). Therefore, we compared the early features of immune responses in groups A and B to those in patients with more severe COVID-19, binned in 7-day intervals to provide finer definition, to seek early immune and inflammatory correlates of disease severity (Figure 5A).

Figure 5

Early immune changes associated with mild or severe disease outcome

(A) Heatmap showing the log2 fold change in median absolute cell counts, CRP or complement measures between COVID-19 cases and HCs by severity and in 7-day time bins post screening (group A) or symptom onset (groups B–E). Wilcoxon test FDR adjusted p value: ∗p < 0.05, ∗∗p < 0.005, ∗∗∗p < 0.0005.

(B and C) Mixed-effect model with quadratic time trend showing cellular trajectories over time in sample groups B, D, and E in non-naive HLA-DR+CD38+ CD8 T cells (B) or plasmablasts (C) (cells/μL) from weeks 1–3 post symptom onset (samples, n = 207).

(D) Number of CD3+ T cells secreting IFN-γ spontaneously or following SARS-CoV-2 antigen stimulation in samples from groups B (n = 22) and D and E combined (n = 25), 1 or 2 weeks post symptom onset. Kruskal-Wallis test p values: ∗p < 0.05.

(E) Log2 fold change (FC) in expression of CD8+ T cell transcripts reflecting T cell activation and surface protein expression (detected by antibody staining) in CITE-seq data from non-naive CD8+ T cells from patients in groups A and B (n = 5) and C, D, and E combined (n = 13), relative to HCs (n = 11).

(F) Normalized gene set enrichment score for gene sets associated with TCR-dependent and bystander T cell activation in single-cell transcriptomic data from non-naive CD8+ T cells from patients in groups A and B versus C, D, and E. FDR adjusted p value shown by circle diameter.

(G) Area under the curve for SARS-CoV-2 spike-specific IgG titers at 1, 2, and 5 weeks post screening (group A) or symptom onset (groups B–E). Groups C, D, and E are combined to increase statistical power, Kruskal-Wallis test p values annotated as in (A).

(H) SARS-CoV-2 antibody titers achieving 50% neutralization (NT50) in patients from groups B–E in the first 2 weeks post symptom onset (n = 102). Samples with no detectable neutralizing activity at the lowest dilution (dotted line) are plotted at an arbitrary NT50 of 1. p value and Pearson’s correlation shown.

(I) Boxplots showing SARS-Cov-2 viral load, taken as first positive swab PCR CT, in severity groups. Wilcoxon test p values annotated as in (A).

(J) Schematic summarizing variation in immune features of SARS-CoV-2 infection across cases of varying disease severity.

See also Figure S6.

A number of key features emerged. First, there was no evidence of systemic inflammation in groups A and B. CRP was normal (Figure 5A), cytokines were not raised (Figure 1E), and there was no RNA evidence of systemic inflammatory gene-related signatures (Figures 3A and 3B). Exceptions were transient increases in C3c and TCC. Most cell types that were profoundly reduced in groups C–E were normal in A and B (Figures 2 and S3), but some showed mild reductions in group B; pDCs and memory B cells have been shown as examples (Figures 5A–5C, S6A, and S6B). The fall in pDCs was modest and brief compared to severe disease and was consistent with tissue localization and local interferon production to support the antiviral response (Cella et al., 1999).

In group B, cytotoxic CD8+ T cells rose earlier than in groups C–E, apparent by day 7 and peaking up to 2 weeks after symptom onset, in contrast to the later and more sustained rise seen in severe COVID-19 (Figures 5A, 5B, and S6A). Early enrichment of a CD8+ cytotoxic RNA signature was also seen in group B compared to C to E (Figure S6C). Consistent with these findings, spontaneous generation of IFN-γ by T cells was more pronounced in group B samples taken 2 weeks post symptom onset (Figure 5D).

Antigen-specific CD8+ T cell responses were determined after stimulation with anti-SARS-CoV-2 peptides and subtraction of background spontaneous IFN-γ-producing T cells. They were similar in all severity groups early (Figures 5D and S6D) and so were not responsible for the increased effector CD8+ T cells seen in mild disease, suggesting instead a role for “bystander” activation. Such bystander-activated cells are important in early anti-viral defense (Maurice et al., 2021) and characteristically express NKG2D, involved in the killing of infected cells, and both IL-7 receptor (IL7R) and CD8, which are downregulated in a T cell receptor (TCR)-dependent fashion (Kaech et al., 2003; Slifka and Whitton, 2000). CITE-seq protein expression analysis of activated CD8+ T cells demonstrated an increase in surface NKG2D, IL7R, and CD8 in group B compared to more severe groups (Figure 5E), observations confirmed by RNA analysis (Figure 5E). Bystander CD8+ T cells expressing CXCR3 rapidly home to sites of inflammation (Maurice et al., 2021), consistent with enrichment for CXCR3 RNA, but not surface protein (Figure 5E). Transcriptional signature-derived bystander-activated CD8+ T cells were enriched in patients with mild disease, while those from TCR-activated cells were associated with severe COVID-19 (Figure 5F), again consistent with widespread early bystander activation in the CD8+ T cell population in patients destined to have good disease outcomes.

An early increase in plasmablasts was seen in groups A and B, occurring up to a week before a larger rise in more severe disease (Figures 5C and S6B). We therefore measured total immunoglobulin concentrations and anti-spike IgG and anti-SARS-CoV-2 neutralising antibody titers (Figures 5G, 5H, S6E, and S6F). Group B patients maintained their serum IgM concentrations, which fell rapidly in those with more severe disease (Figure S6G), and their titers of anti-spike IgG and early neutralisation responses were comparable to patients progressing to more severe COVID-19 (Figures 5G and 5H). This suggests that the B cell difference between patients with mild and more severe disease might lie in more robust non-antigen-specific B cell activation, perhaps impacting via “natural” antibody and/or non-antibody-dependent mechanisms.

Virus at first swab, as assessed by PCR CT value, was comparable in groups B, C, and D (and the lower titer in group A samples was not comparable as they will have been taken later after infection, as described above). Initial viral titer was therefore not associated with an increased risk of hospital admission (being similar in groups B, C, and D) but was higher in group E (Figure 5I). These viral titers were reflected in interferon-related transcription signatures, which are prominent in groups B–E (Figures 3A, 3C, and S6H). Despite the fact that high IFN signatures correlated with severe disease, they were not necessarily driving that severity. The fact that the subgroup with the highest IFN within group E do best (Figure S6I) suggests that the inflammatory pathways that create severe disease are distinct from IFN (which is also consistent with comparable kinetics of reduction of IFN signature regardless of severity; Figure 3C) and that a robust IFN response may be beneficial in this context, though this needs confirmation in a larger dataset.

Taken together, these data suggest that an early adaptive immune response is prominent in individuals who are asymptomatic or have mild disease, characterized by a rapid production of activated bystander CD8+ T cells, plasmablasts, and likely pDC tissue localization before antigen-specific responses become apparent. This appeared a more important correlate of severity than viral titer, which only became relevant in those progressing to ventilation or death (Figure 5J).

In those with more severe disease (groups C–E), evidence of systemic inflammation was present from the first blood test (Figure 5). If we focused on the 16 patients in groups C–E sampled between 2 days before and 4 days after symptom onset, 15 had a CRP >10 mg/L and/or neutrophil activation eigengene >0. All five patients sampled between 2 days before and 2 days after symptom onset met these criteria. In contrast, this was not seen in group A or B despite the fact that they mounted a more prominent early cellular response. It was not clear whether inflammation seen in C–E was causally related to the poor early cellular response seen in these patients. Such inflammation clearly did not, however, develop later from the progression of a non-inflammatory immune response or as a result of failure to clear virus and suggested that the inflammatory die is cast by the time symptoms appear; thus, strategies to prevent it and reduce its clinical impact would need to be established very early (Figure 5J).

Distinct patterns of immune recovery in COVID-19

In contrast to groups A and B, cellular changes in groups C–E were profound and usually most prominent at the first bleed (Figure 2), so determination of change in immune cell subsets over time was likely to be most informative in these groups. Therefore, we explored cell kinetics in groups C–E, assigning patients to two categories based on whether their CRP concentrations remained elevated above 10 mg/L (“persisting CRP”) or fell below 10 mg/L (“resolving CRP”) by their final bleed within 3 months post symptom onset (Figure 6A). The latter group included both individuals with early high CRP that then fell, together with those for which CRP remained low (10 mg/L) throughout. Changes in CRP over time differed between these two groups when assessed using a mixed-effects model, with time modeled as a continuous variable (Figure 6B).

Figure 6

Cellular and transcriptional trajectories in persisting and resolving disease (n = 263)

(A) CRP (mg/L) from groups C, D, and E grouped by persisting and resolving CRP.

(B) Mixed-effect model with quadratic time trend showing log2(CRP) trajectories in both patient groups, and the likelihood-ratio test p value for the time × group interaction term. Gray band, interquartile range in HCs.

(C) Heatmap showing the log2 fold change in median absolute cell count between COVID-19 cases in groups C, D, and E, split according to persisting or resolving CRP, and HCs. 12-day time bins. Wilcoxon test FDR adjusted p value: ∗p < 0.05, ∗∗p < 0.005, ∗∗∗p < 0.0005. Rate of cell number changes shown by lollipop plot; faster rate of recovery, or deviation from normal, is indicated by increasing stem length. Points are colored by log2 fold change in median absolute cell counts from HCs at 0–12 days; black outline indicates failure to recover to HC numbers within 60 days (defined in STAR Methods).

(D) Mixed-effect models showing longitudinal trajectories of gene module eigenvalues capturing neutrophil degranulation, interferon-stimulated genes, heme metabolism, and oxidative phosphorylation in CRP groups, p values reported as (B).

To compare cellular changes over time between persisting and resolving CRP patient groups, a “rate of change” for each cell population was calculated over 60 days post symptom onset. In brief, this rate captured both the initial deviation in cell counts from normal within a window of 0–12 days and the time taken for cells counts to stabilize within a normal range if cellular recovery did occur (see STAR Methods). Five predominant trajectories were observed; populations that did not deviate from heathy numbers over the duration of study (e.g., NKT cells), those that increased progressively from normal over time (e.g., effector CD8+ T cells), those that fell progressively from normal over time (e.g., transitional B cells), those that trended toward recovery after an initial rise in numbers (plasmablasts), and those that tended toward recovery after an initial drop in numbers (e.g., naive CD4+ T cells) (Figure 6C).

The absolute number of most cell populations fell precipitously early and then showed variable recovery. For descriptive purposes, these were arranged into a group of cell subtypes that failed to recover, or recovered, in the persisting CRP group (Figure 6C, quadrants I and III, respectively) and their equivalents in the resolving CRP group (Figure 6C, quadrants II and IV, respectively). Notably, a few populations remained abnormal in both persisting and resolving CRP groups out to 60 days post symptom onset, including pDCs, Tfh, MAIT, and Vg9Vd2 (hi) gd T cells. All other populations showing an early drop in counts (with the exception of naive CD8+ T cells) recovered in those with resolving CRP (II and IV) and at rates more rapid than seen in those with persisting high CRP. In the persisting CRP group, a number of cell types remained markedly abnormal (including memory B cells and various CD4+ T cell subsets: quadrant I), whereas a second group of cell types recovered despite persisting inflammation (including NK cells and some CD8+ T cell subsets: quadrant III).

We then explored the relationship between cell recovery and the inflammatory response. It might be expected that where CRP remained persistently elevated, immune defects might persist, if these defects were secondary to the inflammatory state. Consistent with this, the cohort with persistently raised CRP also has raised TNF-α and IL-6 protein (Figure 1D). Likewise, transcriptional signatures of TNF-α/IL-6 and neutrophil activation were increased in severe disease (Figure 3), particularly in the persistent CRP group (Figure 6D). This ongoing inflammation may contribute to the sustained reduction in cell numbers at late times seen in quadrant I, together with persistently raised HLA-DR+/CD38+ effector T cells and plasmablasts (Figure 6C). Therefore, it was also not unexpected that most cell types reduced in acute disease recover over a few weeks as the CRP falls, as was seen for most cells in quadrants II and IV (Figure 6C).

More intriguing were the cells that recovered rapidly in the face of ongoing inflammation (quadrant III). While the reasons for this are likely to differ between cell types and to be multifactorial, these cell reductions might be driven in part by the viral infection per se and/or virus-induced interferon. It was notable that, after initially rising, IFN-γ (Figures 1E and S1B) and ISG signatures fell to normal values independent of both disease severity group (Figure 3C) and CRP (Figure 6D) but correlated with declining virus titer (Figures 3D and S4E). Thus, cell types known to leave the circulation due to interferon stimulation (Kamphuis et al., 2006), such as T and NK cells (Hirsch and Johnson, 1986; Zafranskaya et al., 2007), may recover as viral infection is controlled and interferon-dependent inflammation falls, independent of ongoing CRP-associated inflammation.

Finally, a small number of cell types remained statistically abnormal after 60 days, even in the resolving CRP group. Thus, effector CD4+ and CD8+ T cells and plasmablasts remained elevated, and pDCs, Tfh, MAIT, and Vg9Vd2-expressing gd T cells remained reduced (Figure 6C) and were among those cells most predictive of poor prognosis (Figure 4C). These abnormalities persist despite resolution not only of CRP but of neutrophil degranulation, TNF-α/IL-6, and glycolysis-related signatures (e.g., Figures 1D and 6D). Possible mechanisms behind these sustained abnormalities are discussed below.

The late appearance of OXPHOS and ROS pathways correlates with differential immune recovery

In late, severe COVID-19, whole-blood transcriptome analysis showed prominent inflammation-related signatures that were distinct from those seen early in disease. These signatures were related to OXPHOS, ROS, and heme metabolism. These were demonstrated in an unbiased fashion using WGCNA, where modules characterized by OXPHOS and heme metabolism were prominent at days 25 to 48 post symptoms, with OXPHOS most prominent in group E and heme metabolism in C-E (Figures 3B and 3C). Enrichment of Hallmark signatures confirmed the association of OXPHOS and heme metabolism in groups C–E and also found association of a ROS signature (Figures 3E and 7A). Consistent with this, the expression of the genes driving the enrichment of each signature was upregulated in the three most severe clinical groups (Figure 7B; Table S3). The late rise in these three correlated signatures occurred irrespective of persisting or resolving CRP-associated inflammation (Figure 6D) and appeared independent of specific cell population recovery (Figure S7A).

Figure 7

Altered transcriptional changes in prolonged disease (n = 183)

(A) Enrichment score for Hallmark genesets capturing heme metabolism, OXPHOS- and ROS-related genes (by GSEA) in groups A–E in samples taken 25–48 days post screening (group A) or symptom onset (groups B–E).

(B) Heatmap showing relative expression of the intersection of GSEA leading edge genes from groups C, D, and E, across severity groups in samples taken 25–48 days post screening (group A) or symptom onset (groups B–E).

(C) Heatmap showing correlation between transcriptional eigengenes and absolute cell counts, at 25–48 days post symptom onset. Boxes are colored by strength of correlation, Pearson correlation p values: ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

(D) Schematic representation of the trajectory of immunological changes in SARS-CoV-2 infection over time in patients with persisting or resolving systemic inflammation.

See also Figure S7.

We then correlated cellular and transcript signature changes in COVID-19. In the first 24 days after symptom onset, there was a strong association between TNF-α/IL-6, neutrophil degranulation, and interferon signatures with most cell subsets reduced in severe disease (Figure S7B). Later, between 24 and 48 days after symptom onset, these associations changed (Figures 7C and S7C). While TNF-α/IL-6 and neutrophil degranulation signatures were still associated with many cell subsets that continue to be reduced, the interferon signature was no longer a significant player. The persistent increase seen in activated effector CD4+ and CD8+ T cells and plasmablasts was now particularly associated with the OXPHOS signature, which, having become more prominent later in disease (Figure 3B), had a much more restricted and specific association with immune dysfunction than other inflammatory signatures.

It was thus clear that, for some cell types, the association between their number and the inflammatory milieu changed over time but, for others, was more consistent. It was the inflammatory signatures that appeared late in disease, in particular OXPHOS, that were specifically associated with persistent derangement of cell types of potential pathological importance, such as increased HLA-DR+CD38+ T cells and plasmablasts and reduced pDCs (Figure 7D).

Discussion

In asymptomatic or mildly symptomatic, non-progressive SARS-CoV-2 infection (groups A and B), there was evidence of an early robust adaptive immune response. Circulating plasmablasts and CD8+HLA-DR+CD38+ activated T cells expanded earlier and in higher numbers than in more severe COVID-19 groups, most notably in the first week after symptom onset. Both of these cell populations then contracted in A and B, as they continued to rise in groups C–E. Despite this, and a prominent early interferon transcription signature, both the antibody and T cell SARS-CoV-2-specific immune responses were comparable between early mild and severe disease. This suggested that the increased plasmablasts and effector CD8+ T cells could reflect an enhanced bystander response in mild disease, something then supported by single-cell CITE-seq analysis.

Bystander CD8+ T cell activation was first described in the context of lymphocytic choriomeningitis virus (LCMV) infection in mice (Tough et al., 1996). It involves activation of memory CD8+ T cells independent of TCR stimulation (whether by the initiating or cross-reactive antigen) (Maurice et al., 2021). It is driven by type I IFN in many inflammatory contexts—in particular, viral infection. Bystander CD8+ T cells are activated well before the antigen-specific response is seen, being detectable within a day of viral infection and often peaking within a week (Berg et al., 2003). Bystander cells can migrate to infected organs using CXCR3 and kill virally infected cells through NKG2D- and granzyme-B-dependent mechanisms (Maurice et al., 2019). They have been shown to be critically important in the early defense against viral infections, in part by direct anti-viral effects and in part by producing IFN-γ and other cytokines that then activate antigen-presenting cells (Soudja et al., 2014) and control the memory/effector balance of subsequent antigen-specific responses (Krummel et al., 2018). Similar, though less well defined, bystander phenomena occur in CD4+, MAIT, and gd T cell responses (Holzapfel et al., 2014; Ribeiro et al., 2015; Ussher et al., 2018). The early and pronounced production of bystander-activated CD8+ T cells in patients groups A and B could be directly related to disease resolution and prevention of progression to more severe disease, a process likely to begin before symptom onset. Understanding the factors that impede this early bystander response could be important for developing strategies to prevent severe disease.

In contrast, persistent bystander CD8+ T cell activation has been associated with inflammatory pathology in the context of both chronic infection and autoimmunity. NKG2D-dependent killing of non-virally infected hepatocytes in hepatitis A exacerbates liver damage (Kim et al., 2018) and may also play a role in inflammatory lung disease (Borchers et al., 2009), and there is evidence that NKG2D ligands are expressed in the lungs in COVID-19 (Wauters et al., 2021), so it could be that a late and persistent expansion of bystander-activated effector T cells in severe COVID-19 could be a driver of lung pathology. The well-documented links between bystander CD8+ T cell activation and autoimmunity (Groh et al., 2003; Meresse et al., 2004) also raise the possibility that they may play a role in the autoimmune manifestations observed in COVID-19.

The early increase in circulating plasmablasts seen in mild COVID-19 is likely to be comprised of SARS-CoV-2-specific cells, perhaps to an extent not yet fully reflected in serum antibodies, together with non-SARS-CoV-2-specific bystander-activated cells. Bystander B cell activation is known to involve memory (Horns et al., 2020) as well as innate-like B cells (including B1 cells in the mouse and marginal zone cells in both mouse and human). Their function is not fully defined but includes an increase in production of “natural antibodies” (Antin et al., 1986) alongside non-antibody-associated functions, which may include antigen transport to secondary lymphoid organs, antigen presentation to T cells, cytokine production, and immune regulation. These bystander B cell functions are known to play a role in early defense against bacterial infection (Boes et al., 1998; Haas et al., 2005; Ochsenbein et al., 1999) and may be involved in autoimmunity (Espéli et al., 2019; Sanderson et al., 2017). Their role in viral infection is less clear, but it may be that early bystander B cell activation helps determine the outcome of COVID-19 in a manner analogous to the role proposed for bystander CD8+ T cells. A closer investigation of the early B cell response in COVID-19 will be required to confirm this.

At the same time as this pronounced early immune response is seen in groups A and B, there is no evidence of systemic inflammation, apart from early, transient complement activation. CRP, circulating TNF-α and IL-6, and transcriptional signatures of a number of inflammatory pathways are not raised in groups A and B but are already prominent in groups C–E. The severe and widespread leucocyte subset depletion seen at the initial bleed in the more severe COVID-19 patient groups C–E, and observed by many others (Carissimo et al., 2020; Laing et al., 2020; Mathew et al., 2020; Su et al., 2020), is not apparent in those with asymptomatic or mildly symptomatic disease, suggesting that this, too, is a feature of a pathological immune response. Immune cell subset numbers correlated strongly with severe and progressive disease. Coupled with evidence of early systemic inflammation in severity groups C–E, our findings suggest that the immune pathology associated with severe COVID-19 either is established immediately post infection or, if there is a transition point from an effective to a pathological response, is likely to occur before the time of symptom onset. This finding may have major implications as to how disease needs to be managed, because intervention to prevent immune pathology would need to be targeted very early in the disease course and perhaps pre-emptively in high-risk groups screened and diagnosed before symptoms develop.

The reason for the failure to mount a robust early B and T cell response in the context of severe COVID-19 was likely to be multifactorial. There was no evidence for a relationship with viral load and progression to inflammatory disease, as initial viral titers were comparable between groups B, C, and D. Once inflammatory disease is established, however, viral titer may be associated with subsequent outcome, as increased viral titer was seen in group E, consistent with reports that high initial viral titer might be associated with mortality (Pujadas et al., 2020). Genetic association studies in severe COVID-19 point to genes that are implicated in driving antiviral responses. The most prominent are associations with genes involved in the type I interferon pathway (Pairo-Castineira et al., 2020; Zhang et al., 2020), known to be the key driver of bystander T cell activation (Maurice et al., 2021). Increasing age and comorbidity such as diabetes and chronic inflammatory disease are known to suppress early CD8+ T and B cell responses (Shen-Orr et al., 2016; Weiskopf et al., 2009). An in-depth understanding of these risk factors may instruct strategies to assess risk of progression before inflammatory responses become self-sustained.

While clear distinctions in immune responsiveness were apparent between groups A and B versus groups C, D, and E, differences between the more severe groups themselves were less obvious. Those with symptomatic disease warranting admission to hospital clustered together using the size of just 13 key cell populations: these clusters correlated strongly with clinical severity, and immune cell subset numbers together with inflammatory mediator concentrations provided prediction of subsequent progression, as well as COVID-19-associated death. Similar observations have been made by others (Chen and John Wherry, 2020). Predicting progression after presentation to hospital with moderate/severe disease could be of limited clinical use, and it would be of more benefit to predict progression in cases with milder COVID-19, but this may not be possible in practice, given that inflammatory immunopathology is likely to be present at symptom onset. A study to address this issue would need to be conducted in particularly high-risk patient groups to ensure an adequate event rate and would require diagnosis through asymptomatic screening to detect changes before symptoms develop.

We can find no equivalent studies in other severe infections that explore the kinetics of the cellular immune response across this range of time and clinical severity. Nonetheless, many of the changes seen in response to SARS-CoV-2 have been reported in the context of other viral infections. Both an early reduction in blood pDC and T cells has been shown in respiratory syncytial virus (RSV) (Russell et al., 2017) and influenza A (Fox et al., 2012; Lichtner et al., 2011), for example. NK cells in the blood also fall in influenza, probably as a result of migration to the lungs. Finally, plasmablast expansion has also been seen, particularly in dengue fever (Wrammert et al., 2012). The nature of these studies makes a direct comparison with the situation in SARS-CoV-2 difficult, but it would seem likely that many immunological changes seen in COVID-19 mirror those seen in responses to other infections, though, in general, the cellular changes we observe appear to be more profound, peak earlier in disease course, and are more persistent in COVID-19 than in studies performed in other infectious contexts.

The recovery of the profound immune dysregulation seen in severe COVID-19 is potentially of major clinical relevance, as such recovery may be required for the resolution of inflammatory disease or to prevent secondary infection or SARS-CoV-2 reinfection. We found that immune cell abnormalities often persist for weeks to months after SARS-CoV-2 infection, and different cell populations exhibited quite different patterns of resolution. Some recovered as systemic inflammation resolved, others did not, and a third group resolved even in the face of persistent systemic inflammation. Understanding the inflammatory drivers or associations of this differential recovery could provide insight into the immune pathology of COVID-19 and potentially of other infections. To begin to explore this, we correlated immune changes with measurements of systemic inflammation throughout the disease course. Patients with severe COVID are characterized by high CRP, and this correlates with evidence of TNF-α- and IL-6-driven processes at both the protein and transcriptome expression levels, as well as with both neutrophil activation and glycolytic metabolism. The fact that many cellular abnormalities persisted while these biologic processes were apparent, while others appeared to resolve alongside them, suggests that the nature of systemic inflammation is important in driving different aspects of immune pathology.

Finally, some cell populations remained markedly abnormal or showed a limited recovery, even once CRP-associated inflammation had resolved and indeed after patients had been discharged from hospital. These persistent changes may reflect a slow intrinsic regenerative capacity of the cell type concerned, but in other situations, such as the continued elevation of effector T and B cells, it is tempting to speculate that there is ongoing abnormal signaling driving such changes. For that reason, we explored late changes that are seen in the inflammatory response in COVID-19. Transcriptional signatures associated with OXPHOS-, ROS-, and heme-related metabolic pathways arose late in those with severe COVID-19 but were not prominent in early disease. Activation of immune cells results in metabolic reprogramming that supports cell growth, proliferation, and differentiation. Disruption of metabolic pathways can, through many machanisms, result in immune dysfunction (Bantug et al., 2018). It is unlikely that the metabolic signatures observed here simply reflect heightened bioenergetic requirements of activated immune cells, as one would expect similar requirements to be also present early in disease. OXPHOS can drive inflammation (Mills et al., 2017), and we note that COVID-19 patients treated with metformin, which inhibits complex I of the respiratory chain, had lower amounts of circulating inflammatory cytokines (Cheng et al., 2020). The ROS transcriptional signature may reflect more abundant production of ROS, inevitably accompanying increased OXPHOS. Alternatively, it may reflect specific mitochondrial pathology and thus contribute to immune cell dysfunction (Nathan and Cunningham-Bussel, 2013). Mitochondria are also critically involved in heme biosynthesis. Heme serves as a prosthetic group for haemoglobin as well as many other proteins, including several that constitute the respiratory chain of mitochondria. While free heme can act as a damage-associated molecular pattern and promote ROS formation, the role of heme biosynthesis versus catabolism in balancing cellular sensitivity to oxidants is complex and context dependent (Prestes et al., 2020). Here, given correlated regulation of heme and OXPHOS pathways in the clinical groups C–E, activity of these modules may be interrelated and possibly together reflect dysfunctional mitochondria. How heme and OXPHOS transcriptional programmes are linked on a molecular level cannot be inferred from our data. Erythroid cell activation has recently been detected in severe COVID-19 (Bernardes et al., 2020) and could also contribute to a heme transcriptional signature. However, the increase in heme metabolism in our cohort correlates strongly with falling haemoglobin, and reticulocytes in patients in groups C, D, and E are low, suggesting suppression rather than activation of erythropoiesis. Understanding the mechanism linking metabolic dysregulation to persistent immune pathology in COVID-19 will require further study over longer disease courses.

Limitation of study

Many of the abnormalities we have observed in COVID-19 might also be features of other severe viral infections. To identify which are COVID specific will require a comparison with an appropriate disease control group. Sample size is critical when studying heterogeneous disease, and severe COVID-19 falls into this category. While the cohorts described here are comparable in size with other large detailed immunophenotyping studies, larger ones would nonetheless increase our power. In particular, validation of the prognosis prediction models described here is required in large independent cohorts. Very early features of the immune response are associated with disease outcome; exploring these during the period between infection and symptom onset will be important for understanding disease progression but presents major practical difficulties. Continued follow-up of patients will be needed to determine the persistence of abnormalities still observed at late time points. Finally, because our patients were recruited during the first pandemic wave, a follow-up study examining the immune response to new SARS-CoV-2 strains with different virulence could be informative.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact Prof Kenneth Smith (kgcs2@cam.ac.uk).

Materials availability

This study did not generate new unique reagents.

Data and code availability

The datasets generated during this study are available at NIHR CITIID COVID-19 Cohort (https://www.covid19cellatlas.org/patient/citiid/). In addition, whole blood RNA-seq data are available at European Genome-phenome Archive (EGA, ID:EGAS00001005332), and flow cytometry data are available at FLOW Repository (IDs: FR-FCM-Z3XQ, FR-FCM-Z3SR, FR-FCM-Z3ST, FR-FCM-Z3SS).

CITE-seq processed data are available to download from Array Express using accession number E-MTAB-10026.

Experimental model and subject details

Human subjects

Study participants were recruited between 31/3/2020 and 20/7/2020 from patients attending Addenbrooke’s Hospital with a suspected or nucleic acid amplification test (NAAT) confirmed diagnosis of COVID-19 (including point of care testing (Collier et al., 2020; Mlcochova et al., 2020)), patients admitted to Royal Papworth Hospital NHS Foundation Trust or Cambridge and Peterborough Foundation Trust with a confirmed diagnosis of COVID-19, together with Health Care Workers identified through staff screening as PCR positive for SARS-CoV-2 (Rivett et al., 2020). Detailed information on the patients (age, gender, and clinical features) can be found in Table S1. Controls were recruited among hospital staff attending Addenbrooke’s serology screening program, and selected to cover the whole age spectrum of COVID-19 positive study participants, across both genders. Only controls with negative serology results (45 out of 47) were subsequently included in the study. Recruitment of inpatients at Addenbrooke’s Hospital and Health Care Workers was undertaken by the NIHR Cambridge Clinical Research Facility outreach team and the NIHR BioResource research nurse team. Ethical approval was obtained from the East of England – Cambridge Central Research Ethics Committee (“NIHR BioResource” REC ref 17/EE/0025, and “Genetic variation AND Altered Leucocyte Function in health and disease - GANDALF” REC ref 08/H0308/176). All participants provided informed consent.

Inpatients were sampled at study entry, and then at regular intervals as long as they remained admitted to hospital (approximately weekly up to 4 weeks, and then every 2 weeks up to 12 weeks). Discharged patients were invited to provide a follow-up sample 4-8 weeks after study enrolment. Health care workers were sampled at study entry, and subsequently after approximately 2 and 4 weeks. At each time-point, blood samples were drawn in EDTA, sodium citrate, serum and PAXgene Blood RNA tubes (BD Biosciences) and processed by the CITIID-NIHR COVID BioResource Collaboration group (Document S1).

Cell lines

HEK293 T cells (a female cell line) were a kind gift from Paul Lehner, authenticated by STR profiling (Menzies et al., 2018; Miles et al., 2017). They were cultured in DMEM supplemented with 10% fetal calf serum (FCS), 100 units/mL penicillin, and 0.1 mg/ml streptomycin at 37 C in 5% CO2.

Virus

The virus used in this study was the clinical isolate SARS-CoV428 2/human/Liverpool/REMRQ0001/2020, a kind gift from Ian Goodfellow (University of Cambridge), isolated by Lance Turtle (University of Liverpool), and David Matthews and Andrew Davidson (University of Bristol) (Daly et al., 2020; Patterson et al., 2020).

Method details

Clinical data collection

Clinical data were retrospectively collected by review of medical charts and entered into spreadsheets or Castor EDC, a cloud-based clinical data management system. Available laboratory test results and administered in-patient medications were extracted from Epic electronic health records (Addenbrooke’s Hospital) and from MetaVision ICU (RPH ITU). Data were merged from the various data sources using R version 3.6 and the R packages readr (1.3.1), openxlsx (4.1.4), dplyr (0.8.3), tidyr (1.0.2) and lubridate (1.7.4).

Health care workers were classified in 2 groups (A and B) according to whether they were asymptomatic (group A) or had possible COVID-19 symptoms (group B) at the time of PCR testing. For this purpose, new-onset fever (> 37.8 C), cough, loss of sense of smell, hoarseness, nasal discharge or congestion, shortness of breath, wheeze, headache, muscle aches, nausea, vomiting and diarrhea were considered as possible COVID-19 symptoms.

Participants in group A were further sub-grouped according to whether they were completely asymptomatic (n = 8), or had had any of the above COVID-19 symptoms before PCR testing (n = 10, median time from symptoms to COVID-19 PCR test 26 days, range 9-42 days).

Group B participants included both staff who were self-isolating because of COVID-19 symptoms (n = 30), and staff members who were reporting fit for duty, but had some symptoms that did not reach the threshold for self-isolation at that time (n = 10).

Hospital patients were assigned to 3 severity groups, mainly reflecting the maximum intensity of respiratory support for COVID-19 received during their hospital stay:

group C: did not receive any supplemental oxygen. Five patients were discharged soon after initial diagnosis and assessment but followed up as part of the study.

group D: received supplemental oxygen using low flow nasal prongs, simple face mask, Venturi mask or non-rebreather face mask.

group E: received any of non-invasive ventilation (NIV), mechanical ventilation or ECMO. Patients who received supplemental oxygen (but no ventilation) and deceased in hospital were also assigned to group E.

In patients who were already established on home NIV for chronic respiratory failure, NIV delivered as per the home prescription (e.g., nocturnal) was not considered for the purpose of classification. Moreover, oxygen requirements that were clearly not related to COVID-19 were also not considered for classification purposes. In particular, 2 cases who received low flow supplemental oxygen for non-COVID-19 indications (ascitic splinting in decompensated cirrhosis in one case, and recovery from anesthesia after orthopedic surgery in the other) were assigned to group C. Cases in group C were further sub-classified according to chest radiology results (X-ray and, if available, CT scan), as:

abnormal radiology: chest X-ray/ CT scan showed changes compatible with COVID-19

normal radiology: chest X-ray/ CT scan did not show any abnormality compatible with COVID-19 (reported as normal or showing lung changes diagnostic of conditions other than COVID-19).

Immunological parameters were analyzed according to time since onset of symptoms, or otherwise time since positive SARS-CoV-2 NAAT (in group A and in 4 asymptomatic patients in group C). Seven cases admitted to hospital for COVID-19 had no date of onset of symptoms documented in the medical records. In these cases, the date of onset of symptoms was estimated as follows: hospital admission date - median time from symptoms to hospital admission in patients admitted for COVID-19.

Following clinician review, 6 cases were considered not classifiable, due to complex concomitant pathologies that coexisted with COVID-19 and dominated the clinical picture, confounding the interpretation of clinical outcome. These cases were not included in any analyses; more details are reported in Table S2. Data S1 summarizes the timing of research samples and clinical trajectories for volunteers in severity groups C, D and E included in the analysis.

Peripheral blood mononuclear cell preparation and flow immunophenotyping

Each participant provided 27 mL of peripheral venous blood collected into 9 mL sodium citrate tube. Peripheral blood mononuclear cells (PBMCs) were isolated using Leucosep tubes (Greiner Bio-One) with Histopaque 1077 (Sigma) by centrifugation at 800x g for 15 min at room temperature. PBMCs at the interface were collected, rinsed twice with autoMACS running buffer (Miltenyi Biotech) and cryopreserved in FBS with 10% DMSO. All samples were processed within 4 h of collection.

Five distinct antibody cocktails (Data S4) were used to label approximately 106 PBMCs using standard methods. T regulatory cells were fixed and permeabilized following surface staining prior to the addition of intracellular antibodies. Samples were stored at 4°C and acquired within 4 h using a 5-laser BD Symphony X-50 flow cytometer. Single color compensation tubes (BD CompBeads) or cells were prepared for each of the fluorophores used and acquired at the start of each flow cytometer run.

For direct enumeration of T, B and NK cells, an aliquot of whole blood (50 μl) was added to BD TruCount tubes with 20μl- BD Multitest 6-color TBNK reagent (BD Biosciences) and processed as per the manufacturer’s instructions.

Samples were gated in FlowJo v10.2 according to the schema set out in Data S4. The number of cells falling within each gate was recorded. For analysis, these were expressed as an absolute concentration of cells per ml, calculated using the proportions of daughter populations present within the parent population determined using the BD TruCountsystem.

CyTOF

The protocol used to isolate PBMCs led to an impaired recovery of the different monocyte populations, specifically intermediate and non-classical monocytes (Figure S2). To extend our analysis to these and other granulocyte populations we performed mass cytometric analysis on a subgroup of patients and healthy controls (249 samples). Briefly, whole blood samples (270μl) were stained using the Fluidigm Maxpar® Direct Immune Profiling Assay according to the manufacturer’s instructions. Samples were cryopreserved at −80°C following staining and thawed for analysis within 4 weeks. Samples were acquired using a Fluidigm Helios mass cytometer and normalized using the CyTOF Software v6.7.1016. FCS files generated were analyzed using the Maxpar® Pathsetter software v2.0.45 (Verity Software House, Topsham, ME). Standard settings were used to generate immune cell frequencies for 37 immune cell populations. Absolute cell numbers were calculated using the proportions of these immune cell populations within the parent populations determined by BD TruCount.

Reticulocyte counts

Reticulocyte numbers were measured using a Sysmex XN-1000 hematology analyzer according to manufacturer instruction and as previously described (Akbari et al., 2020). Briefly, Sysmex technology uses three signals to define the physiological and structural properties of cells and to distinguish reticulocytes from the other blood cells: forward scatter, side scatter and side fluorescent light. These measurements rely on the similar electromagnetic radiation and fluid dynamics concepts of a flow cytometer; reticulocyte specific fluorescent probes are covered by a patent deposited by Sysmex Corporation (i.e., Fluorocell RET, cat# BN-337-547).

Reticulocyte Ratio:

Reticulocyte Count:

Complement

Complement activation was assessed by measuring C3 activation products (C3a and C3c) together with the terminal complement complex (TCC) as an end product of the complement cascade. Concentrations of these complement components were measured in EDTA plasma from patients using commercially available enzyme-linked immunosorbent assays (ELISA) kits (HK354 (C3a), HK368 (C3c), HK328 (TCC), Hycult Biotech, Uden, the Netherlands) according to the manufacturer’s protocols.

CRP

High sensitivity CRP was measured using the standard assay by the Core Biochemical Assay Laboratory (CBAL) at Cambridge University Hospitals NHS Foundation Trust.

Cytokines

IL-6, IL-10, IL-1β, TNF-α and IFN-γ were measured in serum from patients and HCs by high sensitivity Base Kit HS Cytokine A Mag (cat# LHSCM000, R&D Systems / Biotechne) on a Luminex analyzer (Bio-Plex, Bio-Rad, UK) as standard clinical assay performed by the Clinical Immunology Laboratory at the Department of Biochemistry and Immunology, Addenbrooke’s Hospital Cambridge.

IFN-μ FLUOROSPOT assays

Frozen PBMCs were rapidly thawed, and the freezing medium was diluted into 10ml of TexMACS media (Miltenyi Biotech), centrifuged and resuspended in 10ml of fresh media with 10U /ml DNase (Benzonase, Merck-Millipore via Sigma-Aldrich), PBMCs were incubated at 37°C for 1 h, followed by centrifugation and resuspension in fresh media supplemented with 5% Human AB serum (Sigma Aldrich) before being counted. PBMCs were stained with 2μl of each antibody: anti-CD3- fluorescein isothiocyanate (FITC), clone UCHT1; anti-CD4- phycoerythrin (PE), clone RPA-T4; anti-CD8a- peridinin-chlorophyll protein - cyanine 5.5 (PerCP Cy5.5), clone RPA-8a (all BioLegend, London, UK), LIVE/DEAD Fixable Far Red Dead Cell Stain Kit (Thermo Fisher Scientific). PBMC phenotyping was performed on the BD Accuri C6 flow cytometer. Data were analyzed with FlowJo v10 (Becton Dickinson, Wokingham, UK). 1.5 to 2.5 × 105 PBMCs were incubated in pre-coated Fluorospot plates (Human IFN-γ FLUOROSPOT (Mabtech AB, Nacka Strand, Sweden)) in triplicate with peptide mixes specific for Spike, Nucleocapsid and Membrane proteins of SARS-CoV-2 (final peptide concentration 1μg/ml/peptide, Miltenyi Biotech) and an unstimulated and positive control mix (containing anti-CD3 (Mabtech AB), Staphylococcus Enterotoxin B (SEB), Phytohaemagglutinin (PHA) (all Sigma Aldrich)) at 37°C in a humidified CO2 atmosphere for 48 h. The cells and medium were decanted from the plate and the assay developed following the manufacturer’s instructions. Developed plates were read using an AID iSpot reader (Oxford Biosystems, Oxford, UK) and counted using AID EliSpot v7 software (Autoimmun Diagnostika GmbH, Strasberg, Germany). All data were then corrected for background cytokine production and expressed as SFU/Million CD3+T cells.

SARS-CoV-2 serology

Quantification of Spike SARS-CoV-2 specific antibodies was performed by ELISA as described by Xiong X et al. (Xiong et al., 2020). Briefly, serum samples collected at time of enrolment in the study and at the 4-8 week follow-up visit were first screened for positivity and then antibody titers were determined by an end-point analysis. AUC values were calculated in R (3.6.3) using the flux (0.3-0) package. Kruskal–Wallis test was used to calculate p values among the different disease severities.

SARS-CoV-2 neutralisation assay

The virus used in this study was the clinical isolate SARS-CoV-2/human/Liverpool/REMRQ0001/2020, a kind gift from Ian Goodfellow (University of Cambridge), isolated by Lance Turtle (University of Liverpool) and David Matthews and Andrew Davidson (University of Bristol) (Daly et al., 2020; Patterson et al., 2020).

Sera were heat-inactivated at 56°C for 30 min, then frozen in aliquots at −80°C. Neutralising antibody titers at 50% inhibition (NT50s) were measured as previously described (Gerber et al., 2021).

In brief, HEK293T reporter cells expressing Renilla luciferase (Rluc) and SARS-CoV-2 Papain-like protease-activatable circularly permuted firefly luciferase (FFluc) were seeded in flat-bottomed 96-well plates. The next day, SARS-CoV-2 viral stock (MOI = 1) was pre-incubated with a 3-fold dilution series of each serum for 2 h at 37°C, then added to the cells. After 24 h, cells were lysed in Dual-Glo Luciferase Buffer (Promega) diluted 1:1 with PBS and 1% NP-40. Lysates were transferred to white half-area 96-well plates, and infectious virus quantitated as the ratio of FFluc/Rluc activity measured using the Dual-Glo kit (Promega) according to the manufacturer’s instructions.

Experiments were conducted in duplicate. To obtain NT50s, FFluc/Rluc ratios were analyzed using the Sigmoidal, 4PL, X is log(concentration) function in GraphPad Prism.

Whole blood bulk RNA-Seq

Whole blood RNA was extracted from PAXgene Blood RNA tubes (BD Biosciences) of 188 COVID-19 patients, at up to 2 time points, and 42 healthy volunteers. RNA-Sequencing libraries were generated using the SMARTer® Stranded Total RNA-Seq v2 - Pico Input Mammalian kit (Takara) using 10ng RNA as input following the manufacturer’s protocol. Libraries were pooled together (n = 96) and sequenced using 75bp paired-end chemistry across 4 lanes of a Hiseq4000 instrument (Illumina) to achieve 10 million reads per sample. Read quality was assessed using FastQC v.0.11.8 (Babraham Bioinformatics, UK), and SMARTer adaptors trimmed and residual rRNA reads depleted in silico using Trim_galore v.0.6.4 (Babraham Bioinformatics, UK) and BBSplit (BBMap v.38.67(BBMap - Bushnell B. - https://sourceforge.net/projects/bbmap/)), respectively. Alignment was performed using HISAT2 v.2.1.0 (Kim et al., 2019) against the GRCh38 genome, achieving a greater than 95% alignment rate. Count matrices were generated using featureCounts (Rsubreads package - (Liao et al., 2020) and stored as a DGEList object (EdgeR package (Robinson et al., 2010) for further analysis.

All downstream data handling was performed in R (R Core Team, 2015). Counts were filtered using filterByExpr (EdgeR) with a gene count threshold of 10 CPM and the minimum number of samples set at the size of the smallest disease group. Library counts were normalized using calcNormFactors (EdgeR) using the method ‘weighted trimmed mean of M-values’. The function ‘voom’ (Law et al., 2014) was applied to the data to estimate the mean-variance relationship, allowing adjustment for heteroscedasticity.

Single cell RNA-seq

CITE-seq data were generated from frozen PBMCs as described by Stephenson et al. (Stephenson et al., 2021). Briefly, after thawing, pools of 4 samples were generated by combined 500,000 viable cells per individual (total of 2 million cells per pool). TotalSeq-C antibody cocktail (BioLegend 99813) was used to perform cell surface marker staining on 500,000 cells per pool. 50,000 live cells (up to a maximum of 60,000 total cells) for each pool were processed using Single Cell V(D)J 5¢ version 1.1 (1000020) together with Single Cell 5¢ Feature Barcode library kit (1000080), Single Cell V(D)J Enrichment Kit, Human B Cells (1000016) and Single Cell V(D)J Enrichment Kit, Human T Cells (1000005) (10xGenomics) according to the manufacturer’s protocols. Samples were sequenced on NovaSeq 6000 (Illumina) using S1 flowcells. Droplet libraries were processed using Cellranger v4.0. Reads were aligned to the GRCh38 human genome concatenated to the SARS-Cov-2 genome (NCBI SARS-CoV-2 isolate Wuhan-Hu-1) using STAR (Dobin et al., 2013) and unique molecular identifiers (UMIs) deduplicated. CITE-seq UMIs were counted for GEX and ADT libraries simultaneously to generate feature X droplet UMI count matrices.

Quantification and statistical analysis

All statistical analyses were conducted using custom scripts in R (R Core Team, 2015). Absolute cell counts (cells/uL) were offset by +1 to allow subsequent log2 transformation of zero counts. Where shown, time measures represent time from symptom onset (for severity groups B, C, D and E) or first positive COVID-19 swab (group A). Unless otherwise specified, longitudinally collected data was grouped by bins of 7 or 12 days. Pairwise statistical comparison of absolute cell counts, CRP or serum measures between individuals in a given severity group at a given time bin and HCs, or between severity groups, was conducted by Wilcoxon test unless otherwise specified. For analyses involving repeated-measures, false discovery rate corrected (Benjamini & Hochberg) p value were reported. For individuals sampled more than once within a given time bin, data from the earliest blood collection was used.

Cell subset deconvolution of the whole blood RNA-Seq dataset was performed using pathway-level information extractor (PLIER) (http://gobie.csb.pitt.edu/PLIER). Latent factors were generated by leveraging off pre-existing knowledge of cell specific pathways. To better understand the relationship between gene expression and clinical severity, weighted gene co-expression network analysis was carried out using the WGCNA package (Langfelder and Horvath, 2008) in R. Briefly, a signed adjacency matrix was generated and a soft thresholding power was chosen to impose approximate scale-free topology. Modules were identified from the resulting topological overlap matrix with a specified minimum module size of n = 30. Modules were summarized using singular value decomposition, and the resulting module eigengene correlated with clinical traits. Significance of the correlation between a given clinical trait and a modular eigengene was assessed using linear regression with Bonferroni adjustment to correct for multiple testing. Modules were annotated using Enrichr (Chen et al., 2013). Longitudinal mixed modeling of gene module changes over time was conducted using the nlme package in R (Pinheiro et al., 2021), including time () with a quadratic trend and disease severity category or unsupervised cluster ids () as fixed effects, and sampled individuals as random effects ():

i.e., using the lme formula:

Gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was used to identify biological pathways enriched in COVID-19 severity groups relative to healthy controls. Briefly, a list of ranked genes, determined by Signal-To-Noise ratio, was generated. An enrichment score was calculated, determined by how often genes from the geneset of interest appeared at the top or the bottom of the pre-ranked set of genes with the enrichment score representing the maximum deviation from zero. To assess statistical significance, an empirical phenotype- based permutation test was run, where a collection of enrichment scores was generated from the random assignment of phenotype to samples and used to generate a null distribution. To account for multiple testing, an FDR rate q < 0.20 was deemed significant. A leading edge analysis was performed to determine the genes contributing the most to the enrichment of a given pathway and was subsequently illustrated in a heatmap. Hallmark gene sets from the Molecular Signatures Database (https://www.broadinstitute.org/gsea/msigdb) were used in analysis.

Principal component analysis (PCA) of centered and scaled absolute counts for 24 major cell types was conducted using the pca()function from the package mixOmics (Rohart et al., 2017). Unsupervised clustering of log2 transformed absolute cell counts, normalized to the median of the corresponding control population, was conducted using the heatmap.2() function from the package gplots (Warnes et al., 2020), with a Euclidean distance function applied to both rows and columns of the data matrix and hierarchical clustering computed using the ward.D method. Partial least-squares discriminant analysis (PLS-DA) was conducted using the plsda() function from the package mixOmics (Rohart et al., 2017), a supervised method of sample discrimination whereby sample clustering is informed by group membership (here patient clusters 1 and 2). The classification performance of the PLS-DA model was determined using the perf() function via 10 iterations of 5-fold cross-validation, with two components deemed sufficient to minimize the balanced error rate of prediction. Variable selection on components 1 and 2 was conducted using the tune() function, with 13 cell types selected as those most strongly contributing to discrimination of patient clusters. An AUROC curve showing the performance of a predictive model based on these 13 cell types was generated using the auroc() function. To assess whether clinical severity was reflected on a transcriptional level in an unsupervised fashion, K-means clustering was utilized on normalized whole blood RNASeq gene expression counts. Heatmaps were created using the ComplexHeatmap package (Gu et al., 2016), with data scaled and centered prior to visualization.

Cellular recovery rates over 60 days were calculated for each cell type in patients from groups C, D and E, split into those with persistently elevated (> 10mg/L) or resolving CRP (falling below 10mg/L by final bleed), over 60 and 40 days respectively. Using a 12 day sliding window with single day increments, the ‘window of recovery’ for each given cell population was defined as the window in which absolute cell counts for COVID-19 samples no longer differed from controls when assessed by Wilcoxon test, and remained as such for the subsequent 7 windows, and 80% of all windows remaining. Recovery rate was taken as the log2 normalized ratio of test and control absolute counts for patient samples collected within the first time window (0-12 days), subtracted from the equivalent value calculated within the window of recovery, divided by the upper day boundary of the recovery window.

The relationships between immunophenotyping and transcriptional data in the form of gene expression modules were assessed using Pearson’s correlation (Hmisc package) and visualized with corrplot.

Consortia

The members of the Cambridge Institute of Therapeutic Immunology and Infectious Disease-National Institute of Health Research (CITIID-NIHR) COVID BioResource Collaboration are John Allison, Ali Ansaripour, Stephen Baker, Laura Bergamaschi, Ariana Betancourt, Sze-How Bong, Georgie Bower, John R. Bradley, Ashlea Bucke, Ben Bullman, Katherine Bunclark, Helen Butcher, Jo Calder, Laura Canna, Daniela Caputo, Debbie Clapham-Riley, Chiara Cossetti, Jerome D. Coudert, Eckart M.D.D. De Bie, Aloka De Sa, Eleanor Dewhurst, Giovanni di Stefano, Jason Domingo, Gordon Dougan, Benjamin J. Dunmore, Anne Elmer, Madeline Epping, Codie Fahey, Stuart Fawke, Stewart Fuller, Anita Furlong, Nick Gleadall, Ian G Goodfellow, Stefan Gräf, Barbara Graves, Jennifer Gray, Richard Grenfell, Ravindra K. Gupta, Julie Harris, Christoph Hess, Sarah Hewitt, Andrew Hinch, Josh Hodgson, Elaine Holmes, Christopher Huang, Oisín Huhn, Kelvin Hunter, Tasmin Ivers, Sarah Jackson, Isobel Jarvis, Emma Jones, Jane Kennet, Sherly Jose, Masa Josipovic, Mary Kasanicki, Nathalie Kingston,Jenny Kourampa, Elisa Laurenti, Ekaterina Legchenko, Paul J. Lehner, Emma Le Gresley, Daniel Lewis, Rachel Linger, Paul A. Lyons, Michael Mackay, John C. Marioni, Jimmy Marsden, Jennifer Martin, Cecilia Matara, Nicholas J. Matheson, Anne Meadows, Sarah Meloy, Nicole Mende, Federica Mescia, Alice Michael, Rachel Michel,Lucy Mwaura, Francesca Muldoon, Francesca Nice, Criona O’Brien, Ciara O’Donnell, Georgina Okecha, Ommar Omarjee, Nigel Ovington, Willem H. Owehand, Sofia Papadia, Caroline Patterson, Marianne Perera, Isabel Phelan, Linda Pointon, Petra Polgarova, Gary Polwarth, Nicole Pond, Jane Price, Cherry Publico, Rebecca Rastall, Carla Ribeiro, Nathan Richoz, Veronika Romashova, Sabrina Rossi, Jane Rowlands, Valentina Ruffolo, Caroline Saunders, Natalia Savinykh Yarkoni, Rahul Sharma, Joy Shih, Mayurun Selvan, Kenneth G.C. Smith, Sarah Spencer, Luca Stefanucci, Hannah Stark, Jonathan Stephens, Kathleen E Stirrups, Mateusz Strezlecki, Charlotte Summers, Rachel Sutcliffe, James E.D. Thaventhiran,Tobias Tilly, Zhen Tong, Hugo Tordesillas, Carmen Treacy, Mark Toshner, Paul Townsend, Lori Turner, Neil Walker, Jennifer Webster, Michael P. Weekes, Nicola K. Wilson, Jennifer Wood, Marta Wylot, and Cissy Yong.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Anti-human CCR4 PEVio770 (REA279)	Miltenyi	RRID:AB_2655909
Anti-human CCR5 AF647 (HEK/1/85a)	Biolegend	RRID:AB_528760
Anti-human CCR6 PECy7 (G034E3)	BioLegend	RRID:AB_10916518
Anti-human CCR7 BV650 (G043H7)	BioLegend	RRID:AB_2563867
Anti-human CCR7 APC-fire (G043H7)	BioLegend	RRID:AB_2750147
Anti-human CD11c AF700 (B-ly6)	BD	RRID:AB_10612006
Anti-human CD123 BV786 (7G3)	BD	RRID:AB_2738662
Anti-human CD127 PECy5 (eBioRDR5)	Thermo	RRID:AB_2043801
Anti-human CD127 PerCP efluor710 (eBioRDR5)	Thermo	RRID:AB_2762464
Anti-human CD14 BV711 (MjP9)	BD	RRID:AB_2744290
Anti-human CD14 BV510 (63D3)	BioLegend	RRID:AB_2716229
Anti-human CD141 BV605 (1A4)	BD	RRID:AB_2740151
Anti-human CD15 BV510 (W6D3)	BioLegend	RRID:AB_2563400
Anti-human CD16 BUV496 (3G8)	BD	RRID:AB_2870224
Anti-human CD161 PE (HP-3G10)	BioLegend	RRID:AB_1501083
Anti-human CD163 PE-CF594 (GHI/61)	BD	RRID:AB_2737711
Anti-human CD19 BUV496 (SJ25C1)	BD	RRID:AB_2870221
Anti-human CD193 BV510 (5E8)	BioLegend	RRID:AB_2571977
Anti-human CD1c AF647 (F10/21A3)	BD	RRID:AB_2744318
Anti-human CD20 BUV805 (2H7)	BD	RRID:AB_2870192
Anti-human CD24 BB700 (ML5)	BD	RRID:AB_2744333
Anti-human CD25 PE (BC96)	Thermo	RRID:AB_1659682
Anti-human CD27 BV711 (O323)	BioLegend	RRID:AB_2563809
Anti-human CD28 BV785 (CD28.2)	BioLegend	RRID:AB_2632607
Anti-human CD28 BUV805 (L293)	BD	RRID:AB_2872889
Anti-human CD3 BUV395 (SK7)	BD	RRID:AB_2744382
Anti-human CD3 BUV661 (UCHT1)	BD	RRID:AB_2870239
Anti-human CD3 BV510 (UCHT1)	BioLegend	RRID:AB_2563468
Anti-human CD303 PE-Vio770 (REA693)	Miltenyi	RRID:AB_2657198
Anti-human CD304 PE (U21-1283)	BD	RRID:AB_2744361
Anti-human CD32 BB700 (FLI8.26)	BD	RRID:AB_2871430
Anti-human CD38 BUV661 (HIT2)	BD	RRID:AB_2870242
Anti-human CD39 BV421 (A1)	BioLegend	RRID:AB_2564575
Anti-human CD39 APC-fire (A1)	BioLegend	RRID:AB_2650839
Anti-human CD4 BUV496 (SK3)	BD	RRID:AB_2870220
Anti-human CD40 BUV395 (5C3)	BD	RRID:AB_2739110
Anti-human CD45 BUV805 (HI30)	BD	RRID:AB_2870179
Anti-human CD45RA BV570 (HI100)	BioLegend	RRID:AB_2563813
Anti-human CD45RA BUV805 (HI100)	BD	RRID:AB_2871317
Anti-human CD56 FITC (MEM188)	Thermo	RRID:AB_10372519
Anti-human CD69 BV650 (FN50)	BioLegend	RRID:AB_2563158
Anti-human CD71 BV650 (CY1G4)	BioLegend	RRID:AB_2687103
Anti-human CD73 Brilliant Violet 785™ (AD2)	BioLegend	RRID:AB_2687234
Anti-human CD73 BV785 (AD2)	BioLegend	RRID:AB_2687234
Anti-human CD80 PECy5 (L307.4)	BD	RRID:AB_397239
Anti-human CD86 BUV737 (FUN-1)	BD	RRID:AB_2814790
Anti-human CD86 PECy7 (BU63)	BioLegend	RRID:AB_2728392
Anti-human CD8b BV480 (2ST8.5H7)	BD	RRID:AB_2743598
Anti-human CD95 BUV737 (DX2)	BD	RRID:AB_2870117
Anti-human CRTh2 PE-dazzle (BM16)	BioLegend	RRID:AB_2572053
Anti-human CXCR5 APC-R700 (RF8B2)	BD	RRID:AB_2739103
Anti-human FoxP3 APC (236A/E7)	Thermo	RRID:AB_10804651
Anti-human GLUT1 AF647 (202915)	BD	RRID:AB_2869787
Anti-human Helios Pedazzle (22F6)	BioLegend	RRID:AB_2565797
Anti-human HLA-DR APC-H7 (G46-6)	BD	RRID:AB_10611876
Anti-human HLADR BV605 (L243)	BioLegend	RRID:AB_2561913
Anti-human KLRG1 FITC (REA261)	Miltenyi	RRID:AB_2652570
Anti-human PD-1 BV421 (EH12.2H7)	BioLegend	RRID:AB_10960742
Anti-human TCR Vg9 AF700 (B3)	Biolegend	RRID:AB_2814207
Anti-human TCR-DV1 PECy7 (TS8.2)	Thermo	RRID:AB_2762454
Anti-human TCR-DV2 PerCPCy5.5 (B6)	BioLegend	RRID:AB_2687330
Anti-human TCRgd BUV737 (11F2)	BD	RRID:AB_2872944
Anti-human TCRV7.2 BV711 (3C10)	BioLegend	RRID:AB_2629680
Anti-human Vb11 APCVio770 (REA559)	Miltenyi	RRID:AB_2653747
Zombie Aqua	BioLegend	Cat#423101
Zombie Yellow	BioLegend	Cat#423103
Anti-human CD3 FITC (UCHT1)	BioLegend	RRID:AB_2562046
Anti-human CD4 PE (RPA-T4)	BioLegend	RRID:AB_2562053
Anti-human CD8a PerCP-Cy5.5 (RPA-8a)	BioLegend	RRID:AB_1575074
Fixable Far Red Dead Cell Stain Kit	Thermo	Cat#423109
Chemicals, peptides, and recombinant proteins
Spike SARS-CoV-2 peptide prot-S	Miltenyi Biotec	Cat#130-126-701
Spike SARS-CoV-2 peptide Prot-S1	Miltenyi Biotec	Cat#130-127-048
Nucleocapsid SARS-CoV-2 peptide	Miltenyi Biotec	Cat#130-126-699
Membrane SARS-CoV-2 peptide	Miltenyi Biotec	Cat#130-126-703
Staphylococcus Enterotoxin B (SEB)	Sigma Aldrich	Cat#S4881-1MG
Phytohaemagglutinin (PHA)	Sigma Aldrich	Cat#L1668-5MG
anti-CD3	Mabtech AB	Cat#3605-1-1000
Spike SARS-CoV-2 protein	Xiong et al., 2020	N/A
SARS-CoV-2 neutralisation assay	Gerber, et al., 2021	N/A
Critical commercial assays
6-color TBNK Reagent with BD Trucount™ Tubes	BD	RRID:AB_2868707
Maxpar® Direct™ Immune Profiling Assay™	Fluidigm	Cat#201325
Human IFNg FLUOROSPOT	Mabtech AB	Cat#X-01A-10
hIFN-g HS LxPA MAG	R&D systems/ Biotechne	Cat#LHSCM285B
Human IL-1b Mag Bead Set	R&D systems/ Biotechne	Cat#LHSCM201
Human IL-6 Mag Bead Set	R&D systems/ Biotechne	Cat#LHSCM206
Human IL-10 Mag Bead Set	R&D systems/ Biotechne	Cat#LHSCM217
Human TNF-a Mag Bead Set	R&D systems/ Biotechne	Cat#LHSCM210
Base Kit, HS Cytokine A, Mag	R&D systems/ Biotechne	Cat#LHSCM000
Fluorocell™ RET	Sysmex Corporation	Cat#BN-337-547
C3a, Human, ELISA kit	Hycult Biotech	Cat#HK354
C3c, Human, ELISA kit	Hycult Biotech	Cat#HK368
TCC, Human, ELISA kit	Hycult Biotech	Cat#HK328
SMARTer® Stranded Total RNA-Seq v2 - Pico Input Mammalian kit	Takara	Cat#634414
Deposited data
CITE-seq data	Stephenson et al., 2021, Cambridge cohort	Array Express:E-MTAB-10026
Absolute cell count, Whole blood RNA-seq, cytokine and complement measurements, SARS-CoV-2 specific antibody titers	This paper	https://www.covid19cellatlas.org/patient/citiid/
Flow cytometry data: B cell and Treg panels	This paper	Flow Repository:FR-FCM-Z3XQ
Flow cytometry data: Monocyte panel	This paper	Flow Repository:FR-FCM-Z3SR
Flow cytometry data: Conventional T cell panel	This paper	Flow Repository:FR-FCM-Z3ST;
Flow cytometry data: non-conventional T cell panel	This paper	Flow Repository:FR-FCM-Z3SS
Whole blood RNA-seq data	This paper	EGA:EGAS00001005332
Software and algorithms
Maxpar® Pathsetter™ software v2.0.45	Verity Software House, Topsham, ME	N/A
FlowJo v10.2	FlowJo LLC	N/A
CyTOF Software v6.7.1016	Fluidigm	N/A
AID EliSpot v7 software	Autoimmun Diagnostika GmbH, Strasberg, Germany	N/A
FastQC v.0.11.8	Babraham Bioinformatics, UK	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Trim_galore v.0.6.4	Babraham Bioinformatics, UK	https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/
BBMap v.38.67	BBMap - Bushnell B.	https://sourceforge.net/projects/bbmap/
HISAT2 v.2.1.0	Kim et al., 2019	http://daehwankimlab.github.io/hisat2/
R	R Core Team, 2015	N/A
pathway-level information extractor	PLIER	http://gobie.csb.pitt.edu/PLIER
Gene set enrichment analysis (GSEA)	Broad Institute	https://www.gsea-msigdb.org/gsea/index.jsp
Other
XN-1000 hematology analyzer	Sysmex	N/A
Symphony X-50	BD	N/A
Helios mass cytometer	Fluidigm	N/A
AID iSpot reader	Oxford Biosystems	N/A
Luminex	Bio-Plex, Bio-Rad, UK	N/A