SARS-CoV-2 sculpts the immune system to induce sustained virus-specific naïve-like and memory B cell responses

SARS-CoV-2 triggers a broad antibody response in the early stage of infection and a robust and long-lasting RBD-specific IgG1 response

To dissect the dynamics of humoral immune responses to SARS-CoV-2, we collected blood and serum samples from a cohort of 25 hospitalized COVID-19 patients in the acute phase of infection (COVT1). Of these patients, 24% (n = 6; darker symbols in Figures) were admitted to intensive care unit (ICU). On average, COVT1 samples were collected 11 days post-symptom onset (PSO). A follow-up blood sample was drawn from 20 of 25 patients (COVT2), including 4 ICU patients, averaging 70 days PSO (Fig. 1A). Blood and serum samples were also collected from healthy controls (n = 21). Aside from being a major risk factor for fatal outcome in COVID-19 patients, age profoundly influences the immune system of healthy individuals. To minimize these age-related effects, we mostly selected patients younger than 65 (median age = 51) and included age-matched healthy controls (median age = 50). A summary of demographic and clinical data of all the individuals included in this study is provided (Table S1).

• Download figure
• Open in new tab

Fig. 1. SARS-CoV-2 induces a broad antibody response in the early stage of infection and a long-lasting IgG response.

(A) Schematic diagram of the study timeline and cohort characteristics. Range of days post-symptom onset (PSO) is indicated with a box and median is indicated with a line for COVID-19 patients in the acute (COVT1) and convalescent (COVT2) phase of infection. (B) Area under the curve (AUC) for each of the RBD-specific and (C) NP-specific antibody classes and subclasses analyzed from COVT1 and COVT2 sera samples. Sera from healthy controls (HCs) were analyzed in parallel to establish negative threshold values defined as the HC AUC mean plus 2 times the standard deviation of the mean. Data are presented as paired individual dots. Dashed line indicates negative threshold. Wilcoxon matched pairs test (*P < 0.05, **P < 0.01, and ***P < 0.001). (D) Relative change in antibody levels between COVT1 and COVT2 plotted against the corresponding antibody levels at COVT1. r stands for Spearman’s rank-order correlation. (A-D) Dark-colored dots show ICU patients. COVT1, n=20; COVT2, n=20.

To elucidate whether COVID-19 patients can mount a broad and long-lasting humoral response to SARS-CoV-2, we performed isotype-specific (IgM, IgA and IgG) and subtype-specific (IgA1, IgA2, IgG1, IgG2, IgG3 and IgG4) enzyme-linked immunosorbent assays (ELISAs) that quantified antibody responses to the receptor-binding domain (RBD) of the spike protein or nucleocapsid protein (NP) from SARS-CoV-2. Compared to healthy controls, sera from COVID-19 patients in acute infection (COVT1) and convalescence (COVT2) showed significantly higher titers of all RBD- and NP-specific antibody classes and subclasses analyzed (Fig. S1). According to our longitudinal study and as previously reported^18,27, RBD- and NP-specific IgM and IgA titers significantly decreased two months PSO (Fig. 1B, C), presumably because they mainly originated from short-lived plasmablasts (PBs). Interestingly, both IgA1 and IgA2 subclasses contributed to the decline of virus-specific IgA, suggesting a limited involvement of mucosal immunity during convalescence. Unlike RBD- and NP-specific IgM and IgA, virus-specific IgG titers were maintained over the first two months PSO (Fig. 1B, C). Remarkably, RBD-specific IgG1 was the only IgG subclass that significantly increased overtime (Fig. 1B).

We further studied the dynamics of the humoral response in COVID-19 throughout the study’s time course and found that the magnitude of relative changes in virus-specific IgM or IgG titers directly correlated with their initial levels at COVT1. Indeed, patients with initial higher IgM or lower IgG titers showed greater relative decrease or increase, respectively (Fig. 1D), similarly to what has been previously reported for convalescent individuals at later time points¹⁸. Thus, SARS-CoV-2 infection triggers a broad antibody response in terms of antigen-specificity and antibody isotypes in early stages and an expansion of an RBD-specific and IgG1-dominated response during convalescence.

COVID-19 is associated with transitory expansion of circulating plasmablasts and naïve B cells during the acute phase of infection

To further explore the dynamics of the B cell response to SARS-CoV-2, we analyzed circulating B cells from healthy controls and COVID-19 patients in the acute and convalescent phase of infection by high-dimensional spectral flow cytometry (Table S2). As reported by others²⁸, the frequency of CD19⁺ B cells within live peripheral blood mononuclear cells (PBMCs) was significantly higher in COVT1 compared to healthy controls or COVT2 (Fig. 2A). T-distributed stochastic neighbour embedding (tSNE) defined the major B cell populations (Fig. 2B), the contribution of each group of samples to these populations (Fig. 2C), and the differential expression of multiple surface B cell proteins within these populations (Fig. 2D). In parallel, we queried the data by traditional gating (Fig. S2A, C). This analysis revealed that, for most COVID-19 patients, the acute response to SARS-CoV-2 was dominated by a transitory expansion of CD38⁺⁺CD10−CD27^high plasma cells (PCs; Fig. 2B-F), the majority of which expressed HLA-DR⁺ (Fig. 2G), an antigen-presenting protein typically expressed by newly generated short-lived PBs²⁹. Next, the analysis of antibody classes and subclasses confirmed the induction of unswitched as well as IgA1 or IgG class-switched PCs soon after infection (Fig. 2H and Fig. S2B). We then assessed the homing potential of these PCs through the analysis of CXCR3 and CXCR4 chemokine receptors, which guide PCs to inflamed tissues or bone marrow, respectively³⁰. We found that the majority of PCs from COVT1 had a CXCR3⁺CXCR4^- phenotype, suggesting their targeted migration into inflamed tissues (Fig. 2I, J). Compared to healthy controls and COVT2, COVT1 samples also showed a significantly increased proportion of CXCR3⁺CXCR4⁺ PCs (Fig. 2I, K).

• Download figure
• Open in new tab

Fig. 2. Deep profiling of B cell subsets in COVID-19 patients reveals transitory expansion of circulating plasmablasts, immature transitional and naïve B cells during the acute phase of infection.

