Cohort description

Our laboratory previously profiled 30 participants with COVID-19 from across the severity spectrum and 7 healthy controls by scRNA-seq as well as 26 participants (a subset of the 30 in scRNA-seq) and 12 healthy controls by CyTOF^29,30 (Fig. 1a). This cohort was well controlled for both viral variant and prior exposure to SARS-CoV-2, as these samples were collected during the first 4 months of the COVID-19 pandemic. Thus, all individuals experienced primary infection with ancestral SARS-CoV-2 because variants of concern had not yet evolved⁴. Peripheral blood mononuclear cells (PBMCs) and red blood cell-lysed whole blood were profiled using scRNA-seq, and PBMCs and NK cells were profiled using CyTOF (Supplementary Table 3). This cohort contained participants with COVID-19 from across the full range of the World Health Organization (WHO) disease severity spectrum, including participants with mild disease who remained unhospitalized (WHO 1–3), hospitalized participants with moderate disease (WHO 4–5), hospitalized participants with severe disease who required intubation (WHO 6–7) and participants who were later deceased (WHO 8). Samples were collected 0–66 days after a positive COVID-19 polymerase chain reaction (PCR) test (median 4 days), with the majority of samples drawn during the acute phase of disease (Extended Data Fig. 1a).

a, Study outline summarizing number of samples profiled and peak WHO severity score for each modality. VOC, variant of concern b, Boxplot of NT50 against five SARS-CoV-2 variant pseudoviruses. Each point represents the average of three technical replicates. Boxplots are drawn as median (center line), interquartile range (IQR; box) and 1.5× IQR (whiskers). P values calculated using two-sided, paired Wilcoxon signed-rank test with Bonferroni’s correction for multiple hypothesis testing. Lines connect unique donors. c, Heatmap of NT50 against five variants for each participant (column). Each cell represents the average of three technical replicates. d, UpSet plot of all unique combinations of SARS-CoV-2 variants neutralized, colored by peak WHO severity score. e, Scatter plot of breadth score and WHO severity score. f, Scatter plot breadth score and days after positive test. P value for e and f calculated using two-sided Spearman rank correlation. Each point represents one donor and is shaped by sample acuity. g, UMAPs of complete scRNA-seq (top) and CyTOF (bottom) datasets colored by cell types, breadth groups and peak WHO severity score. PB, plasmablast; Eos, eosinophil; Prog, progenitor; Prolif lymph, proliferating lymphocyte. Panel a created with BioRender.com.

Source data

Variability in antibody breadth during ancestral SARS-CoV-2 infection is independent of severity

To determine immunological correlates of antibody breadth in SARS-CoV-2 infection, we evaluated antibody neutralization breadth against variants of concern using matched plasma from the Stanford COVID-19 Biobank. Matched plasma samples were from the same or closest blood draw (one donor’s closest plasma sample was +3 days and the remainder were from the same draw). All individuals were infected with ancestral SARS-CoV-2, equalizing our measurement of breadth against variants of concern. To stratify participants by breadth against variants of concern, we measured serum neutralization against Wu-1 (ancestral) and four variant SARS-CoV-2 pseudoviruses that emerged after patient samples were collected (alpha, beta, delta and original omicron B.1.1.529) in a HeLa cell line stably expressing ACE2 and TMPRSS2³¹ (Fig. 1a–d). We found that there was heterogeneity in both the neutralization titer against each variant and in the number of variants neutralized (Fig. 1b–d). The median 50% neutralization titer (NT50; the reciprocal of the highest dilution to obtain <50% infection) against Wu-1 was the highest with ~1 log of variation in NT50 across individuals. The median NT50 of the cohort was reduced compared with Wu-1 in alpha, beta, delta and omicron B.1.1.529; this was significant in beta, delta and omicron (Fig. 1b). As our cohort showed little neutralizing activity against original omicron, as has been previously reported⁶, we did not include variants that emerged later. We defined each participant’s breadth score as the number of variants neutralized; breadth scores ranged from 0 to 4 (Fig. 1c). All individuals with a breadth score of 4 neutralized Wu-1, alpha, beta and delta; this neutralization profile represents the maximum breadth observed in our cohort of ancestral infection. To compare optimal breadth with all other individuals, those with a breadth score of 4 were classified as broad neutralizers, whereas those that neutralized different combinations of 0–3 variants were classified as narrow neutralizers (Fig. 1c,d).

As it is possible that severity could drive breadth via immune activation, we validated that the breadth score did not correlate with the WHO severity score (Fig. 1e). The severity score was also well distributed over different variant neutralization combinations (Fig. 1d). Therefore, we can be confident that immunological signatures of breadth cannot be primarily attributed to severity, and, by including participants across the entire severity spectrum, we can analyze immune correlates of antibody breadth regardless of immunological signatures of disease severity. In addition, the number of days between a positive test and sample collection did not correlate with breadth scores; therefore, breadth cannot be primarily attributed to acute or convalescent sample identity (Fig. 1f). Breadth score and breadth groups were assigned to each cell. As shown in Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) projections of scRNA-seq and CyTOF data, major immune cell types are represented in both broad and narrow neutralizers, allowing us to thoroughly investigate systemic contributions to a productive immune response with regard to antibody breadth (Fig. 1g).

