Natural killer cell immunotypes related to COVID-19 disease severity

Study design, clinical cohort, and general NK cell activation in COVID-19

Twenty-seven in-hospital patients with COVID-19 disease (10 moderate and 17 severe) were prospectively recruited after admission to the Karolinska University Hospital as part of a larger immune cell atlas effort (Karolinska COVID-19 Immune Atlas). Strict inclusion and exclusion criteria were used to ensure the setup of homogeneous patient cohorts, including time since symptom debut in relation to hospital admission, and disease severity staging (Fig. 1, A and B). Seventeen healthy controls that were SARS-CoV-2 IgG seronegative and symptom-free at time of sampling were included as controls. To minimize inter-experimental variability and batch effects between patients and controls, all samples (n = 44) were acquired, processed, and analyzed fresh on three consecutive weeks in April and May of 2020 at the peak of the COVID-19 pandemic in Stockholm, Sweden (Fig. 1B). More detailed patient characteristics are provided in the Materials and Methods and tables S1 and S2.

To study the NK cell response, peripheral blood mononuclear cells were analyzed with a 28-color NK cell-focused flow cytometry panel (table S3, gating strategy in fig S1A). No change in total NK cell percentages, or the percentage of CD56^dim and CD56^bright NK cells was noted in patients compared to controls. However, the absolute counts of total NK cells, CD56^dim, and CD56^bright NK cells were severely reduced in patients compared to controls (Fig. 1C, fig. S1A). NK cell activation was assessed by analyzing expression of K_i-67, HLA-DR, and CD69 (Fig. 1D). A robust NK cell activation was noted in COVID-19 patients, where the degree of activation was similar in moderate and severe patients (Fig. 1, E and F). CD56^dim NK cells undergo continuous differentiation starting from less differentiated NKG2A⁺ and CD62L⁺ cells transitioning toward more terminally differentiated CD57⁺ NK cells (26, 27). A higher percentage of both NKG2A⁺ and CD62L⁺ CD56^dim NK cells expressed Ki67 compared to NKG2A^– and CD62L^– cells in patients with COVID-19 (Fig. 1, G and H). A subsequent Boolean analysis, simultaneously taking co-expression patterns of NKG2A, CD62L, and CD57 into account, revealed that the NK cell response in patients with COVID-19 occurred primarily among less differentiated CD62L⁺ NK cells (fig. S1B).

Taken together, using strict inclusion and exclusion criteria and patient recruitment during a defined period of time, we obtained a well-controlled cohort of patients with moderate and severe COVID-19. Peripheral blood NK cells in these patients were robustly activated compared to healthy controls. The general degree of activation was however not directly associated with disease severity.

Phenotypic assessment of NK cells in COVID-19

We next performed a detailed phenotypic assessment of CD56^bright and CD56^dim NK cells in the COVID-19 patients. By principal component analysis (PCA) of CD56^bright (including 26 phenotypic parameters) and CD56^dim NK cells (27 phenotypic parameters), patients clustered separately from healthy controls (Fig. 2, A and B). This was, among other phenotypic markers, driven by changes in expression of CD98, K_i-67, and Ksp37 in CD56^bright (fig. S2A), and Tim-3, CD98, CD38, CD69, and Ksp37 in CD56^dim NK cells (fig. S2B). Flow cytometry analysis by manual gating further revealed increased expression of perforin, Ksp37, MIP-1β, CD98, Tim-3, and granzyme B in CD56^bright NK cells from COVID-19 patients, as well as CD98, Tim-3, NKG2C, MIP-1β, and CD62L, among other markers, in CD56^dim NK cells (Fig. 2, C and D, fig. S2, C and D). In line with the activation profile (Fig. 1, E and F), PCA did not highlight an obvious separation between moderate and severe COVID-19 patients (Fig. 2B), and similar results were obtained when unsupervised hierarchical clustering was performed using expression of all studied phenotypic markers within CD56^bright and CD56^dim NK cells (Fig. 2, E and F). Of note, the clustering clearly distinguished healthy controls from COVID-19 patients. Finally, the phenotypes of responding (K_i-67 positive) compared to non-responding (K_i-67 negative) NK cells were compared. This analysis revealed that responding CD56^bright NK cells co-expressed higher levels of granzyme B, CD25, HLA-DR, and Ksp37 when compared to non-responding cells and that responding CD56^dim NK cells expressed higher levels of Tim-3 and HLA-DR (fig. S2E).