(A) Frequency of CD19⁺ B cells from live peripheral blood mononuclear cells (PBMCs) in healthy controls (HCs) and COVID-19 patients at T1 (COVT1) and at T2 (COVT2). (B) Merged tSNE projection of CD19⁺ B cells for HCs (n=11), COVT1 (n=16) and COVT2 (n=16) samples concatenated and overlaid, with main B cell populations indicated by color. PCs, plasma cells. ME, memory. DN2, CD38^dullCD10⁻IgD^-CD27^-CD21^-CD11c⁺ double negative type 2. (C) Merged tSNE projection of CD19⁺ B cells with each sample group indicated by color. (D) tSNE projections of expression of the indicated cell surface markers. (E) Frequency of CD38⁺⁺CD10⁻CD27⁺ PCs from total CD19⁺ B cells and (F) representative dot plots in HCs, COVT1 and COVT2. Numbers indicate percentages in the drawn gates. (G) Frequency of HLA-DR⁺ plasmablasts (PBs) from total circulating CD38⁺⁺CD10⁻CD27⁺ cells in HCs, COVT1 and COVT2. (H) Frequencies of IgM⁺ PCs (left) and IgG/E⁺ PCs (right) from total CD19⁺ B cells. (I) Representative flow cytometry plots (left) and relative percentage (right) of CXCR3^-CXCR4^- (R3^-R4^-), CXCR3^-CXCR4⁺ (R3^-R4⁺), CXCR3⁺CXCR4⁺ (R3⁺R4⁺), CXCR3⁺CXCR4^- (R3⁺R4^-) PCs measured in each group of samples. Numbers indicate percentages in the drawn gates. (J) Frequency of CXCR3⁺CXCR4^- cells within total circulating PCs. (K) Frequency of CXCR3⁺CXCR4⁺ cells within total circulating PCs. (L) Representative flow cytometry plots (left) and frequency (right) of immature cells (CD38^intCD10⁺IgD⁺CD27^-CD21^low) from total transitional B cells, gated with thick black line. Numbers indicate percentages in the drawn gates. (M) Frequency of naïve B cells (non-PC CD19⁺CD27^-IgD⁺) from CD19⁺CD38^dull B cells in each group of samples. (N) Paired analysis of naïve B cell frequency in COVT1 (n=20) and COVT2 (n=20). Data are presented as individual dots. Bars represent mean ± SEM. Dark-colored dots show ICU patients. Two-tailed Mann-Whitney U test was performed to compare HCs with COVT1 and HCs with COVT2. Wilcoxon matched pairs test was performed to compare COVT1 with COVT2 (*P < 0.05, **P < 0.01, and ***P < 0.001). Unless mentioned otherwise, HCs, n=19; COVT1, n=25; COVT2, n=20.

The tSNE projections also highlighted prominent differences in the transitional and naïve B cell compartments among healthy controls, COVT1 and COVT2 (Fig. 2B-D). Using conventional gating strategies (Fig. S2C) and pairwise longitudinal comparisons, we identified an increased proportion of immature CD21^low transitional B cells (Fig. 2L) as well as a highly significant transient expansion of naïve B cells during the acute phase of infection (Fig. 2M, N). Furthermore, our data revealed a positive correlation between the frequency of naïve B cells and the frequency of CD19⁺ B cells (Fig. S2D), which implies that naïve B cell expansion is responsible for the increased frequency of CD19⁺ B cells. These changes likely reflect an increased inflammation-dependent mobilization of developing B cells from the bone marrow into the periphery, a phenomenon also observed in other infections³¹. Compared to similar cells from healthy controls, naïve B cells from COVT1 exhibited significantly lower expression of CD21, HLA-DR, IgD and the chemokine receptors CCR6, CXCR5 and CXCR3 but not CCR7, presumably due to their activation by signals from the inflammatory environment (Fig. S2E-K). Thus, aside from promoting a transient induction of CXCR3⁺ unswitched and class-switched PBs, SARS-CoV-2 may induce the mobilization of precursors of mature B cells together with the transient expansion and activation of naïve B cells.

SARS-CoV-2 infection drives transient increases of extrafollicular memory switched B cells and long-lasting expansion of IgG1 memory B cells

To gain insight into the dynamics of B cell memory (ME) responses to SARS-CoV-2, we analyzed the frequency of different antigen-experienced B cell populations in healthy controls and COVID-19 patients during acute infection and convalescence. First, we found that the frequency of ME B cells within total CD19⁺ B cells was significantly lower in acute COVT1 compared to COVT2 and healthy controls (Fig. S3A and Fig. 3A), presumably due to the transient expansion of newly formed PBs and naïve B cells. Unlike PCs but similar to naïve B cells, the frequency of CXCR3⁺ ME B cells diminished during the acute phase of infection (Fig. 3B), suggesting a possible impairment in their recruitment into inflamed tissues and lymphoid follicles³². Within the ME B cell compartment, we observed a significant and persistent reduction in the proportion of IgM⁺IgD⁺CD27⁺ ME unswitched B cells in COVID-19 patients (Fig. S3A and Fig. 3C) and a progressive expansion of IgG1-expressing ME B cells (Fig. 3D) within the memory compartment. Interestingly, the frequency of unswitched ME B cells positively correlated with the serum concentration of inflammatory marker ferritin (Fig. 3E) as well as with NP-specific IgM titers (Fig. 3F). Moreover, COVID-19 patients with higher plasma concentration of ferritin and other markers of inflammation (data not shown) also showed a lower frequency of IgG1 ME B cells during acute infection (Fig. 3G). Remarkably, ME B cells expressing other IgG subclasses did not increase after infection (Fig. S3B). These results may reflect the consolidation of efficient anti-viral IgG1-mediated immunity in less severe patients.

• Download figure
• Open in new tab

Fig. 3. COVID-19 is associated with temporary expansion of extrafollicular switched memory B cells, long-lasting contraction of IgM⁺IgD⁺CD27⁺ memory B cells and late increase in IgG1⁺ memory B cells.