NK cells in narrow neutralizers of SARS-CoV-2 highly express ISGs and markers of activation

To identify the most robust immunological signal that is correlated with breadth, we used pseudobulk whole PBMC samples from each patient to identify genes that significantly correlated with breadth score using DESeq2 (ref. ³²). In whole PBMCs, we identified 25 genes whose expression significantly correlated with low breadth score, 24 of which were ISGs, verified by the interferome database including classical ISGs such as MX family, OAS family and IFI family genes³³ (Fig. 2a). We next examined which cell types in narrow and broad breadth groups expressed these genes and found that NK cells from narrow neutralizers were among the highest expressors of the genes correlated with a low neutralization score (Fig. 2b). We also found NK cells to be one of only two cell types with differences in cell type proportion between broad and narrow neutralizers: broad neutralizers have significantly fewer NK cells when compared with narrow neutralizers and healthy controls (Fig. 2c). This trend further applied to breadth score, where NK cells as a proportion of PBMCs are inversely correlated with breadth score in both scRNA-seq and CyTOF modalities (Fig. 2d). Plasmacytoid dendritic cells (pDCs) (<1% of PBMCs in most participants) were significantly more abundant in narrow neutralizers in scRNA-seq data (Extended Data Fig. 2b,c). This perturbation of NK cells between breadth groups and breadth scores motivated us to further analyze the role of the NK cell compartment in modulating antibody breadth.

a, Heatmap of normalized expression of genes significantly correlated with low breadth score using DEseq2. Genes (columns) are ordered by hierarchical clustering; each row represents one donor. b, Heatmap of normalized expression of genes in a for each breadth group and cell type, ordered by hierarchical clustering for rows. c, Boxplots of proportion of NK cells in breadth groups from scRNA-seq and CyTOF datasets. Boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers), colored by peak WHO severity score and shaped by acuity. P values calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. Each point represents one donor. d, Correlation between proportion of NK cells and breadth score in scRNA-seq and CyTOF datasets. Line of best fit and 95% confidence interval (CI) are shown. P values calculated using two-sided Spearman rank correlation. Each point represents one donor.

In UMAP space, there was separation of NK cells by both severity and breadth groups (Fig. 3a,b). We identified differentially expressed genes (DEGs) in NK cells between broad and narrow breadth groups (log₂fold change (FC) >0.25 and adjusted P value <0.05). NK cells from narrow neutralizers exhibit higher expression of multiple ISGs such as IFI44L, RSAD2, XAF1, IFIT3, MX1, IFIT1 and so on as well as perforin (PRF1), a marker of NK cell cytotoxicity (Fig. 3c and Supplementary Table 2). NK cells from narrow neutralizers also upregulated CX3CR1, a marker of migration to lymphoid tissues where NK cells could influence antibody responses and that is associated with cytotoxic function³⁴.

a, UMAP of NK cells from scRNA-seq dataset colored by breadth group. b, UMAP of NK cells from scRNA-seq dataset colored by peak WHO severity score. c, Volcano plot of DEGs between broad and narrow neutralizers in NK cells from scRNA-seq dataset. P values calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. d, Protein–protein interaction graph depicting minimal significant interaction graph of NK cell DEGs (FDR <1 × 10⁻⁶). e, Dot plot of KEGG pathways significantly enriched in narrow neutralizer NK cells, showing their expression in each cell type. Plot colored by FDR and sized by number of genes present in KEGG pathway, with rows ordered by hierarchical clustering. f, Heatmap of normalized expression of NK cell genes significantly correlated with low breadth score using DEseq2. Each column represents one donor. g, Boxplots quantifying arcsinh-transformed mean signal intensity (MSI) of markers from CyTOF dataset in NK cells. Boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers), colored by peak WHO severity score and shaped by acuity. P values calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. Each point represents one donor.

Source data

To visualize the interaction landscape and underlying regulators of DEGs in NK cells, we leveraged the Bionet package to find the highest-scoring subgraph (false discovery rate (FDR) <1 × 10⁻⁶) of DEGs from a background protein–protein interaction network derived from the human STRING database^35,36,37. The minimal significant network of DEGs from NK cells revealed a hub of highly differentially expressed and interacting ISGs (MX1, IFIT1, ISG15 and so on) in narrow neutralizers with high predicted functional node scores relative to other nodes (Fig. 3d). This confirms a high level of ISG-mediated activation in NK cells from narrow neutralizers. CX3CR1 and perforin (PRF1) were also included in the minimally significant graph, highlighting both the migratory and cytotoxic potential of NK cells from narrow neutralizers.