The activated status of NK cells was also confirmed in bronchoalveolar lavage (BAL) from COVID-19 patients by analysis of publicly available single cell RNA sequencing data (scRNA-seq) (Fig. 2, G and H, fig. S2F) (24). Clustering of differentially expressed genes (DEGs) in moderate and severe patients followed by gene ontology enrichment analysis of DEGs revealed six distinct gene modules where effector functions/chemotaxis, activation/proliferation, and interferon response were enriched in patients compared to controls (Fig. 2I, fig. S2, G and H, table S4). Intriguingly, the interferon response and activation/proliferation signatures were most highly enriched in BAL NK cells from patients with moderate disease whereas severe patients had a higher effector function gene signature (Fig. 2I, fig. S2, G and H). Gene set enrichment analysis further highlighted NK cells from BAL fluid of severe patients to display an activated and inflamed profile (Fig. 2J, fig. S2I). Finally, several of the proteins that were increased in peripheral blood NK cells from COVID-19 patients (Fig. 2, C-F, fig. S2, C and D) were also identified as DEGs in BAL NK cells of COVID-19 patients as compared to controls, including GZMB (granzyme B), PRF1 (perforin), HAVCR2 (Tim-3), and CCL4 (MIP-1β) (Fig. 2K).

In summary, peripheral blood CD56^bright and CD56^dim NK cells displayed an activated effector phenotype that is similar in nature in moderate and severe COVID-19 disease, and that could be corroborated by analysis of NK cells from BAL fluid of COVID-19 patients.

Impact of KIR expression and NK cell education on the NK cell COVID-19 response

CD56^dim NK cell function is regulated by inhibitory receptors. Subsets of these cells expressing inhibitory receptors in the presence of self-HLA class I ligands are educated and thus more functional (28, 29). We next evaluated the impact of expression of inhibitory KIRs and NKG2A as well as of NK cell education on the acute COVID-19 response (fig. S3A, table S5). No significant differences were found in the fractions of NK cells expressing different combinations of inhibitory KIRs and NKG2A in COVID-19 patients as compared to controls (fig. S3B). Furthermore, the size of CD56^dim NK cell subsets educated via KIRs or NKG2A was also similar in patients and controls (fig. S3C). Finally, we assessed the NK cell response by measuring K_i-67 up-regulation in CD56^dim NK cell subsets, either solely based on KIR and NKG2A expression or also by integrating information on MHC class I ligand haplotypes (fig. S3, D-F). This comparison indicated that the NK cell response in acute COVID-19 occurred independently of inhibitory KIR expression and NK cell education status.

Increased frequency of adaptive NK cells in severe COVID-19

Our analyses so far revealed robust activation of NK cells with an effector phenotype in COVID-19 patients. The response appeared largely independent of disease severity as few differences were observed when comparing moderate and severe patients. However, we noticed that severe patients had increased frequencies of NKG2C⁺ NK cells (fig. S2D). This indicated an increased frequency of adaptive NK cells, which are often associated with CMV infection and can expand when CMV-seropositive patients undergo other severe acute viral infections (13, 30, 31). Adaptive NK cells are characterized by a terminally differentiated phenotypic signature, and the vast majority of them co-expresses NKG2C and CD57 (30). In severe, but not in moderate COVID-19 patients, we detected higher frequencies of NKG2C⁺CD57⁺ CD56^dim NK cells, as compared to healthy controls (Fig. 3, A and B). Moreover, while the absolute number of NKG2C^–CD57^– cells decreased in COVID-19 patients, NKG2C⁺CD57⁺ NK cells remained stable even in patients with severe disease (Fig. 3, C and D). A deeper phenotypic analysis of NKG2C⁺CD57⁺ NK cells, as compared to their NKG2C^–CD57^– counterparts, highlighted the modulation of several maturation and activation-related markers, for instance NKG2A and CD38 (Fig. 3E, fig S4A). By applying an arbitrary threshold based on the percentage of NKG2C⁺CD57⁺ NK cells over the CD56^dim population (here set to >5%, fig. S4B), we were able to identify patients where NKG2C⁺CD57⁺ NK cells displayed a phenotypic pattern largely overlapping with that observed in adaptive NK cells (Fig. 3E) (30, 32). The expansion of adaptive NK cells was confined to CMV seropositive individuals (Fig. 3F, table S5). Our strategy to identify expansions in the smaller healthy control group showed good correspondence with a previous study assessing similar expansions in a large cohort of more than 150 healthy individuals (Fig. 3G) (30). Strikingly, approximately two-thirds of CMV seropositive severe COVID-19 patients had such an expanded adaptive NK cell population compared to one third in controls and even fewer in the moderate COVID-19 patients (Fig. 3, F and G, fig S4C). Thirteen out of 16 severe COVID-19 patients were negative for CMV DNA in serum (table S5), indicating that the increased proportion of adaptive NK cells in severe patients did not depend on CMV reactivation secondary to COVID-19. In line with this hypothesis, neither the percentage nor absolute counts, or the fraction of proliferating NKG2C⁺CD57⁺ NK cells showed any correlation with CMV IgG levels detected in patient sera (fig. S4D and table S5). Notably, the clinical picture was similar in severe patients with and without adaptive NK cell expansions (fig. S4, E-F, table S6, S7).