(A) Frequency of memory (ME) B cells defined as depicted in Fig. S3 from total CD19⁺ B cells in healthy controls (HCs), COVT1 and COVT2. (B) Representative flow cytometry histogram and frequency of CXCR3⁺ cells from total ME B cells in HCs, COVT1 and COVT2. (C) Frequency of CD38^dullCD10⁻IgM⁺IgD⁺CD27⁺ ME B cells within total ME B cells in HCs, COVT1 and COVT2. (D) Frequency of IgG1⁺ ME B cells within total ME B cells in HCs, COVT1 and COVT2 (left) and paired analysis in COVT1 (n=20) and COVT2 (n=20, right). (E) Spearman correlation analysis of ferritin levels in COVT1 plotted against CD38^dullCD10⁻ IgM⁺IgD⁺CD27⁺ ME B cells within total ME B cells in COVT1. (F) Spearman correlation analysis of NP-specific IgM antibody titers in COVT1 plotted against CD38^dullCD10⁻ IgM⁺IgD⁺CD27⁺ ME B cells within total ME B cells in COVT1. (G) Spearman correlation analysis of IgG1 ME B cells within total ME B cells in COVT1 plotted against ferritin levels in COVT1. (E-G) r stands for Spearman’s rank-order correlation. (H) Representative flow cytometry plots and frequency of DN memory B cells defined as CD38^dullCD10⁻IgD^-CD27^- out of total IgD^- ME B cells. Numbers indicate percentages in the drawn gates. (I) Representative flow cytometry plots and frequency of extrafollicular DN2 ME B cells defined as CD38^dullCD10⁻IgD^-CD27^-CD21^-CD11c⁺ out of total ME B cells. Numbers indicate percentages in the drawn gates. (J) Frequency of DN3 ME B cells defined as CD38^dullCD10⁻IgD^-CD27^-CD21^-CD11c^- out of total ME B cells. (K) Frequency of DN2 ME B cells from total IgG1⁺ ME B cells. Data are presented as individual dots. Dark-colored dots show ICU patients. Error bars represent mean ± SEM. Two-tailed Mann-Whitney U test was performed to compare HCs with COVT1 and HCs with COVT2. Wilcoxon matched pairs test was performed to compare COVT1 with COVT2 (*P < 0.05, **P < 0.01, and ***P < 0.001). Unless mentioned otherwise, HCs, n=19; COVT1, n=25; COVT2, n=20.

In agreement with findings published previously²⁹, we categorized IgD^-CD27^- double negative (DN) ME B cells as early activated CD27^-CD21⁺CD11c^- ME or DN1, extrafollicular CD27^-CD21^- CD11c⁺ PB precursors or DN2, and DN3 CD27^-CD21^-CD11c^- B cells (Fig. S3C). The analysis of these ME B cell subsets revealed a significant and transient increase in the frequency of CD27^-IgD^- DN cells in COVT1 compared to COVT2 and healthy controls (Fig. 3H), as reported by others²⁸. Among DN B cells, DN2 cells are recognized as primed precursors of antibody-secreting cells that differentiate from newly activated naïve B cells through an extra-follicular pathway. In contrast, DN3 has unknown origin and function²⁸. First, we quantified the frequency of DN2 cells within different ME B cells subsets distinguished on the basis of their expressed isotype or subclass. We found a significantly increased proportion of DN2 in ME B cells expressing IgG3, IgG1 and, to lesser extent, IgA1 (Fig. S3D). The frequency of extrafollicular ME DN2 B cells and DN3 cells were higher in COVT1 compared to COVT2 and healthy controls (Fig. 3I, J). Interestingly, the most significant increase in DN2 was observed within the IgG1⁺ ME B cell subset (Fig. 3K). Thus, SARS-CoV-2 infection induces the transient expansion of extrafollicular IgG1⁺ ME B cells, which is followed by a long-lasting IgG1⁺ ME B cell response.

SARS-CoV-2 induces RBD-specific PBs during acute infection followed by sustained expansion of circulating RBD-specific naïve and switched memory B cells in convalescence

Next, a fluorescently labeled recombinant RBD probe was used to identify SARS-CoV-2-specific B cells capable of producing potentially neutralizing antibodies (Fig. S4A). We confirmed the specificity of our assay by comparing our results to those obtained using a double discrimination strategy via the inclusion of two fluorescently labeled RBD probes (Fig. S4B). Antigen-binding CD19⁺ B cells were further characterized according to the expression of CD27, CD21, CD11c, HLADR, IgM, IgD, IgA and the Ig light chain λ (Fig. S4C and Table S3). In agreement with published studies^14,17,33, SARS-CoV-2 infection induced a rapid increase in the frequency of RBD-specific CD19⁺ B cells at T1 that persisted 2 months after viral exposure (Fig. 4A). Interestingly, the phenotypic characterization of these antigen-binding B cells revealed significant changes in virus-specific B cell populations between the acute and convalescent phase of infection (Fig. 4B and Fig. S4D). Consistent with the observed expansion of PCs, COVT1 showed an increased frequency of RBD-specific CD27^highCD21^- PCs compared to healthy controls and COVT2 (Fig. 4C and data not shown). Remarkably, virtually 100% of these RBD-specific cells were newly generated HLADR⁺ PBs (Fig. 4D).

• Download figure
• Open in new tab

Fig. 4. SARS-CoV-2 induces a temporal enrichment of RBD-specific plasmablasts in the early stage of infection and a sustained expansion of RBD-specific naïve and memory B cell subsets.

(A) Representative flow cytometry staining of CD19⁺ RBD-specific B cells (left) and frequency of RBD⁺ cells within total CD19⁺ B cells (right) from each group of samples using a fluorescently labeled RBD probe. Numbers indicate the percentage of RBD-specific cells within total CD19⁺ B cells. Black large dots represent RBD-specific B cells. (B) Relative percentage of B cell subsets within total RBD⁺ CD19⁺ B cells in healthy controls (HCs), COVT1 and COVT2. (C) Representative flow cytometry plots and paired analysis of the frequency of RBD-specific PCs (CD19⁺RBD⁺CD27⁺⁺CD21^-) within total CD19⁺ B cells. In the flow cytometry plot, red large dots represent RBD-specific PCs and numbers indicate percentage of RBD-specific PCs for the drawn sample. (D) Representative flow cytometry plot and relative percentage of total and RBD-specific PCs expressing HLA-DR in COVT1. In the flow cytometry plot, red large dots represent RBD-specific PCs expressing HLA-DR. (E) Frequency and paired analysis of RBD-specific naïve B cells (non-PC CD19⁺RBD⁺CD27^- IgD⁺) within total CD19⁺ in HCs, COVT1 and COVT2. (F) Relative percentage of total or RBD-specific naïve B cells expressing CD11c in COVT1 samples. (G) Representative flow cytometry plots and paired analysis of the frequency of CD11c⁺ cells from RBD-specific naïve B cells. Numbers indicate percentage of RBD-specific CD11c⁺ naïve cells in the drawn gates. (H) Spearman correlation analysis of CD38^dull naïve B cells in COVT1 and RBD⁺ naïve B cells within CD19⁺ cells in COVT2. r stands for Spearman’s rank-order correlation. (I) Relative percentage of total and RBD⁺ naïve B cells expressing lambda or kappa light chains in COVT1. (J) Frequency of RBD-specific ME B cells defined as non-PC CD19⁺RBD⁺IgD^- from total CD19⁺ B cells in HCs, COVT1 and COVT2. (K) Paired analysis of the frequency of RBD-specific IgG⁺ ME B cells defined as IgA^-IgM^-IgD^- within total CD19⁺ in COVT1 and COVT2. (L) Paired analysis of the frequency of RBD-specific IgM⁺IgD⁺ ME B cells defined as non-PC CD19⁺CD27⁺IgM⁺IgD⁺ within total CD19⁺ in COVT1 and COVT2. (M) Representative flow cytometry plots and paired analysis of the frequency of DN2 ME B cells (non-PC CD19⁺RBD⁺CD27^-IgD^-CD11c⁺) from total IgD^- RBD⁺ ME B cells. Numbers indicate percentages of RBD-specific DN2 ME in the drawn gates. Dark-colored dots show ICU patients. Bars represent mean ± SEM. Two-tailed Mann-Whitney U test was performed to compare HCs with COVT1 and HCs with COVT2. Wilcoxon matched pairs test was performed to compare COVT1 with COVT2 (*P < 0.05, **P < 0.01, and ***P < 0.001). HCs, n=11; COVT1, n=19; COVT2, n=19. (B, D, F, I) Data are presented as mean ± SEM.