Using DEGs upregulated in narrow neutralizers across all cell types, we performed pathway analysis to identify Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in narrow neutralizers. NK cells from narrow neutralizers were enriched for multiple pathways related to viral infections and pathological inflammation. These pathways were also significantly enriched in multiple immune cell types such as T cells and classical (CD14) monocytes in narrow neutralizers (Fig. 3e). This indicates that narrow neutralizers have heightened inflammation and ISG expression across the immune system. We also repeated our pseudobulk correlation analysis using DESeq2 to find genes significantly correlated with breadth score in NK cells (Fig. 3f). We found that all genes significantly correlated with a low breadth score were ISGs (verified by interferome database)³³; the majority overlapped with scRNA-seq DEGs from breadth groups (Fig. 3c,f). On the protein level, NK cells from all participants with COVID-19 exhibited upregulation of human leukocyte antigen (HLA) and death receptors (DR4–5); however, narrow neutralizers expressed significantly greater levels of CD8A and LLT-1 (markers of cytotoxicity and activation, respectively) when compared with broad neutralizers^38,39 (Fig. 3g). This reveals specific activation markers present only in narrow neutralizer NK cells in the setting of other signals of disease-driven activation.

NK cells are immature and proliferating in broad neutralizers of SARS-CoV-2

Distinct from narrow neutralizers, in broad neutralizers, NK cells were significantly less abundant, but they express higher transcript levels of MKI67 (Ki-67), KLRC1 (NKG2A) and IL7R (CD127), indicating a proliferative, inhibitory and immature phenotype^40,41 (Figs. 3c and 4a). Differential expression analysis indicates upregulation of some genes involved in immune activation such as HLA-DRA but no evidence of canonical ISG activation (Fig. 3c). This was substantiated on the protein level, where NK cells from broad neutralizers expressed significantly higher levels of NKG2A/CD94. CD56 expression was elevated in broad neutralizers, also supporting an immature phenotype (Fig. 4b). HLA-DR was upregulated on the protein but not the transcript level. Similarly, KLRD1 (CD94) and NCAM1 (CD56) were not differentially expressed in RNA but can have poor correlation between transcript and protein expression⁴² (Extended Data Fig. 2a). In addition, HLA-E, the ligand for NKG2A/CD94 heterodimer, is not differentially expressed on any cell type between broad and narrow neutralizers at the protein level⁴³ (Extended Data Fig. 2b).

a, Boxplots quantifying expression of selected DEGs in NK cells from scRNA-seq dataset. b, Boxplots quantifying arcsinh-transformed MSI of selected markers in NK cells from CyTOF dataset. c, Boxplots quantifying arcsinh-transformed MSI of HLA-bw4 across cell types from CyTOF dataset. d, Boxplots quantifying average expression of KEGG modules in NK cells from scRNA-seq data. All boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers), colored by peak WHO severity score and shaped by acuity. P values calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. Each point represents one donor.

All cell types in broad neutralizers express higher levels of HLA-Bw4 (significant in all cell types except B cells), whereas expression of HLA-Bw6 is not different between broad and narrow neutralizers (Fig. 4c and Extended Data Fig. 2c). HLA-Bw4 and HLA-Bw6 represent mutually exclusive epitopes on class I HLA-B, and HLA-Bw4 but not HLA-Bw6 is directly recognized by NK cells and is correlated with protection from COVID-19 and control of human immunodeficiency virus⁴⁴. This may support the partial contributions of genetics to productive antibody responses to SARS-CoV-2. A KIR3DL1 + HLA-Bw4+ genotype is also associated with protection from severe COVID-19 (ref. ⁴⁵); however, neither KIR3DL1 nor any other KIR measured in CyTOF or detected in scRNA-seq was differentially expressed between any group (Extended Data Fig. 2d,e). Using genes significantly upregulated in NK cells from broad neutralizers, we performed pathway analysis using the KEGG pathway database⁴⁶. We found upregulation of multiple pathways indicating normal cell functions and immaturity in NK cells from broad neutralizers (Fig. 4d and Supplementary Table 2).

NK cell transcriptomic clusters enriched in narrow neutralizers express cytotoxic proteins and differentiation markers