NKG2C⁺CD57⁺ NK cells found in COVID-19 patients displayed signs of proliferation, although at lower levels compared to NKG2C^–CD57^– NK cells (Fig. 3, H and I, fig. S4, H-J). Responding NKG2C⁺CD57⁺ NK cells expressed higher levels of HLA-DR, CD38, CD62L, and MIP-1β, but showed lower expression of NKG2A, NKG2D, TIGIT, and CD25 as compared to non-responding NKG2C⁺CD57⁺ NK cells (Fig. 3J, fig. S4K). Next, we set out to more specifically address what was driving the emergence of adaptive NK cell expansions. Reanalysis of published scRNA-seq data revealed that both immune and non-immune cell types displayed up-regulated expression of genes encoding HLA-E and HLA-A, B, and C in BAL fluid of moderate and severe COVID-19 patients as compared to controls (Fig. 3K and fig. S4L). The number of NKG2C⁺CD57⁺ NK cells was inversely correlated with several soluble factors, including LT-α, IL-12β, IFN-γ, and amphiregulin (AREG) (Fig. 3L). However, when specifically comparing soluble factors in severe patients with and without expansions, only a limited number of factors were differentially expressed (fig. S4M). In line with this, associations between phenotypic features and soluble factors were in most cases similar for adaptive and non-adaptive NK cells (fig. S4N).

Taken together, severe, but not moderate, COVID-19 disease is associated with higher frequencies of adaptive NK cells that display signs of proliferation and activation without detectable concurrent CMV reactivation.

Dimensionality reduction analysis of NK cells separates severe from moderate patients

Except for adaptive NK cell expansions being more prevalent in severe COVID-19 patients, manual gating-based flow cytometry analysis did not identify features specific to either moderate or severe COVID-19 patients. In an additional attempt to identify such features through an unsupervised approach, we performed Uniform Manifold Approximation and Projection (UMAP) analysis on all patients and controls (Fig. 4, A and B, fig. S5). This revealed distinct topological regions between patients and controls (Fig. 4A). We applied PhenoGraph to our samples to further describe NK cell subpopulations and to quantify differences in their relative abundance between COVID-19 patient groups. PhenoGraph clustering identified 36 distinct clusters (Fig. 4C). Approximately half of these clusters were equally shared between controls, and moderate and severe patients (Fig. 4D), and some clusters were patient specific. Furthermore, this analysis also revealed clusters that contained NK cells predominantly from controls, moderate, or severe COVID-19 patients (Fig. 4, C and D). The phenotype of NK cells in clusters with similar relative abundance across the three groups was homogenous, but greater differences were observed when assessing clusters uniquely and significantly enriched in moderate or severe patients (Fig. 4, E and F, fig. S5D).