In addition, RBD-targeting IgM⁺IgD⁺CD27^- naïve B cells increased progressively in COVID-19 patients (Fig. 4E). These cells were enriched in CD11c⁺ early activated B cells in COVT1 and even more in COVT2 (Fig. 4F-G). Of note, the frequency of RBD-specific naïve B cells during convalescence strongly correlated with the frequency of total naïve B cells in the acute phase of infection (Fig. 4H), suggesting that the early mobilization of developing B cells to the periphery and the expansion of naïve B cells could contribute to humoral defense by increasing germline RBD-specific B cells. To elucidate the mechanisms leading to the sustained enrichment of RBD-targeting naïve B cells, we measured the plasma concentration of IL-7, a cytokine involved in B and T lymphopoiesis³⁴. Compared to healthy controls, COVID-19 patients showed more IL-7 during acute infection (Fig. S4E) and its concentration correlated with the frequency of RBD-specific naïve B cells in COVT2 (Fig. S4F), supporting a direct role of this cytokine in enhancing the expansion and recruitment of relatively immature transitional and naïve B cells with useful specificities. Additionally, RBD-specific naïve B cells expressed more Igλ compared to total naïve B cells (Fig. 4I). Given that Igλ light chain has a more varied structural repertoire than Igκ³⁵, Igλ enrichment may provide the host with the advantage of developing a broader humoral response to SARS-CoV-2.

We then analyzed whether RBD-specific ME B cells were induced and maintained throughout the study’s time course. The number of RBD⁺ ME B cells in COVT1 was significantly greater than in healthy controls and increased even further in COVT2 (Fig. 4J). Similarly, there was a sustained increase in the number of RBD-specific IgG (IgM^- IgD^-IgA^-) class-switched ME and, to a lesser extent, in RBD-specific IgM⁺IgD⁺CD27⁺ unswitched ME B cells (Fig. 4K-L). The RBD-specific ME compartment included only a few IgA⁺ and unswitched IgM⁺IgD^- ME cells and these subsets did not increase overtime (Fig. S4G). Consistent with our characterization of B cell subsets, RBD-specific extrafollicular ME B cells, or DN2, expanded during the acute phase of infection but their frequency dramatically decreased in convalescence (Fig. 4M).

To assess the persistency of virus-specific ME B cell responses, RBD-targeting B cell subsets were analyzed by flow cytometry in a small cohort of non-hospitalized convalescent individuals three months (COVT3) and six months (COVT6) PSO (Fig. S5). In agreement with previous reports^14,17,18, RBD-specific B cells persisted up to 6 months after infection (Fig. S5A) and mainly consisted of IgG⁺ ME B cells (Fig. S5B). Accordingly, RBD-specific IgG and, to a lesser extent, IgA were significantly higher in sera from COVT3 and COVT6 compared to healthy controls (Fig. S5C). Thus, SARS-CoV-2 infection promotes a rapid and transient induction of RBD-specific PBs and extrafollicular ME B cells as well as a long-lasting expansion of RBD-targeting activated naïve B cells and IgG class-switched ME B cells.

COVID-19 is associated with a sustained activation of CD4⁺ and CD8⁺ T cells

To further characterize the dynamics of immune responses over time in COVID-19 patients, we simultaneously analyzed by multicolor spectral flow cytometry the phenotypic landscape of other circulating lymphoid and myeloid populations using an in-house developed 23-color antibody panel (Table S4). The projection of the data for CD19^-CD3⁺ cells in tSNE space allowed the definition of the main T cell subsets, the contribution of each group of samples to T cell populations, and the expression patterns of relevant T cell surface proteins (Fig. 5A-C). In parallel, we queried the data by traditional gating (Fig. S6). As reported previously^12,36, the acute phase of infection was characterized by a striking reduction in the frequency of CD3⁺ T cells (Fig. 5D), which was driven by a decrease in the frequency of CD4⁺ and, to an even larger extent, CD8⁺ T cells (Fig. S7A, B). Compared to healthy controls, COVT2 convalescent individuals still had a significantly decreased frequency of CD3⁺ cells (Fig. 5D) and in particular of CD4⁺ T helper cells (Fig. S7A), which suggested their persistent loss from the periphery. We then evaluated how SARS-CoV-2 infection associated with temporary or long-lasting perturbations in both CD4 and CD8 T cell compartments. Consistent with findings published earlier¹², the frequency of naïve T cells within total CD4⁺ T cells in acutely infected patients was significantly lower compared to age-matched healthy controls and further declined during convalescence (Fig. 5E). This reduction was associated with a prominent expansion of activated CD38⁺HLA-DR⁺ CD4⁺ T cells during the acute and, to a lesser extent, convalescent phase of infection (Fig. 5F), as observed in other viral infections^37,38.

• Download figure
• Open in new tab

Fig. 5. COVID-19 is associated with long-lasting contraction of naïve T cell compartment, CD4⁺ and CD8⁺ T cell activation, and expansion of circulating T follicular helper cells.