To further investigate the phenotype of NK cells enriched in broad and narrow neutralizers, we applied unsupervised clustering to scRNA-seq data of all NK cells from participants with COVID-19 (Fig. 5a). We found five biologically relevant clusters, of which C0 was significantly enriched in narrow and C1 in broad neutralizers (Fig. 5b). NK cells in C0 are CD56^dimCD16^high mature NK cells with higher co-expression of B3GAT1 (CD57) and KLRC2 (NKG2C) than in other clusters. C0 also expresses high transcript levels of KLRK1 (NKG2D), NCR3 (NKp30), FASLG (Fas-L), LAMP1 (CD107a), PRF1 (perforin) and GZMB (granzyme B), indicating cytotoxic potential, and had intermediate expression of CD69 and high levels of CD38, demonstrating activation. Owing to the expression of B3GAT1 and KLRC2, we investigated C0 and other clusters’ expression levels of receptors and transcription factors involved in the function of human cytomegalovirus (HCMV)-associated adaptive NK cells, which have been shown to lack cytotoxic function against autologous T cells⁴⁷. Distinct from HCMV-associated adaptive NK cells, NK cells in C0 express intermediate levels of SYK and high levels of FCGR3A (CD16), SH2D1B (EAT-2) and ZBTB16 (PZLF), all of which have low or absent expression in HCMV-associated adaptive NK cells. They also expressed the highest levels of NCR3 (NKp30), which is lacking on HCMV-associated adaptive NK cells. Thus, the C0 cluster enriched in narrow neutralizers is not representative of HCMV-associated adaptive NK cells but appears to have a phenotype consistent with high cytotoxic function and maturity. NK cells in the C1 cluster (enriched in broad neutralizers) are CD56^dimCD16^low but exhibit lower expression of activation markers CD38 and CD69 and cytotoxic markers KLRK1, NCR3, LAMP1 and PRF1. However, the C1 cluster still highly expresses some cytotoxicity-related genes (FASLG and GZMB) and the cytokine IFNG (IFNγ). Thus, the cluster of NK cells enriched in broad neutralizers represents NK cells with diminished activation and more potential for cytokine secretion. The C1 cluster also expressed high levels of MKI67 (Ki-67), further validating the increased proliferation of NK cells in broad neutralizers as compared with narrow neutralizers. Similar to HCMV-associated adaptive NK cells, C1 NK cells expressed low levels of NCR3, FCGR3A, SYK, SH2D1B and ZBTB16, while not expressing B3GAT1 or KLRC2. C1 NK cells illustrate that lack of NKp30 along with low expression of SYK, EAT-2 and PZLF may prevent killing of autologous cells in nonadaptive NK cell populations. We did not observe a clear population in C0–C4 of bona fide HCMV-associated adaptive NK cells, potentially owing to the low prevalence of HCMV in Santa Clara county where these samples originated^48,49. The remaining clusters, C2–C4, are not differentially expressed between broad and narrow neutralizers. C2 consists of CD56^dimCD16^high NK cells expressing low levels of activation markers. C3 comprises CD56^bright NK cells expressing cytotoxic NK cell receptors but not markers of degranulation. NK cells in C4 are CD56^dimCD16^high and express markers of degranulation but low levels of cytotoxic receptors (Fig. 5c).

a, UMAP of NK cells from scRNA-seq data colored by Seurat cluster. b, Boxplot NK cell cluster frequencies in each breadth group colored by peak WHO severity score and shaped by acuity. P value calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. Each point represents one donor. c, Heatmap of normalized gene expression of NK cell phenotypic and functional markers for each cluster. d, Violin plot of NK cell exhaustion (defined as average expression of LAG3, PDCD1 and HAVCR2; Methods) for each cluster. e, UMAP of NK cells from scRNA-seq data colored by relative differentiation calculated by cytoTRACE2. f, Boxplots of potency score of each NK cell cluster calculated by cytoTRACE2 ordered by potency. All boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers). Diff., differentiation.

In the setting of severe COVID-19, NK cells can become functionally exhausted while still expressing cytotoxic molecules^{50,51,52,53,54}. Thus, we evaluated whether the highly activated NK cells in narrow neutralizers might be exhausted. We analyzed aggregated expression of canonical exhaustion markers (LAG3, PDCD1 and HAVCR2) to generate an exhaustion score. Although all clusters exhibited some evidence of exhausted cells, the C0 cluster enriched in narrow neutralizers was not predominantly composed of cells expressing exhaustion genes (Fig. 5d). Thus, the expression of cytotoxic proteins observed in C0 is distinct from the exhaustion observed in severe COVID-19, in which NK cells do not retain intact cytotoxic function. Finally, we applied CytoTRACE2, a predictive method to evaluate absolute developmental potential in scRNA-seq data, to all NK cells from participants with COVID-19 (ref. ⁵⁵) (Fig. 5e). We found that C1 contained the most immature cells with the highest developmental potential, while C0 was the most differentiated, consistent with their expression of CD57 and NKG2C and further supporting the observation that broad neutralizers are associated with immature NK cells and narrow neutralizers are associated with differentiated NK cells (Fig. 5f).

Cell–cell interaction predictions indicate differences in cell signaling between broad and narrow neutralizers