Next, each PhenoGraph cluster was stratified within patients based on six clinical categorical parameters (viremia, seroconversion, oxygen need, steroid treatment, underlying comorbidities, and outcome) (Fig. 4G) in the same manner as moderate and severe patient categories were defined (Materials and Methods). Here, most clusters showed a significantly different relative abundance between patient categories and healthy controls in at least one of the clinical parameter-defined patient groups (Fig. 4G, fig. S6A, table S8). Importantly, for some clusters, the relative abundances were different depending on the clinical parameter. This approach was critical for our understanding, as it enabled us to step away from looking at disease severity on the whole and to instead look into the defining clinical parameters separately. For example, cluster 21 with high CD25 and CD98 expression (Fig. 4H) was almost completely absent from patients who were on mechanical ventilation and the ones who died from COVID-19. Conversely, CD69-expressing cluster 22, and the four most highly K_i-67-expressing clusters accounted for a similar percentage across all patient categories in our dataset (Fig. 4H, fig. S6, B-F). Hierarchical clustering of PhenoGraph clusters’ relative abundance across clinical categorical parameters revealed the presence of two “immunotypes” in patients. One of these contained categorical parameters attributed to milder disease course, and the other one was enriched for viremia, chronic underlying diseases, steroid treatment, mechanical ventilation, and fatal outcome (Fig. 4I). Indeed, a similar clustering was obtained when moderate and severe disease status was included (fig. S6G). Finally, the NK cell phenotypes associated with moderate (immunotype 1) versus severe (immunotype 2) disease were assessed. This revealed up-regulated MIP-1β, CD98, and TIGIT in the moderate immunotype and conversely high expression of perforin, NKG2C, and Ksp37 in the severe immunotype, while the activation markers Ki67 and CD69 remained unchanged (Fig. 4J). Interestingly, high perforin, NKG2C, and Ksp37 are features of adaptive NK cell expansions which were also associated with severe disease (Fig. 3).

In summary, unsupervised high-dimensional analysis of the NK cell response in COVID-19 disentangles NK cell phenotypes associated with disease severity.

CD56^bright NK cell arming is associated with disease severity

Finally, we performed a more detailed analysis of the NK cell response in relation to clinical laboratory parameters and clinical scoring systems. PCA of moderate and severe patients based on clinical laboratory parameters revealed separation between the two groups (Fig. 5A). This was primarily driven by high systemic levels of neutrophils, IL-6, D-dimer, and C-reactive protein, clinical laboratory data typically associated with severe disease (Fig. 5B). Out of the altered NK cell phenotypic parameters (Fig. 2, fig. S2), the expression levels of perforin and granzyme B in CD56^bright NK cells correlated with IL-6 levels in the patients (Fig. 5C). Perforin in CD56^bright NK cells in the COVID-19 patients also positively correlated with Sequential Organ Failure Assessment (SOFA) and SOFA-R scores as well as neutrophil count and negatively with PaO2/FiO2 ratio (Fig. 5D). Similarly, patients with ongoing SARS-CoV-2 viremia in serum also presented with higher CD56^bright NK cell perforin expression, and these viremic patients had more severe disease (fig S7A). The up-regulation of perforin and granzyme B on CD56^bright NK cells (Fig. 2) was further positively associated with a general activation and up-regulation of effector molecules within CD56^bright and CD56^dim NK cells and inversely correlated with expression of the inhibitory check-point molecule TIGIT (Fig. 5E, fig. S7, B and C). As such, the association with high perforin was in line with the severe immunotype identified through unsupervised high-dimensional analysis (Fig. 4).

Thus, integration of clinical laboratory parameters into phenotypical NK cell analysis revealed that arming of CD56^bright cells with cytotoxic molecules correlated with parameters associated with a severe disease course.

Protein-protein interaction network driving CD56^bright NK cell arming in COVID-19

To better understand the significance of the activated CD56^bright NK cell phenotype and its association with disease severity, we took advantage of an extensive characterization of soluble serum proteins that had been performed on the same 27 patients through the Karolinska COVID-19 Immune Atlas effort. Of the 276 quantified soluble factors, the concentrations of 58 were significantly associated with perforin and granzyme B expression on CD56^bright cells and 9 showed strong, primarily positive, correlations (fig. S7D). To gain insight into the biological basis of the peripheral protein signature that correlated with CD56^bright NK cell activation, we analyzed protein-protein interaction networks. We identified 337 interactions (table S9) that were integrated into a network of 959 elements (nodes, Fig. 5F, table S10). Many of these related to viral infections, suggesting redundancy between different viral infections in the NK cell activation mode (Fig. 5G, table S11). Additionally, several distinct signaling pathways potentially involved in driving the CD56^bright NK cell response in COVID-19 together with the soluble factors were identified (Fig. 5H, fig. S7F, table S11). This included NF-κB, PI3K-Akt, and FoxO signaling (Fig. 5H). Altogether, our findings indicate that NK activation pathways could be shared with other viral infections, connecting NK cell-specific phenotypic traits to components of the systemic inflammatory milieu related to COVID-19 disease pathogenesis.