(A) Merged tSNE projection of CD3⁺ cells for healthy controls (HCs; n=21), COVT1 (n=25) and COVT2 (n=20) donors concatenated and overlaid, with main T cell populations indicated by color. CM, central memory. EM, effector memory. cTfh, circulating T follicular helper T cells. (B) Merged tSNE projection of CD3⁺ T cells with each sample group indicated by color. (C) tSNE projections of expression of the indicated cell surface markers. (D) Frequency of CD3⁺ T cells within live lymphocytes. (E) Frequency of naïve CD4⁺ T cells (CD27⁺CD45RA⁺CCR7⁺) from total CD4⁺ T cells. (F) Representative flow cytometry plots (left) and frequency of CD38⁺HLA-DR⁺ activated CD4⁺ T cells within total CD4⁺ T cells. Numbers indicate percentages in the drawn gates. (G) Frequency of circulating cTfh cells (non-naïve CXCR5⁺PD-1⁺) within total CD4⁺ T cells. (H) Frequency of activated cTfh cells (CD38⁺ICOS⁺) within total cTfh cells. (I) Frequency of naïve CD8⁺ T cells (CD27⁺CD45RA⁺CCR7⁺) within total CD8⁺ T cells. (J) Frequency of CM CD8⁺ (CD27⁺CD45RA^-CCR7⁺) and (K) EM3 (CD27^-CD45RA^-CCR7^-) T cell subsets within total CD8⁺ T cells. (L) Representative flow cytometry plots and frequency of CD38⁺HLA-DR⁺ activated CD8⁺ T cells within total CD8⁺ T cells. Numbers indicate percentages in the drawn gates. (M) Spearman correlation analysis of activated CD4⁺ T and CD8⁺ T cells in all samples analyzed. r stands for Spearman’s rank-order correlation. (A-M) Data are presented as individual dots. Dark-colored dots show ICU patients. Bars represent mean ± SEM. Two-tailed Mann-Whitney U test was performed to compare HCs with COVT1 and HCs with COVT2. Wilcoxon matched pairs test was performed to compare COVT1 with COVT2 (*P < 0.05, **P < 0.01, and ***P < 0.001). HCs, n=21; COVT1, n=25; COVT2, n=20.

We then analyzed the frequency of circulating T follicular helper (cTFH) cells. Tfh cells are a specialized subset of CD4⁺ T cells that provide cognate help to antigen-specific B cells in the germinal center (GC) to initiate and maintain humoral immune responses³⁹. These cells include a circulating counterpart, cTfh cells, which co-express PD-1 and CXCR5 as their GC-based equivalents do (Fig. S6). In general, cTfh cells can be used as a surrogate to evaluate the Tfh cell activity in lymphoid tissues⁴⁰. Of note, activated CD38⁺ICOS⁺ cTfh cells likely reflect a recent exit from the GC immediately after antigen encounter⁴¹. The frequency of cTfh cells was significantly higher in COVT2 patients compared to healthy controls and COVT1 patients, which is consistent with a model of antigen persistency and long-lasting GC reaction previously suggested by others²² (Fig. 5G). Interestingly, recently activated CD38⁺ICOS⁺ cTfh cells were more profoundly increased during the acute phase, which probably reflects the peak of GC responses in the early phase of the infection (Fig. 5H).

We further explored whether SARS-CoV-2 infection promoted perturbations of other CD4⁺ helper T cells subsets, including central memory (CM) T cells, effector memory (EM) 1, 2 and 3 T cells, terminally differentiated effector memory (EMRA) T cells, T regulatory (Treg) cells as well as T helper 1 (Th1), Th1/17, Th2 and Th17 cells (Fig. S7C). We detected an increase in the frequency of Treg cells during the acute phase of infection and a small but significant increase in the frequencies of Th1 and Th1/17 cells in COVT2 compared to COVT1 upon pairwise comparison, whereas Th2 cells showed a decreased frequency over time (Fig. S7C).

Given the role of CD8⁺ T cells in viral infection and consistent with previous reports^12,36, COVID-19 was associated with several changes in the CD8⁺ T cell compartment. We observed a significantly reduced frequency of naïve CD8⁺ T cells during convalescence as well as an expansion of CM and EM3 CD8⁺ T cells (Fig. 5I-K). On the other hand, our analysis showed no differences in the frequency of EM1, EM2 and EMRA CD8⁺ T cells between groups (Fig. S7D). We also found that activated CD38⁺HLA-DR⁺CD8⁺ T cells strongly expanded during the acute phase of infection and still persisted at a higher frequency compared to healthy controls during convalescence (Fig. 5L). As expected, patients with a higher increase in activated CD4⁺ T cells showed a higher frequency of activated CD8⁺ T cells as well (Fig. 5M). Thus, SARS-CoV-2 causes a marked T cell loss during the acute phase of infection, which is associated with long-lasting and coordinated activation of CD4⁺ and CD8⁺ T cells as well as expansion and activation of professional B cell-helping T cells.

SARS-CoV-2 infection drives profound changes in the innate immune compartment

Acute SARS-CoV-2 infection triggers alterations in circulating innate immune cell subsets, some of which have been associated with COVID-19 severity⁴². Nevertheless, how these perturbations persist during convalescence after viral clearance is unclear. We combined global high-dimensional mapping via tSNE of CD19^-CD3^- cells with traditional gating analysis to study significant alterations in the phenotype and frequency of peripheral blood innate immune cell subsets in COVT1 and COVT2 and healthy controls (Fig. 6A-C and Fig. S8A). Although the frequency of total circulating CD11c⁺HLA-DR⁺ myeloid cells among live PBMCs in COVT1 and COVT2 was comparable to that of healthy controls (Fig. S8B), we observed significant differences in the phenotype as well as subset composition of these myeloid cells. We initially explored the expression of the activation-induced molecule CD38 and HLA-DR on circulating myeloid cells. We found a striking induction of CD38 and a reduction in the expression of HLA-DR during acute infection (Fig. 6D-E). Of note, the decrease in HLA-DR may interfere with proper antigen presentation and has been directly linked to an immunosuppressive phenotype of monocytes during sepsis⁵. Indeed, HLA-DR levels strongly correlated with disease severity defined according to the patient’s oxygen requirement (Fig. S8C). However, these changes were only transient, as CD38 and HLA-DR returned to their homeostatic expression levels in the convalescent phase (Fig. 6D-E).

• Download figure
• Open in new tab

Fig. 6. SARS-CoV-2 infection associates with transitory and long-lasting changes in the innate immune compartment.

(A) Merged tSNE projection of CD19^-CD3^- cells for healthy controls (HCs; n=21), COVT1 (n=25) and COVT2 (n=20) samples concatenated and overlaid, with main innate immune cell populations indicated by color. pDCs, plasmacytoid dendritic cells. mDCs, myeloid dendritic cells. Mo, monocytes. NK, natural killer cells. (B) Merged tSNE projection of CD19^-CD3^- innate immune cells with each sample group indicated by color. (C) tSNE projections of indicated cell surface markers. (D) Geometric mean fluorescence intensity (GeoMFI) of CD38 and (E) HLA-DR expression in CD19^-CD3^-CD11c⁺ myeloid cells from HCs, COVT1 and COVT2. (F) Representative flow cytometry plots and frequency of non-classical CD14^dimCD16⁺ Mo as a proportion of live cells. Numbers indicate percentages in the drawn gates. (G) GeoMFI of CD38 in CD14^dimCD16⁺ Mo. (H) Frequency of mDCs, (I) basophils, (J) pDCs, and (K) NK cells within total live cells. (L) GeoMFI of CD38 in NK cells. (M) Representative flow cytometry plots and frequency of CD56^brightCD16^- NK cells. Numbers indicate percentages in the drawn gates. Data are presented as individual dots. Dark-colored plots show ICU patients. Bars represent mean ± SEM. Two-tailed Mann-Whitney U test was performed to compare HCs with COVT1 and HCs with COVT2. Wilcoxon matched pairs test was performed to compare COVT1 with COVT2 (*P < 0.05, **P < 0.01, and ***P < 0.001). HCs, n=21; COVT1, n=25; COVT2, n=20. (D,E,G and L) Data are presented in boxplots. Box boundaries represent the 1^st and 3^rd quartile of the distribution, while the center line represents the 2nd quartile (median). Whiskers go down to the smallest value and up to the largest.