Given that multiple immune cell types in narrow neutralizers demonstrated inflammatory phenotypes, we next investigated how cell–cell interactions may contribute to differential NK cell phenotypes in broad and narrow neutralizers. Our laboratory has previously found that monocytes contribute to NK activation and dysfunction in severe COVID-19 by interacting with NK cells both through direct receptor–ligand interaction and via cytokine secretion⁵⁶. On the protein level, we found that monocytes in narrow neutralizers express significantly higher levels of LLT-1 compared with broad neutralizers (Fig. 6a). LLT-1 can activate NK cells and is downregulated in moderate and severe COVID-19 (ref. ³⁰). LLT-1 binds CD161 on NK cells that is not differentially expressed^57,58 (Fig. 6b). ULBPs 1, 2, 5 and 6—ligands for the NK cell activating receptor NKG2D—were elevated in monocytes from broad neutralizers, while CD112 and CD155 (ligands for DNAM-1, TIGIT and/or CD96 on NK cells) were not differentially expressed between breadth groups^59,60,61,62 (Fig. 6a). Cognate NK cell receptors NKG2D, DNAM-1, TIGIT and CD96 were also not differentially expressed between breadth groups on NK cells (Fig. 6b).

a, Boxplots quantifying arcsinh-transformed MSI of NK cell ligands on monocytes in CyTOF data. b, Boxplots quantifying arcsinh-transformed MSI of cognate NK cell receptors in CyTOF data. All boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers), colored by peak WHO severity score and shaped by acuity. P values calculated using two-sided Wilcoxon rank sum test with Bonferroni’s correction for multiple hypothesis testing. Each point represents one donor. c,d, Circos plot showing top 50 predicted cell–cell communication pairs sent to NK cells colored by sender cell in narrow (c) and broad (d) neutralizers. e, Correlation between B cell subset and NK cell frequency colored by breadth group and shaped by acuity. Line of best fit and 95% CI are shown. P values calculated using two-sided Spearman rank correlation. Each point represents one donor.

Source data

To predict cell–cell communication from scRNA-seq data, we used the Multinichenetr package to identify the top differentially expressed receptor–ligand interaction pairs between breadth groups. The MultiNicheNet method leverages the downstream signaling database NicheNet-v2 which integrates ligand–receptor, signaling and gene regulatory data into a network model that allows linkage of ligand–receptor interactions with active downstream signaling⁶³ (Fig. 6c,d). There were distinct patterns of predicted cell–cell communication received by NK cells in broad and narrow neutralizers. In narrow neutralizers, multiple cell types (CD14 monocytes, CD16 monocytes, pDCs, platelets and dendritic cells (DCs)) were predicted to signal via TGFBI (TGFβ induced protein) through integrins (ITGA5 and ITGA6), which inhibits adhesion and migration⁶⁴. We also investigated predicted downstream signaling targets of both TGFβI and TGFβ in the top 50 predicted cell–cell communication pairs in narrow neutralizers. TGFβI (TGFBI) had low regulator potential to drive ISG expression, but TGFβ (TGFB, upstream of TGFBI) is predicted to drive expression of many ISGs including MX1 and ISG15 in narrow but not broad neutralizers (Extended Data Fig. 3a,b). TGFB and IRF7 co-expression in the context of viral infection can drive type I interferon and ISG expression^65,66. Here, IRF7 was significantly upregulated in narrow neutralizers, which may work in conjunction with TGFβ signaling to drive interferon-stimulated gene expression in narrow neutralizers. Other sources of activation included CD48 interaction with 2B4 (CD244)⁶⁷ on CD16 monocytes and CD7-mediated activation by CD14 monocytes through SECTM1 (ref. ⁶⁸). NK cells themselves were also predicted to send auto-activation signals via SLAMF7–SLAMF7 (refs. ^69,70), CD56–CD56 (NCAM1)⁷¹ and CD8–HLA-F⁷² interactions. There were also predicted interactions between SPON2 and integrin alpha 5 (ITGA5) in NK cells from narrow neutralizers, which can inhibit migration⁷³. In addition, signals involved in the inhibition of angiogenesis (SEMA4A–PLXND1⁷⁴ and SERPINF1–PLXDC1⁷⁵) were predicted from pDCs and plasmablasts (Fig. 6c).

CD4⁺ T cells, CD8⁺ T cells, proliferating lymphocytes, platelets and B cells were only predicted to communicate with NK cells in broad neutralizers. Here, signaling was dominated by inhibitory receptors; HLA-E, HLA-A, HLA-F and HLA-G were predicted to signal in multiple cell types through NKG2A (KLRC1) and/or LILRB1^{40,43,76,77,78,79,80}. Receptor–ligand pairs involved in adhesion and lymphocyte homing were also predicted with signals from galectin 1 (LGALS1) from proliferating lymphocytes and CD16 monocytes; selectin P ligand (SELPLG) from NK cells, CD4⁺ T cells and CD8⁺ T cells; versican (VCAN) from monocytes; selectin P (SELP) from platelets and other signals through integrins across multiple cell types^81,82,83,84. Some of the integrin-mediated signaling predicted in broad neutralizers are known to be inhibitory, such as TIMP2 interaction with integrin beta 1 (ITGB1)⁸⁵. Platelets are predicted to send inhibitory signals in broad neutralizers via selectin P (SELP), HLA-A, HLA-E and THBS1. Additional relevant signaling included immune synapse proteins CD58 and CD2 as well as TIMP2–CD44 that can also be involved in migration and activation^86,87 (Fig. 6d). Overall, these findings suggest that signals sent to NK cells in narrow neutralizers drive activation, including ISG expression as well as inhibition of migration and adhesion, and broad neutralizers are characterized by dominant inhibitory signals as well as positive regulation of migration and immune synapse formation.