We then analysed the frequency of CD14⁺CD16^- conventional monocytes, CD14⁺CD16⁺ intermediate monocytes and CD14^dimCD16⁺ non-conventional monocytes (Fig. 6A-C and Fig. S8A). Classical CD14⁺CD16^- monocytes remained unaltered, whereas intermediate CD14⁺CD16⁺ monocytes showed a higher frequency in COVT2 compared to COVT1 and healthy controls (Fig. S8D-E). Furthermore, we observed a robust but transient decrease in the frequency of non-conventional CD14^dimCD16⁺ monocytes (Fig. 6B, C and F), which appeared to be strongly activated during acute infection (Fig. 6G). Non-conventional monocytes have a patrolling function and contribute to the antiviral response⁴³, which may point to their selective recruitment to the inflamed lung upon FcγRIIIA-mediated activation by virus-IgG immune complexes. Accordingly, non-conventional monocytes have been shown to be enriched in the lungs of critical COVID-19 patients⁴⁴. The frequency of circulating myeloid dendritic cells (mDCs) and basophils were also strongly reduced in COVT1 compared to healthy controls, but both returned to homeostatic levels two months PSO (Fig. 6H-I). Interestingly, we detected increased CD25 expression on basophils from COVT1 compared to COVT2, which suggests basophil activation during the acute phase of infection (Fig. S8F).

Finally, we analyzed the frequency of circulating plasmacytoid DCs (pDCs) and natural killer (NK) cells, two subsets of the innate immune system that play crucial protective roles in viral infections. As reported previously⁴⁵, the frequency of pDCs significantly decreased during acute infection (Fig. 6J). Remarkably, convalescent individuals maintained significantly lower frequency of pDCs compared to healthy controls, even after viral clearance (Fig. 6J), suggesting persistent impaired type I IFN responses in infected individuals.

Although our data did not reveal any significant changes in the frequencies of natural killer (NK) cells within live cells, the analysis of CD38 expression indicated an activated phenotype during the acute phase of the infection (Fig. 6K, L). Moreover, the characterization of NK cell subsets based on the relative expression of CD56 and CD16 showed a persistent decrease in the frequency of circulating immature CD56^brightCD16^- cells and an increase in the proportion of naturally cytotoxic CD56^dimCD16^high cells within total NK cells (Fig. 6M and Fig. S8G). Thus, SARS-CoV-2 infection drives profound changes in the circulating innate immune compartment, some of which persist after viral clearance during convalescence.

High dimensional analysis reveals coordinated immune responses and their association with disease severity

The immune system is composed of layered defense mechanisms of increasing specificity that protect the host from infection. To investigate the putative associations between innate and adaptive immune trajectories developed upon SARS-CoV-2 infection, we performed pairwise correlations across 41 variables that were significantly different among the groups analyzed, including virus-specific antibody titers measured by ELISA and selected immune parameters identified by high-dimensional flow cytometry (Data File S1). We found that the global antibody response to SARS-CoV-2 positively correlated with the frequency of circulating PCs, extrafollicular DN2 ME B cells, activated CD4⁺ T cells, cTfh cells and CD8⁺ T cells. It also positively correlated with CD38 expression on circulating myeloid cells, which suggests a close coordination between innate and adaptive immune responses. On the other hand, the level of virus-specific antibody titers, the percentage of circulating PCs and the frequency of CD19⁺ cells inversely correlated with the frequency of CD3⁺ T cells, non-conventional monocytes and pDCs, which suggests virus-induced interference with different components of both adaptive and innate immune responses in a fraction of COVID-19 patients (Fig. 7A). Next, principal component analysis (PCA) was performed to examine the general distribution of healthy controls and infected individuals during acute and convalescent stage according to these immune parameters. The analysis revealed that patients during acute infection were notably more distant from healthy controls than same patients during convalescence in PCA space (Fig. 7B). Virus-specific IgA titers, as well as the frequency of IgG and IgM PCs, activated CD4 T cells and pDCs, were among the variables that mostly contributed to the observed clustering in PC1 space. The frequency of naïve B cells and, to a lesser extent, IgM⁺IgD⁺ ME B cells, IgA1 PCs and NP-specific IgM titers were the immune features that mostly contributed to the clustering in PC2 space (Fig. 7C, D). Remarkably, when we computed the Euclidean distance from each COVID-19 patient at T1 to the centroid of healthy controls, we observed a positive correlation between Euclidean distance and disease severity defined according to the patient’s oxygen requirement and serum concentration of ferritin and lactate dehydrogenase (Fig. 7E). Of note, the same analysis performed on samples from COVT2 did not show equivalent correlation patterns (data not shown), which suggests that immunological recovery trajectories occur independently of disease severity status. Thus, cross-dataset correlation and principal component analysis reveals coordinated immune responses in COVID-19 patients as well as immune parameters associated with disease severity and inflammation.

• Download figure
• Open in new tab

Fig. 7. High-dimensional analysis of all variables studied reveals coordinated immune responses and their association with disease severity