Although predicted cell–cell communication between NK cells and T_FH is of particular interest, we re-clustered CD4⁺ T cells and did not identify T_FH in our scRNA-seq dataset. This is expected, as T_FH are extremely rare in blood^88,89 and too few cells are evaluated by scRNA-seq for T_FH to be well represented (Extended Data Fig. 3c,d). However, we hypothesized that we could capture the indirect effect of NK cells on T_FH by analyzing the correlation between NK cell and B cell prevalence. The proportion of memory B cells had a significant negative correlation with the proportion of NK cells (Fig. 6e). The proportion of intermediate B cells and plasmablasts also exhibited a negative correlation with the prevalence of NK cells, and this trend approached significance. Total B cells and naive B cells did not exhibit this same correlation. This suggests that NK cell activity is correlated with diminished frequency of the B cell subsets that are critical for antibody production and durable antibody-based memory.

Interferon α-activated NK cells exhibit enhanced killing and upregulate cytotoxic markers when cocultured with iT_FH-like cells

The predicted cell–cell communication in the peripheral immune system did not identify a source of interferons. We hypothesize that the source of interferon-mediated activation was in tissue. Using publicly available data of ancestral SARS-CoV-2 infection from nasopharyngeal swabs, we found that CD8⁺ T cells express IFNG in nasal tissue^90,91 (Extended Data Fig. 4a,b). IFNA transcripts were not detected in scRNA-seq data across all datasets but are thought to be upstream of IFNG expression in T cells^92,93,94. While we cannot look at differences in broad and narrow neutralizers in tissue expression of interferons, as no such data or samples exists, this shows that interferons are expressed in immune tissue during ancestral SARS-CoV-2 infection, which could disseminate to the blood, causing ISG expression observed in narrow neutralizers.

The activation and cytotoxic potential of NK cells in narrow neutralizers raises the hypothesis that NK cells may target T_FH in lymphoid tissue to limit antibody breadth. To directly assess the functionality of NK cells with narrow neutralizer phenotype against T_FH, we developed an in vitro coculture system with healthy human donor NK cells and iT_FH-like cells. Bona fide T_FH cells cannot be differentiated in vitro, but activation of primary, naive CD4⁺ T cells with TGFβ and interleukin (IL)-12 has been shown to induce high expression of CXCR5, ICOS and BCL6, which recapitulates key features of T_FH cells^95,96. After isolation of naive CD4⁺ T cells from cryopreserved, healthy PBMCs, cells were activated with staphylococcal enterotoxin B (SEB; a T cell superantigen) for 4 h; then TGFβ and IL-12 were added to the culture for 72 h and CXCR5+ cells were sorted from this population. Sorted cells are referred to here as iT_FH-like cells (Fig. 7a and Extended Data Fig. 4c). These iT_FH-like cells were confirmed to express high levels of ICOS, PD-1, CD40L and BCL6 when compared with resting (freshly isolated) CD4⁺ T cells or CD4⁺ T cells activated for 72 h in the absence of cytokines, indicating a T_FH-like phenotype (Extended Data Fig. 4d).

a, Experimental design of in vitro coculture experiments with healthy IFNα-activated NK cells cocultured with iT_FH-like cells. b, Boxplots of ΔCt relative to 18S for IFIT1, ISG15, MX1 and OAS1 in NK cells from healthy donors unstimulated or activated with IFNα for 24 h. Each point represents the average of three technical replicates (n = 10 across three experiments). c, Boxplots of percentage of NK cells expressing cytotoxic and activation markers after 16-h coculture +/− iT_FH-like cells (E:T = 1:5) and +/− NK activation with IFNα before coculture (n = 11 across three experiments). d, Boxplots of background-subtracted percentage of dead target cells in 3-h killing assay with iT_FH-like, activated or resting CD4⁺ T target cells +/− NK activation with IFNα before killing assay (E:T = 5:1) (n = 18 across seven experiments for iT_FH-like target cells, n = 10 across four experiments for activated and resting CD4 target cells). e, Boxplot of CXCR5 mean fluorescence intensity (MFI) on CD56^dim CD16^high NK cells +/− stimulation with IFNα for 24 h (n = 12 across two experiments). All boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers). P values calculated using two-sided, paired Wilcoxon signed-rank test with Bonferroni’s correction for multiple hypothesis testing. Lines connect unique donors. f, Correlation between Δkilling (% killing IFNα pre-activated NK cells − % killing unstimulated NK cells) and ΔΔCt for expression of ISGs relative to unstimulated NK cells. Line of best fit and 95% CI are shown. P values calculated using Spearman rank correlation. Panel a created with BioRender.com.