(A) Spearman correlation mapping of indicated parameters for healthy controls (HCs), COVT1 and COVT2. Spearman’s rank correlation coefficient (ρ) was indicated by heat scale; Spearman p-value significance levels were corrected using Benjamini-Hochberg method significance (*P < 0.05, **P < 0.01, and ***P < 0.001). PCs, plasma cells. R3, CXCR3. R4, CXCR4. ME, memory. cTfh, circulating T follicular helper cells. CM, central memory. EM3, effector memory type 3. pDCs, plasmacytoid dendritic cells. mDCs, myeloid dendritic cells. Baso, basophils. NK, natural killer cells. trans, transitional B cells. (B) Principal component analysis (PCA) based on antibody titers, frequencies and GeoMFI of all measured markers (centered and scaled). RBD-specific parameters were excluded from the analysis. The first component explains 25% of the variation, whereas the second component explains 12.3% of the variation. Each dot represents a donor and each color represents a donor group. Dark-colored dots show ICU patients. Confidence ellipses for each group are plotted (confidence level set to 95%). (C) Variable correlation plot showing eigenvector-based coordinates of the top 20 variables in the two-dimensional space defined by the first two principal components. The relative position of the clinical variables reflects their relationship (positive correlated variables point to the same side of the plot; negative correlated variables point to opposite sides of the plot), while the length of the arrow is proportional to their contribution to the principal components. (D) Relative contribution of the top 10 variables to PC1 (top panel) and PC2 (bottom panel). (A-D) HCs, n=16; COVT1, n=25; COVT2, n=20. (E) Spearman correlation analysis of oxygen flow, ferritin and lactate dehydrogenase (LDH) levels in COVID-19 patients at COVT1 (n=25) plotted against the Euclidean distance from each COVT1 patient to the centroid of the HC group. The centroids by group were computed by averaging PC1 and PC2 for each member of group. The oxygen flow ranges from 0, meaning no need for supplemental oxygen, to 8, meaning the need for orotracheal intubation and mechanical ventilation. Data are presented as individual dots. Dark-colored dots show ICU patients. r stands for Spearman’s rank-order correlation.

Experimental Model and Subject Details Study Cohort

Blood samples were collected from COVID-19 patients (n = 25) in acute phase of infection (COVT1) hospitalized at the Hospital del Mar (Barcelona, Spain), with patient informed consent. Only patients with confirmed SARS-CoV-2 infection by reverse transcription quantitative polymerase chain reaction (RT-qPCR) of nasopharyngeal swab were included. All COVT1 samples were collected within the first 14 days following symptom onset, being the median time from symptom onset to blood testing 11 days [8-14]. A follow-up blood sample was collected from 20 out of the 25 initial participants two months after first sample collection (COVT2), being the median time from symptom onset to COVT2 blood testing 70 days [59-84]. The median age of patients was 51 years old. 52 % were males. The median length of stay from hospitalization to discharge was 16 days [5-90]. Clinical laboratory data were collected from the date closest to the first day of hospitalization. Blood samples from age-matched healthy donors (n = 21) were collected. All research subjects were pseudonymized by replacing directly identifiable variables with an ID immediately when collecting data. Patients ID were not known to anyone outside the research group. Demographic and clinical data of patients included in the study are summarized in Table S1. All procedures followed were approved by the Ethical Committee for Clinical Investigation of the Institut Hospital del Mar d’Investigacions Mèdiques (Number 2020/9189/I).

Sample collection and processing

For all COVID-19 patients and healthy volunteers, sera were collected from whole blood in silica-treated tubes (BD Biosciences) where the blood was incubated for 30 min without movement to trigger coagulation. Next, samples were centrifuged for 10 min at 1300 g at room temperature (RT), heat-inactivated at 56°C for 1 hour and stored at −80°C prior to use. PBMCs were isolated from whole blood collected with EDTA anticoagulant via Ficoll-Paque Premium (Cytiva) following manufacturer’s instructions. Briefly, each sample was diluted with an equal volume of phosphate buffered saline (PBS). The diluted blood was carefully layered over an equal volume of Ficoll solution in a Falcon tube, followed by centrifugation at 400 g for 30 min at 20°C without brake. The upper layer containing plasma was collected and stored at −80°C. The layer of mononuclear cells was gently removed and washed twice with PBS. Pelleted cells were counted using Turk solution, resuspended in fetal bovine serum (FBS, Gibco) with 10% Dimethyl sulfoxide (DMSO, Sigma) and stored in liquid nitrogen prior to use.

Production of recombinant SARS-CoV-2 proteins

The pCAGGS RBD construct, encoding for the receptor-binding domain of the SARS-CoV-2 Spike protein (amino acids 319-541 of the Spike protein) along with the signal peptide plus a hexahistidine tag was provided by Dr Krammer (Mount Sinai School of Medicine, NY USA). The pLVX-EF1alpha-nCoV2019-N-2xStrep-IRES-Puro construct, encoding for the full-length SARS-CoV-2 nucleocapsid protein (NP) fused to a double Strep-tag at the C-terminus was a gift from Dr Krogan (University of California, San Francisco USA). Recombinant proteins were expressed in-house in Expi293F human cells (Thermo Fisher Scientific) by transfection of the cells with purified DNA and polyethylenimine (PEI). For secreted RBD proteins, cells were harvested 3 days post transfection and RBD-containing supernatants were collected by centrifugation at 13000rpm for 15min. RBD proteins were purified in Hitrap-ni Columns in an automated Fast Protein Liquid Chromatography (FPLC; Äkta avant), concentrated through 10 kDa Amicon centrifugal filter units (EMD Millipore) and resuspended in phosphate buffered saline (PBS). For NP, cell lysates from transfected cells were centrifuged at 13000 rpm for 20 min and the supernatant was loaded into a 5 ml Strep-Tactin XT column. The eluted nucleocapsid was concentrated through a 10 kDa Amicon centrifugal filter unit and purified with SEC (Sephadex 10/300) in PBS.

Enzyme-linked immunosorbent assay (ELISA)

ELISAs performed in this study were adapted from previously established protocols⁴⁶. Serum samples were heat-inactivated at 56°C for 1 hour and stored at −80°C prior to use. 96-well half-area flat bottom high-bind microplates (Corning) were coated overnight at 4°C with each respective recombinant viral protein at 2µg/ml in PBS (30 µl per well) or with PBS alone. Plates were washed with PBS 0.05% Tween 20 (PBS-T) and blocked with blocking buffer (PBS containing 1.5% Bovine serum albumin, BSA) for 2 hours at RT. Serum samples were serially diluted in PBS supplemented with 0.05% Tween 20 and 1% BSA and added to the viral protein- or PBS-coated plates for 2 hours at RT. After washing, plates were incubated with horseradish peroxidase (HRP)-conjugated anti-human Ig secondary antibodies diluted in PBS containing 0.05% Tween 20 1% BSA for 45 minutes at RT. Plates were washed 5 times with PBS-T and developed with TMB substrate reagent set (BD bioscience) with development reaction stopped with 1M H2SO4. Absorbance was measured at 450nm on a microplate reader (Infinite 200 PRO, Tecan). To detect RBD-specific and NP-specific IgM, IgA and IgG, HRP-conjugated anti-human IgA, IgG and IgM (Southern Biotech) were used at a 1:4000 dilution. To assess the distribution of the different IgG antibody subclasses, HRP-conjugated anti-human IgG1, IgG2, IgG3 and IgG4 (Southern Biotech) were used at a 1:3000 dilution. To detect SARS-CoV-2-specific IgA1 and IgA2, HRP-conjugated anti-human IgA1 or IgA2 (Southern Biotech) were used at a dilution of 1:2000 and 1:4000, respectively. To quantitate the level of each viral antigen-specific antibody class or subclasses optical density (OD) values were calculated after subtraction of background (OD450 of serum dilutions on PBS-coated plates) and the area under the curve (AUC) derived from optical density measurements of six serial dilutions was determined using Prism 8 (GraphPad). AUC values below an established cutoff were replaced by the AUC value of the cutoff for plotting and calculation purposes. Negative threshold values were set using healthy control AUC levels plus 2 times the standard deviations of the mean.