Source data

In order to mimic the high ISG expression in narrow neutralizer NK cells, we activated healthy, isolated NK cells with interferon α (IFNα) to induce ISG expression⁵¹. We found that in healthy donor NK cells, 24 h of activation with IFNα led to significantly greater expression (lower ΔCt relative to 18S) of hallmark ISGs: MX1, OAS1, ISG15 and IFIT1 measured by reverse transcription quantitative PCR (qPCR) (Fig. 7b). We then investigated how NK cells pre-activated with IFNα responded to iT_FH-like cells by coculturing them at an effector:target (E:T) ratio of 1:5 overnight (Fig. 7a). Following coculture, IFNα-activated NK cells expressed significantly more CD107a, perforin and granzyme B compared with unstimulated NK cells in coculture. Expression was also increased in coculture compared with IFNα-activated NK cells in the absence of target cells. This further indicates that NK cells respond specifically to iT_FH-like target cells rather than IFNα activation alone. This trend was also observed in the percentage of activated (CD38+CD69+) NK cells. IFNα-activated NK cells in coculture also expressed more IFNγ than IFNα-activated cells cultured alone, but this difference was not significant between IFNα-activated and unstimulated NK cells in coculture. The proportion of NKG2D+ NK cells was slightly decreased in IFNα-activated cocultured NK cells, although this difference was not significant (Fig. 7c and Extended Data Fig. 5a). When we split NK cells by subtypes (that is CD56^dimCD16^high, CD56^dimCD16^low and CD56^bright), there were no major differences in responses to IFNα or iT_FH-like target cells, with subsets following similar trends to bulk NK cells (Extended Data Fig. 5b–g).

We also performed killing assays with IFNα-activated NK cells against iT_FH-like cells at an E:T ratio of 5:1 for 3 h. We found that IFNα-activated NK cells killed significantly more iT_FH-like cells than unstimulated NK cells (median difference 2.6%) (Fig. 7d). While NK cells have been reported to kill activated CD4⁺ T cells, we determined that interferon-driven killing is iT_FH-like specific because IFNα-activated NK cells did not kill significantly more activated or resting CD4⁺ T cells. We further investigated the specificity of interferon-driven response to iT_FH-like cells by repeating coculture experiments overnight at E:T 1:5 with activated CD4⁺ T cells. We similarly observed that NK cell functional responses to activated CD4⁺ T cells were not IFNα dependent; for example, some markers such as perforin were upregulated by both unstimulated and IFNα-stimulated NK cells in response to activated CD4⁺ T cells and other markers such as CD38 were upregulated in response to IFNα but in both target cell and NK-only wells (Extended Data Fig. 6a,b). This shows that the IFNα-driven NK cell functional and cytotoxic response is specific to iT_FH-like target cells and not target cell activation.

We also found evidence that IFNα activation and ISG expression may increase NK cell migratory potential to lymphoid tissues where T_FH reside. In CD56^dimCD16^high cytotoxic NK cells, treatment with IFNα increased CXCR5 but not CX3CR1 expression, which could permit NK migration to lymph nodes (Fig. 7e). This in vitro IFNα activation did not recapitulate the higher expression of CX3CR1 in narrow neutralizers, possibly because in vivo activation signals were not fully captured in our system (Extended Data Fig. 6c,d). CXCR5 expression decreased in CD56^bright NK cells after IFNα activation in vitro. This subset is the most immature, similar to the NK cell phenotype in broad neutralizers, indicating that the NK cell phenotype enriched in broad neutralizers has lower migratory potential (Extended Data Fig. 6d). Finally, we found that the magnitude of upregulation of ISGs after IFNα activation (the absolute value of the difference in ΔCt values between IFNα-activated and unstimulated NK cells) correlated with the difference in killing of iT_FH-like cells between IFNα-activated and unstimulated NK cells (ΔKilling), with greater cytotoxicity associated with greater increase in ISG gene expression. This correlation was only statistically significant for OAS1, and other genes approached significance (Fig. 7e). Overall, these data suggest that NK cells with narrow neutralizer phenotype have the capacity to specifically target autologous T_FH driven by ISG expression.

NKG2D and NKp30 contribute to NK cell killing of iT_FH-like cells

To understand the mechanism of NK cell targeting of iT_FH-like cells, we profiled NK cell ligand expression on iT_FH-like cells compared with resting and activated CD4⁺ T cells. Ligand expression was also assessed on CXCR5+ (iT_FH-Like) and CXCR5− (not iT_FH-like) populations derived from the unsorted pool of differentiated naive CD4⁺ T cells (Extended Data Fig. 6e). We quantified expression of B7-H6, ULBP1, ULBP2-5-6, ULBP3, MICA, MICB and CD48, ligands of activating receptors NKG2D, NKp30 and 2B4 (refs. ^{59,60,97,98,99}) (Fig. 8a and Extended Data Fig. 6f). We found increased expression of B7-H6, ULBP1, ULBP2-5-6, ULBP3, MICA and MICB on iT_FH-like cells compared with activated and resting cells. Although activation increased the expression of most ligands on CD4⁺ T cells compared with resting cells, iT_FH-like cells showed even higher levels of B7-H6, ULBPs and MICA/B than those seen in activated CD4⁺ T cells or CXCR5− cells. CD48 was highly expressed on all CD4⁺ T cell types (Fig. 8a). To determine if NKp30 (binds B7-H6) or NKG2D (binds ULBPs and MICA/B) contribute to killing of iT_FH-like cells, we blocked these receptors during killing assays against iT_FH-Like cells (Fig. 8b). Blocking NKG2D or NKp30 alone did not significantly reduce IFNα-driven killing of iT_FH-like cells. However, when we blocked both NKG2D and NKp30, IFNα-activated NK cells killed significantly less iT_FH-like cells, comparable to the level of killing of unstimulated NK cells (Fig. 8c). This suggests that both NKG2D and NKp30 play a role in interferon-driven NK cell killing of T_FH and suppression of antibody responses. Notably, NK cell transcriptomic clusters enriched in narrow neutralizers expressed higher levels of NKG2D and NKp30, and IFNα-activated NK cells downregulated NKG2D in culture with iT_FH-like cells, which can indicate use of this receptor that is internalized after ligand engagement¹⁰⁰, further supporting a role for these activating receptors in NK-cell-mediated suppression of T_FH in our cohort (Figs. 5c and 7c).

a, Boxplots of percentage cells expressing NK ligands on iT_FH-like (CXCR5+), CXCR5−, activated and resting CD4⁺ T cells. b, Diagram of NK cell ligands profiled on iT_FH-like cells and their cognate NK cell receptors. c, Boxplots of background-subtracted percentage of dead target cells in 3-h killing assay +/− NKp30- and/or NKG2D-blocking antibodies with iT_FH-like target cells +/− NK activation with IFNα before killing assay (E:T = 5:1) (n = 8 across two experiments). d, Experimental design of in vitro coculture system with healthy IFNα-activated NK cells, iT_FH-like cells and B cells with analysis modalities. e, Boxplots quantifying cell numbers after 6 days of coculture of iT_FH-like and B cells with NK cells +/− stimulation with IFNα (n = 8 across two experiments). f, Boxplots quantifying immunoglobulin concentration in cell culture supernatants after 6 days of coculture of iT_FH-like and B cells with NK cells +/− stimulation with IFNα (n = 8 across two experiments). All boxplots are drawn as median (center line), IQR (box) and 1.5× IQR (whiskers). P values calculated using two-sided, paired Wilcoxon signed-rank test with Bonferroni’s correction for multiple hypothesis testing. Lines connect unique donors. Panels b and d created with BioRender.com.

Source data

NK cell restriction of iT_FH-like cells directly suppresses B cell responses in in vitro cocultures

Finally, to determine whether interferon-driven activation of NK cells can modulate antibody responses and B cell function via restriction of T_FH cells, we used an in vitro coculture system with NK, B and iT_FH-like cells. Unstimulated or IFNα-activated NK cells were incubated with iT_FH-like cells for 1 h at an E:T of 5:1 before addition of autologous B cells. B, iT_FH-like and NK cells were cocultured for 6 days at a ratio of 2:1:5, and cell counts and supernatant antibody levels were evaluated using counting beads and enzyme-linked immunosorbent assay (ELISA), respectively (Fig. 8d). IFNα-activated NK cells suppressed survival and proliferation of iT_FH-like cells in this long-term coculture system. This confirms that low-level differences in cytotoxicity against iT_FH-like cells in short-term killing assays are amplified in long-term cocultures. After 6 days, there were also significantly fewer B cells surviving in coculture with IFNα-activated NK cells than with unstimulated NK cells. Moreover, the numbers of class-switched (IgM−) B cells and plasmablasts (CD38+CD27+) were also significantly reduced in coculture with IFNα-activated NK cells (Fig. 8e). The lower level of B cell survival and activation impacted antibody responses and secretion, with reduced levels of IgG and a significantly lower concentration of IgM in the supernatants from cocultures with IFNα-activated NK cells (Fig. 8f). Here, we confirm that direct effects of NK cell activity on iT_FH-like cells impact specific B cell subsets such as plasmablasts (validates Fig. 6e) and also affect antibody production in vitro.