The concentration of interleukin-7 (IL-7) was measured by a sandwich ELISA kit (Peprotech) following the manufacturer’s instructions on 50 µl plasma samples. The levels of IL-7 were expressed as pg/mL of plasma.

Spectral flow cytometry

For spectral flow cytometry, frozen Peripheral Blood Mononuclear Cells (PBMCs) were thawed in warm complete RPMI 1640 medium (Biowest) and centrifuged at 1450 rpm for 5 min at RT. Cells were then resuspended in 2 mL of LIVE/DEAD Fixable Yellow Dead Cell Stain Kit (Thermo Fisher Scientific) 1:20000 in PBS, incubated for 30 min at RT and stained with two different fluorophore-conjugated antibody cocktails (Table S2 and S4). Only for MIX 3, cells were washed and resuspended in 50 µL cool staining buffer with Fc-receptor blocking reagent for 10 min at 4 °C. For staining, a sequential approach was followed. Anti-CCR7 antibody was added first and incubated for 10 min. Then, all other anti-chemokine receptor antibodies (CXCR3, CXCR4, CXCR5, CCR4 and CCR6 for MIX3, and CXCR3 and CXCR4 for MIX 1) were added to the corresponding tubes and incubated for 10 min. After that, the remaining antibodies were added and incubated for 20 min (MIX 1) or 30 min (MIX 3). Stained cells were washed and resuspended in 4% paraformaldehyde in PBS while vortexing, and incubated for 12 min in ice. Fixed cells were then washed, resuspended in PBS and acquired using the Aurora spectral analyzer (Cytek). Data was analyzed using FlowJo V10.6.2 software (TreeStar).

Flow cytometric detection of SARS-CoV-2-RBD-specific B cells

1.5 million frozen PBMCs were thawed, centrifuged and resuspended in a 1:20000 PBS dilution of LIVE/DEAD Fixable Yellow Dead Cell Stain Kit (Thermo Fisher Scientific) for 30 min at RT to exclude dead cells. Cells were then washed with PBS and incubated with 13.2 pmol SARS-CoV-2 RBD-HIS Biotinylated Recombinant Protein (Sino Biological) diluted in 50 µl 1X PBS supplemented with 0.2% BSA and 2 mM EDTA for 30 min in ice. Cells were washed again and stained with the MIX 2 antibody cocktail (Table S3) plus 1.32 pmol Streptavidin Alexa Fluor 647 (Thermo Fisher Scientific) for 30 min. Stained cells were washed, fixed with 4% paraformaldehyde in PBS for 12 min and acquired with LSR Fortessa (BD Biosciences). Data were further analyzed with FlowJo V10.6.2 software.

Alternatively, 6.6 pmol SARS-CoV-2 RBD-HIS Biotinylated Recombinant Protein was incubated for 1 h in PBS with 0.94 pmol Streptavidin Alexa Fluor 647 and 0.94 pmol Streptavidin Alexa Fluor 488 (Thermo Fisher Scientific), separately. Meanwhile, 1.5 million frozen PBMCs were thawed, centrifuged and resuspended in 50 µl PBS. Next, cells were incubated in the same staining tube with the labeled RBD probes for 20 min in ice. Cells were then washed and stained with the MIX 2 antibody cocktail for 10 min using anti-human IgA AmCyan instead of anti-human IgA FITC (Table S3) and DAPI fluorescent dye (Sigma Aldrich). Stained cells were washed and resuspended in PBS 2%FBS and acquired on LSR Fortessa (BD Bioscience).

High-dimensional data analysis of flow cytometry data

t-Distributed Stochastic Neighbor Embedding (tSNE) analyses were performed with Flowjo V10.6.2 software. tSNE analysis was performed using equal sampling from each FCS file, with 1000 iterations, a perplexity of 30, a Barnes-Hut gradient algorithm and exact KNN algorithm. For CD19⁺ cells, the following samples were concatenated: COVT1/T2 (n = 16) and healthy controls (n = 11); and the following markers from MIX 1 were used to generate the tSNE maps: CXCR4, CXCR3, CD10, CD19, IgM, CD27, CD21, IgA1+IgA2, IgD, HLA-DR, CD38, CD24, CD43, CD11c, CD45RB and LIVE/DEAD YELLOW. For both CD19^-CD3⁺ and CD19^-CD3^- cells, all samples were included in the analysis and the following markers from MIX 2 were used: CCR7, CXCR5, CCR6, CXCR3, CCR4, CD16, CD3, CD8a, CD11c, CD123, CD56, PD-1, CD14, CD45RA, HLA-DR, ICOS, CD4, CD25, CD19, CD27, CCR4, CD127, CD38 and LIVE/DEAD YELLOW.

Data analysis and visualization

GraphPad Prism (version 8.0) and R (version 3.6.3, R Core Team (2019). R: A Language and Environment for Statistical Computing) were used to conduct statistical analyses. For each experiment, the type of statistical test, summary statistics and levels of significance were specified in the figures and corresponding legends. The correlation analysis was performed in R. We performed data imputation of the missing variables using predictive mean matching implemented in the mice package⁵⁹. Pairwise Spearman’s rank correlations between 41 variables (including antibody titers measured by ELISA and immune parameters defined by high-dimensional flow cytometry; Data File S1) for COVT1 (n = 25), COVT2 (n = 20) and healthy controls (n = 16) were calculated and visualized as a correlogram using R package corrplot (Taiyun Wei and Viliam Simko (2017). R package “corrplot": Visualization of a Correlation Matrix). Spearman’s rank correlation coefficient (ρ) was indicated by heat scale; significance was indicated by *P < 0.05, **P < 0.01, and ***P < 0.001; FDR correction was performed using the Benjamini-Hochberg procedure at the FDR < 0.05 significance threshold.

Principal component analysis (PCA) was used to identify the most important features from 41 variables (including antibody titers and immune parameters; Data File S1) using COVT1 (n = 25), COVT2 (n = 20) and healthy controls (n = 16). The PCA was conducted using the “prcomp” function in R and visualized using the “factoextra” package (Kassambara A, and Mundt F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. 2020).