Diversity of peripheral blood human NK cells identified by single-cell RNA sequencing

Natural killer (NK) cells are cytotoxic lymphocytes that are crucial for viral defense¹  as well as detecting and killing cancer cells. This cell type is therefore of great interest for cancer immunotherapy.²  NK cells have also been implicated in a wide range of diseases: those caused by microbial pathogens,^3-5  autoimmune diseases,^6-9  reproductive complications,¹⁰  and transplantation.¹¹  Understanding how NK cells are involved in these diseases will be informed by better understanding of their different subsets.

Human peripheral blood NK cells, based on surface expression of CD56 and lack of CD3, have long been classified into 2 developmentally related but functionally distinct subsets: CD56^dim and CD56^bright. CD56^dim cells, the largest circulating population, express high levels of activating receptor CD16, lytic granules, and killer-cell immunoglobulin-like receptors (KIRs), called cytolytic cells. CD56^bright cells, sometimes termed cytokine-producing cells because they secrete high levels of cytokines and undergo robust proliferation after cytokine stimulation,¹²  express high levels of NKG2A, CD2, CD62L, and CCR7. An intermediate phenotype between CD56^bright and CD56^dim has also been proposed; these CD56^dimCD94^high cells express intermediate levels of CD62L, CD2, and KIR.^13,14 

Several other phenotypically and functionally distinct subsets have been described. CD56^neg cells, found at low frequencies within healthy individuals, are expanded in chronic^15,16  and acute viral infections.¹⁷  These cells display severely impaired natural cytotoxicity and antibody-dependent cellular cytotoxicity.¹⁸  Human cytomegalovirus (CMV) induces the expression of NKG2C¹⁹  and CD57^20-22  on NK cells, thereby skewing the NK cell repertoire and contributing toward age-associated changes.^23,24  These cells express a high frequency of self-KIR and exhibit a highly differentiated phenotype (CD56^dimCD16⁺LILRB1⁺KIR⁺NKG2A⁻) with hallmarks of adaptive immunity; they are therefore commonly termed adaptive or memory NK cells.^25,26  In the absence of CMV, CD57 expression marks “mature” NK cells in most human studies, because historical antigenic stimuli are rarely identified. Cytokine-induced memory-like (CIML) NK cells are a recently described subset. First identified in mice²⁷  and later in humans,^28,29  these cells show enhanced responsiveness to restimulation, including increased interferon-γ (IFN-γ) and granzyme B production after preactivation with cytokines (interleukin [IL]-12, -15, and -18). CIML NK cells are now being tested in preclinical trials,³⁰  because of their antileukemic properties and in vivo persistence after adoptive transfer.

Thus, several subsets of blood NK cells have been identified in health and disease. By focusing on blood NK cells from healthy human donors, we set out to determine which subsets could be identified by single-cell RNA sequencing (scRNA-seq). Our analysis identified unexpected diversity, encompassing both known and novel types of blood NK cell.

A full description of experimental and analytical methods can be found in the supplemental Methods.

Peripheral blood was obtained with informed consent from 2 healthy donors of Caucasian ancestry (one 60-year-old man and one 60-year-old woman) who were recruited to the National Repository of Healthy Volunteers, The University of Manchester. Ethical approval was obtained from North West Centre for Research Ethics Committee (REC: 99/8/84). Before donation, CMV status was determined with an enzyme-linked immunosorbent assay (human anti-cytomegalovirus IgG ELISA kit, ab108724; Abcam) according to the manufacturer’s instructions. Peripheral blood was similarly acquired by venipuncture from a CMV⁻ female 40-year-old donor from Texas Children’s Hospital (Protocol H-30487) and from 2 anonymous CMV seropositive donors (New York Blood Center, New York, NY).

Peripheral blood mononuclear cells were isolated from fresh blood collected into EDTA-coated tubes (Fisher Scientific) by density gradient centrifugation (Ficoll-Paque Plus; Amersham Pharmacia Biotech) or from leukocyte reduction system chambers through Ficoll density centrifugation. Leukocyte reduction system chambers were purchased from the New York Blood Center. Primary human NK cells were isolated by negative magnetic bead selection (Miltenyi Biotec or Stemcell Technologies), into medium alone (Dulbecco’s modified Eagle medium, 30% F12 Ham, 10% human serum, 1% nonessential amino acids, 1 mM sodium pyruvate [Sigma], 2 mM l-glutamine, 50 U/mL penicillin streptomycin, and 50 µM 2-mercaptoethanol [Gibco]), or into medium supplemented with 200 U/mL recombinant human IL-2 (Roche) for 4 hours at 37°C in 5% CO₂.

To assess purity, freshly isolated NK cells were washed in 1% fetal bovine serum/phosphate-buffered saline (PBS) and blocked with 1% human AB serum (Sigma-Aldrich) for 15 minutes at 4°C before staining. NK cells were stained for 20 minutes at 4°C with viability dye (Zombie NIR; Biolegend), phycoerythrin-labeled anti-CD3 monoclonal antibody (mAb; RRID:AB_314062), fluorescein isothiocyanate–labeled anti-CD14 mAb (RRID:AB_314186), AF647-labeled anti-CD19 mAb (RRID:AB_389335), and BV421-labeled anti-CD56 mAb (RRID:AB_11218798) or isotype-matched control mAbs (phycoerythrin-labeled IgG1k (RRID:AB_326435), fluorescein isothiocyanate–labeled IgG2a (RRID:AB_326458), AF647-labeled IgG1k (MOPC-21; Biolegend), and BV421-labeled IgG1k (MOPC-21; Biolegend). Finally, samples were washed in 1% fetal bovine serum/PBS, fixed in 2% paraformaldehyde/PBS for 10 min, assessed by flow cytometry (BD FACS Canto II flow cytometer, BD Biosciences), and then analyzed (FlowJo v 10 software).

Individually barcoded scRNA-seq libraries were prepared using the Chromium Single Cell Controller and Chromium Single Cell 3′ Library and Gel Bead Kit v2 or v1 kit (10X Genomics, Pleasanton, CA). All steps were performed according to the manufacturers’ instructions. Final libraries were sequenced on 1 flow cell of an Illumina HiSeq 4000 system (Illumina, San Diego, CA) with a read length of 26 bp for read 1 (cell barcode and unique molecule identifier [UMI]), 8-bp i7 index read (sample barcode), and 98 bp for read 2 (RNA read), to yield ∼1.25 billion reads per sample (208 000 reads per cell), or on the Illumina NextSeq500.

Raw FASTQ files were processed using CellRanger (v 2.0.1) to map reads against human genome 38 as a reference, filter out unexpressed genes, and count barcodes and UMIs. After alignment, more in-depth quality control, normalization, and analyses were conducted with Seurat (v 2.3.0)³¹  in R (v 3.3). Only data from high-quality cells were used for further analysis.

For comparative analysis, we used the Seurat alignment strategy³¹  to identify common cell types across the different donors and conditions. To do this, we calculated highly variable genes within each data set (supplemental Table 2), took the top 2000 per data set and then further selected 1501 genes, which were highly variable in at least 2 data sets; these 1501 genes were used in the canonical correlation analysis (CCA), to identify shared structures across data sets. Next, cells with expression profiles that were not explained by low-dimensional CCA vs low-dimensional principal components analysis (greater than twofold difference) were excluded from the data set and outputted for further investigation. The CCA subspace was then aligned and the cells were clustered by using a shared nearest neighbor modularity, optimization-based clustering algorithm³²  and visualized with t-distributed stochastic neighbor embedding (tSNE). Clustering was achieved by using 25 components from the CCA dimensionality reduction, identified by means of the shared correlation strength and using a resolution of 0.8.

To identify cluster-specific markers conserved across stimulation conditions, we determined differential expression for each cluster vs all remaining cells and within each cluster, and the P values were combined by the Wilkinson’s method from the metap R package. Differential expression within a cluster, determined by the Wilcoxon rank sum test, was used to identify markers between stimulatory conditions.

For comparison with published data, each cell was scored based on the expression of genes associated with known NK cell populations or signatures. In brief, mean expression values for each gene of interest within a signature were calculated for each cell, and the aggregated expression of control genes was subtracted, with the result indicating whether a gene’s expression was upregulated or downregulated vs the average. Genes within the signature were binned according to average expression levels, and the control genes were randomly selected from each bin.

Enrichment analysis was conducted by using GeneOntology (GO; http://www.geneontology.org/).^33,34  All conserved markers, either upregulated or downregulated within a cluster, with an adjusted P < .05 (in unstimulated cells) were selected for analysis. Enrichment scores (P values) for the selected GO terms were plotted as −log₁₀(P).

Transcriptional heterogeneity was investigated in blood NK cells from 2 CMV⁻ healthy donors (Figure 1A). NK cells from both donors were treated with IL-2 for 4 hours or left untreated, to gain insight into how different subsets respond to this commonly used cytokine stimulation. Confirming their expected cytolytic potential and cytokine production, cells degranulated in response to HLA-deficient target cells and less to HLA⁺ cells (Figure 1B) and produced IFN-γ (supplemental Figure 1).

In total, we sequenced 8905 cells (4992 unstimulated and 3913 IL-2 stimulated). On average, a median of 935 genes and 305 823 reads per unstimulated cell and 932 genes and 388 361 reads per stimulated cell were detected. Using CCA resulted in a loss of 257 cells that had a lower median gene count compared with the remaining cells (897 vs 940). This left 8462 individual cells (4747 unstimulated and 3715 stimulated) for further analysis.

Unsupervised clustering identified 10 subsets in total (Figure 1C). Cells from all 4 conditions contributed to each cluster (supplemental Figure 2A-B), suggesting that clustering was independent of donor-to-donor variability and stimulatory condition. Furthermore, clustering was not significantly influenced by the cell cycle (supplemental Figure 2C).

Differential expression analysis identified conserved markers per cluster, irrespective of stimulatory condition (supplemental Table 3). A small cluster, representing 1% of the total (cluster 9; Figure 1C), uniquely expressed IFI30, TNFAIP2, CLEC7A, and CD86 (supplemental Table 3), identifying them as antigen-presenting cells (APCs). Another small population (1.97% of the total) expressed CD3E, CD3D, and CD3G and major histocompatibility complex class II molecules (HLA-DPB1, HLA-DQB1, and HLA-DRB1), indicative of T cells or NKT cells (cluster 7; Figure 1C). Of note, CD3 staining by flow cytometry was negligible (Figure 1A), suggesting that this population may represent activated NKT cells that have downregulated the T-cell receptor complex.³⁵  In addition, a population (cluster 5; Figure 1C) representing 5.1% of the total expressed a disproportionate number of heat shock proteins and DNA repair genes, including HSPA1B, HSPA1A, and DNAJB1 (supplemental Table 3). These stressed cells are perhaps a result of the isolation process. These 3 populations (clusters 5, 7, and 9) were excluded from further analysis, leaving 7 NK cell clusters for investigation.

For each remaining cluster, the top 20 defining markers were visualized (Figure 1D). Apart from the 2 largest, denoted 0 and 1 (each having a similar core profile with some distinguishing features), the remaining clusters are clearly distinct. The top 2 markers associated with each cluster were further visualized (Figure 1E). The majority were expressed in a cluster-specific manner, with a few examples of markers shared between clusters. CREM, for example, was highly expressed in cluster 0, but expression was even higher in cluster 1. Together, these data identified 7 NK cell subsets.

Next, we related these clusters to known NK cell subsets, using canonical markers related to NK cell activation, inhibition, cytotoxicity, and cytokine secretion (Figure 2A). Clusters 0, 1, and 2 represent the largest populations, comprising 34.6%, 28.7%, and 12.2% of the total, respectively. They all share relatively high expression of FCGR3A (CD16) and cytolytic molecules, variable expression of activating/inhibitory receptors, and high levels of chemokines (Figure 2A). Of the assigned KIRs (KIR2DL1, KIR2DL3, KIR3DL1, KIR3DL2, KIR3DL3, and KIR2DL4), clusters 0, 1, and 2 possessed a high proportion of KIR⁺ cells (Figure 2B), which suggests that those clusters comprise the well-established bulk population of CD56^dimCD16⁺ NK cells.

Clusters 0 and 1 were closely related, evidenced by the high correlation of 0.99 (supplemental Figure 3). In a direct comparison of unstimulated cells (Figure 2C), cluster 0 was found to have significantly higher expression of chemokines (CCL4, CCL3, CCL4L2, and CCL3L3), whereas cluster 1 expressed slightly more TNFRSF9 (4-1BB). Thus, these 2 clusters of CD56^dimCD16⁺ NK cells most likely represent populations of subtly different maturity, rather than distinct cellular subsets.

There was also little to differentiate how these 2 clusters responded to IL-2 stimulation (Figure 2D). Many genes were upregulated, which would be expected with cytokine treatment, including cell-cycle markers (CCND1 and CDK6), IFNG and IL2RA, and regulators of cytokine signaling (SOCS1). In terms of cluster-specific changes, cluster 0 showed increased expression of several genes. Cluster 1 also displayed upregulation, but the fold-change was just below our labeling threshold. GO analysis highlighted enrichment in chemotaxis for cluster 0 and cell activation; whereas cluster 1 was enriched for a broader range of biological processes (Figure 2E). Altogether the data support that clusters 0 and 1 were very similar, each with a subtly distinct CD56^dimCD16⁺ phenotype, such that cluster 1 exhibited a more activated phenotype.

Compared with clusters 0 and 1, cells in cluster 2 expressed more lytic granule genes, FCGR3A (CD16; Figure 2F) and a higher proportion of KIR⁺ cells (Figure 2B), which suggests that these cells are terminally differentiated CD57⁺ NK cells that are highly cytotoxic, express CD57, frequently express inhibitory KIR, and restrict proliferative activity to cytokines.^36,37  Another defining marker is S100A4 (supplemental Table 3), which increases with NK cell maturation, with highest expression seen on CD57⁺ NK cells.³⁸  This cluster also highly expressed genes associated with cytoskeletal remodeling, such as ACTB, ARPC3/4, and CFL1, and negative cytolytic regulators, such as CST7 and PFN1, which are needed for optimal cytotoxicity³⁹  (Figure 2F). Although CD57 is underrepresented in our data set, strong evidence that these cells are terminally differentiated is that these cells do not respond to IL-2 stimulation^36,37  (Figure 2G). GO analysis further supported these observations, with enrichment for cytolysis and Fc receptor signaling, whereas IL-2–mediated signaling and NK cell differentiation were downregulated (Figure 2H). Thus, cluster 2 represents terminally differentiated CD56^dimCD16⁺CD57⁺ NK cells.

Cluster 4 comprised 5.5% of the total and had a distinct transcriptional profile at an extreme end of the spectrum, revealed in diffusion map analysis (supplemental Figure 4). This cluster exclusively expressed GPR183, IL7R, LTB, GZMK, CD62L, and CCR7 (supplemental Table 3). Furthermore, these cells expressed high levels of CD2 and KLRC1 (NKG2A), low levels of FCGR3A (CD16), and little or no PRDM1 (which decreases cytokine production and limits proliferation in CD56^dim cells)^40,41  (Figure 3A). This cluster also displayed a diminished cytolytic signature (Figure 2A). Together, these markers identify cluster 4 as CD56^bright NK cells.⁴² 

Human CD56^bright NK cells are transcriptionally similar to innate lymphoid cells (ILCs).⁴³  To determine whether ILCs contributed to cluster 4, negatively isolated NK cells were stained for CD56, CD127, and perforin in separate healthy donors. ILCs in healthy blood express CD127 (IL7R) and generally lack CD56 and perforin.⁴⁴  In contrast, here, most of the CD127⁺ cells stained positively for both CD56 and perforin (supplemental Figure 5A). Furthermore, CD127⁺ non-NK cell ILC markers reported by scRNA-seq⁴⁵  were expressed in a very small fraction of cells and were expressed across all clusters and at low levels of expression. Thus, this subset represents a specific subset of bona fide NK cells, rather than other ILCs (supplemental Figure 5B-C).

This cluster responded the most potently to IL-2 stimulation. Some of the genes influenced were those also effected in the 2 CD56^dimCD16⁺ populations (clusters 0 and 1) (Figure 3B). In addition, however, upregulation of β-actin (ACTB), lymphotoxin α/β (LTA/LTB), and TNFSF10 (TRAIL) was also observed. Transcripts for IL7R decreased (as previously shown⁴⁶ ), as well as CXCR4, ZFP36L2, and ATM. Consistent with the observation that cluster 4 represents CD56^bright NK cells, GO analysis indicated enrichment for response to cytokine and cytokine-mediated signaling pathway terms, while concurrently displaying a diminished cytolytic profile (Figure 3C).

Module scores calculated at a single-cell level using bulk expression profiles of sorted CD56^dim and CD56^bright NK cells⁴⁷  clearly showed that cluster 4 most resembled CD56^bright NK cells, whereas the remaining clusters resembled CD56^dim NK cells (Figure 3D: left and middle plots). Surprisingly, terminally differentiated NK cells (cluster 2) displayed features of CD56^bright NK cells (Figure 3D; middle plot), pinpointed to 5 markers (CTSW, DUSP1, FOS, JUN, and CD69), which were upregulated in a proportion of cells in cluster 2 (supplemental Figure 6). CD69, DUSP1, FOS, and JUN are either early activation markers or immediate early response genes (IEGs), and this analysis indicates that these markers overlap for activated versions of both CD56^dim and CD56^bright NK cells.

So far, this analysis has identified the more commonly described NK cell subsets: 2 populations of CD56^dimCD16⁺ cells, a population of terminally differentiated NK cells, and a CD56^brightCD16⁻ population, leaving 3 clusters (clusters 3, 6, and 8). Of these, cluster 3 was the largest, representing 5.9% of the total. These cells expressed less FCGR3A (CD16) and cytolytic markers than did the CD56^dim cells (clusters 0, 1, and 2; Figure 4A), but more than the CD56^bright cells (cluster 4), suggesting a diminished cytolytic capacity. Such a phenotype is attributed to nonclassic CD56^negCD16⁺CD7⁺ NK cells, which also lack SIGLEC7 expression, display reduced expression of KLRC1 (NKG2A) and natural cytotoxicity receptors, but still release substantial amounts of chemokines.¹⁸  In our data, cluster 3 expressed substantial amounts of CCL3 and CCL4 and significantly more CCL5 than all other clusters (Figure 4A).

Proteomic analysis of CD56^neg and CD56^dim blood NK cells from healthy donors⁴⁸  found marked similarities between the 2 (principal components analysis correlation 0.92), as also observed at the transcript level (supplemental Figure 3). Across 14 protein markers reported to have decreased expression on CD56^neg cells, we observed that 12 of these also showed a decrease in transcript abundance (Figure 4B). These CD56^neg cells respond to IL-2, similarly to how CD56^dim populations respond (Figure 4C).

Strikingly, the top 20 defining markers for CD56^neg cells (cluster 3), are almost all exclusively type I IFN-stimulated genes (Figure 4D). A strong IFN signature has been reported in another single-cell analysis of NK cells published recently.⁴⁹  Those authors termed these NK cells “inflamed” and found them in the bone marrow. Our data imply that such NK cells are also present in blood. Consistent with this finding, GO analysis indicates enrichment for pathways involved in defense response to viruses, response to type I IFN, and diminished cytolysis (Figure 4E). A subunit of the type I IFN receptor, IFNAR2, is also more highly expressed on a subset of cluster 3 compared with all other cells (supplemental Table 3). Upregulation of IFNAR can occur as result of IL-12 and polyIC stimulation.⁵⁰  Since CD56^neg cells (cluster 3) express IFNAR2 and large amounts of transcripts for CCL3, CCL4, and CCL5, it is also possible that these cells are particularly important for cross-talk with plasmacytoid dendritic cells, which express the ligand for CCL3, CCL4, and CCL5.

Cluster 6, representing 3.7% of the total, displays a phenotype that is something of a hybrid between CD56^dim and CD56^bright NK cells. Akin to CD56^bright cells, these cells express high levels of KLRC1 (NKG2A), IL2RA, XCL1, and XCL2 (Figure 5A). However, they also express FCGR3A (CD16) (to a lesser extent than CD56^dim cells), GZMA, PRF1, and GNLY (at levels similar to those expressed by CD56^dim cells); with even higher GZMB. This cluster also contains a high proportion of KIR⁺ cells (Figure 2B), which is reminiscent of CD56^dimCD94^high intermediary NK cells.¹⁴  However, cluster 6 did not express transcripts at a level to match previously identified protein markers¹⁴  (supplemental Figure 7).

The most prominent feature of this cluster is significant upregulation of GZMB after IL-2 treatment (Figure 5B). This is an important feature of enhanced-recall responses observed with in vitro differentiated CIML NK cells. Module score analysis using markers reported by time of flight cytometry, comparing CIML to control cells maintained in IL-15,²⁸  showed that cluster 6 most resembles CIML NK cells (Figure 5C). Upon IL-2 stimulation, several markers are upregulated that are common to other NK cell subsets (Figure 5D). Those specifically upregulated include CMTM6 (which positively regulates PD-L1 expression^51,52 ), RAB11FIB1, and SAMSN1. Although CIML NK cells are not new to the NK cell field, it is new to suggest that cells with features of CIML NK cells are constitutively present in circulating human blood.

GO analysis indicates enrichment of markers within a variety of pathways, including regulation of cell activation and cell communication and regulation of IL-10 and -12 production and granzyme B production (Figure 5E). Overall, this analysis demonstrated that cells within cluster 6 have a strongly activated phenotype and are primed to respond to IL-2 stimulation.

Cluster 8 makes up 1.6% of the total population and displays a phenotype not documented previously. It resembles terminally differentiated cells (cluster 2; Figure 6A; blue boxes) and both are remarkably similar according to correlation analysis (supplemental Figure 3). However, the definitive feature of this cluster is a striking loss of ribosomal expression (Figure 6A; red box, Figure 6B). Under nutrient starvation, ribosomes are targeted for autophagy-mediated degradation, as they are a major source of amino acids and nucleotides.^53,54  Here, we observed modest enrichment for autophagy-related markers within cluster 8 (Figure 6C). In particular, a proportion of cells showed increased expression of RB1CC1 (ATG17/FIP200), GABARAPL11 (ATG8), MAP1LC3b (LC3B), SH3GLB1 (BIF-1), BECN1 (BECLIN 1), and SQSTM1 (p62) (Figure 6D). The number of genes expressed (nGene) and number of unique molecular identifiers (nUMIs) were significantly lower in cluster 8 (supplemental Figure 8A). This is reminiscent of cells undergoing senescence or quiescence, but we did not observe any enrichment for genes specifically defining these processes (supplemental Figure 8B). The relative percentage of mitochondrial RNA did not differ across any of the subsets (supplemental Figure 8C; left) and module score analysis using proapoptotic markers were similar across all NK cell subsets, showing that these cells were not simply dying (supplemental Figure 8C; right). There were also no significant changes in this population upon IL-2 stimulation, but this may be related to insufficient power, given the small number of cells (supplemental Figure 8D). GO analysis indicated downregulation of oxidative phosphorylation (Figure 6E). Upon lymphocyte activation, cells switch from oxidative phosphorylation to glycolysis.⁵⁵  Interestingly, a proportion of cells in cluster 8 overexpressed the same IEG markers as seen in cluster 2, indicating that cells had been recently activated (supplemental Figure 6), which suggests that low ribosomal expression is linked to increased energy demands because of recent cellular activation.

A recent scRNA-seq analysis of human blood NK cells identified only 2 clusters.⁵⁶  In that study, analysis was performed per ∼1000 cells. Although there are differences in read depths between the previous study and this one, the number of cells has been shown to affect clustering to a greater degree.⁵⁷  We therefore sought to determine whether the extra clusters observed in this study were related to an increase in the number of cells. To achieve this, the analysis was limited to 1000 random cells. Following unsupervised clustering using 15 principal components and a low resolution of 0.4 (as used in the previous study⁵⁶ ), we observed 4 clusters (Figure 6F). The smallest, cluster 3, comprised the contaminating APCs, leaving 3 NK cell clusters. We overlaid the original cluster identities assigned using the full data set and found that the contaminating NKT cells (cluster 7), the cells that resemble CIML cells (cluster 6), and the type I IFN–responding cells (cluster 3) clustered with the 2 CD56^dimCD16⁺CD57⁻ populations (clusters 0 and 1) (Figure 6G), whereas the cells highly expressing heat shock proteins (cluster 5) clustered with the terminally differentiated CD57⁺ population (cluster 2). Overall, this analysis shows that the number of cells profiled is an important consideration of single-cell studies and establishes that our data are consistent with those of a previous work.⁵⁶  This study revealed additional NK cell subsets because a greater number of cells were analyzed.

Having defined 7 populations of NK cells in healthy human blood, we sought to use their transcriptional signatures to understand NK cell diversity in the context of latent CMV infection. Data obtained on blood NK cells isolated from 2 CMV⁺ but otherwise healthy donors and 3 CMV⁻ healthy donors were combined, looking for populations common across donors and limiting the study to 8000 random cells for comparability to the previous analyses. After unsupervised clustering, using 20 components from the CCA dimensionality reduction and a resolution of 0.7 resulted in 11 clusters (Figure 7A) that included a small population of APCs and a cluster of stressed cells that were removed, leaving 9 NK cell subsets for analysis.

We observed a very small, novel population of cells (Figure 7A; cluster J) which, upon further investigation, came primarily from 1 CMV⁻ donor. These cells resembled CD56^bright NK cells (supplemental Figure 9A-C), but lacked expression of KLRC1 (NKG2A) and some other CD56^bright NK cell markers, suggesting that they were transitioning or intermediate CD56^bright NK cells.¹⁴  Strikingly, we identified a population of adaptive NK cells found at substantially increased frequency in CMV⁺ individuals (Figure 7B; cluster C). A proportion of this cluster displayed higher expression of adaptive NK cell markers: KLRC2 (NKG2C) and KIR genes (Figure 7C). This population also expressed lower levels of FCER1G (FсεR1γ), SH2D1B (EAT-2), and ZBTB16 (PLZF), as shown previously,⁵⁸  and high expression of cytolytic molecules (GZMB and PRF1) and CST7 (a negative cytolytic regulators), in line with their potent cytotoxic potential. In addition, we observed increased expression of CD3, CD52, and IL-32, which have also been recently observed in adaptive NK cells⁴⁹  (Figure 7C). GO highlighted enrichment for the regulation of IFN-γ production, again supporting the strong functionality attributed to adaptive NK cells (Figure 7D).

For other NK cell subsets, frequencies altered slightly between the viral states (particularly for the larger populations), most likely to account for the change in adaptive NK cells. Terminally differentiated NK cells were less frequent in this combined analysis, perhaps because of a proportional shift caused by the adaptive NK cells. Inflamed type I IFN NK cells, CIML NK cells, and the low-ribosomal cells (supplemental Figure 9D-F) were found in similar proportions in both CMV⁻ and CMV⁺ individuals (Figure 7B). Interestingly, across all clusters, there was an upregulation of type I IFN genes and IL-32 in the CMV⁺ donors; as shown for cluster A (Figure 7E). Of note, IFIT1, OAS1, IFI6, IFITM1, MX2, and IL-32 were increased in expression across all NK cell subsets in CMV⁺ donors compared with CMV⁻ donors. Thus, CMV impacts both the frequencies of different NK cell subsets and their transcriptional profiles.

scRNA-seq is a powerful technique for discovering transcriptional diversity. Toward unraveling NK cell diversity, scRNA-seq studies have been performed in human tonsil NK cells, in comparison with developmentally related ILCs,⁴⁵  and on matched spleen and blood human NK cells⁵⁶  or NK cells from human bone marrow and blood.⁴⁹  However, these studies either did not resolve subsets, only resolved the CD56^dim and CD56^bright subsets, or were focused on unhealthy individuals or were performed with a relatively small number of blood NK cells. Thus, our objective was to define baseline blood NK cell heterogeneity from healthy individuals in a large number of individual cells, including an assessment of how these populations differentially respond to IL-2 stimulation. By analyzing ∼8000 individual cells, we identified an additional 4 subsets to the commonly identified CD56^dim and CD56^bright populations, while also observing subdivision within the CD56^dimCD16⁺CD57⁻ compartment.

CD56^bright cells have a distinctive phenotype, with their robust capacity to respond to cytokines. Closely associated with this subset was another cytokine-sensitive population: CD56^dimNKG2A^high cells, which showed enhanced cytotoxicity to IL-2, reminiscent of CIML NK cells. Cytokine-activated responses generated in vivo are functionally distinct from in vitro CIML cells generated with exogenous cytokines.⁵⁹  Our study highlighted similarities and differences from in vitro generated CIML cells. For example, NKp80 was more highly expressed in our study, but is downregulated on in vitro CIML cells,²⁸  whereas very little NKp30 was observed. Nevertheless, our data support the existence of CIML NK cells, or closely related cells, in peripheral human blood.

At the other end of the spectrum were terminally differentiated CD57⁺ NK cells, distinguished by their minimal response to IL-2, expressing many transcripts relating to cytolysis and several surface receptors for membrane-bound ligands, such as NKp30 and KIR. This subset seems to be specialized in recognition of pathogens via specific antigens. Multiple subpopulations of CD57⁺ have been distinguished using time of flight cytometry.⁶⁰  Our results suggest that, whereas individual CD57⁺ NK cells each have a specific repertoire of activating and inhibitory receptors, transcriptionally these cells form 1 population.

A subset closely related to CD56^dim cells was defined, by its response, not to IL-2, but to type I IFNs. This subset may be specialized for communication with plasmacytoid dendritic cells and T cells through their secretion of type I IFNs and CCL3, CCL4, and CCL5. Recent work in cancer has shown that NK cell–mediated recruitment of dendritic cells to the tumor microenvironment, through the chemokine CCL5 among others, is associated with patient survival.⁶¹  It remains to be determined whether these IFN-activated cells are the key to switching the immune response from a tolerogenic response to a cytotoxic one. Overall, this subset looks to have a specialized role in immune responses.

Cluster 8 represents a type of NK cell not previously characterized, resembling terminally differentiated CD57⁺ cells, but with significantly reduced ribosomal expression. Reduced ribosomal expression is a feature of a certain type of autophagy. First described in yeast cells⁵³  and more recently in mammalian cells,⁵⁴  under nutrient starvation, ribosomes are targeted for degradation and recycling in a process termed “ribophagy.” GO analysis highlighted a reduction of oxidative phosphorylation, which, when coupled with expression of IEGs, suggests a switch in metabolic states to support cellular activation. This would be consistent with the notion that these cells degrade ribosomes as a source of biosynthetic material to support cellular activation and function.

We also report subdivision within the CD56^dimCD16⁺CD57⁻ population. There were higher levels of 4-1BB on one subset and greater amounts of chemokines and a higher proportion of KIR⁺ cells on the other. However, their responses to IL-2 were very similar, suggesting that these represent subtle subdivisions. Education may play a role in setting the transcriptome of blood NK cells, but in a recent study, Goodridge et al.⁶²  concluded that education was controlled at the protein level and that transcriptional differences between educated and noneducated cells could not be observed using sorted populations and microarray technology.

Our data showed that CMV altered the frequencies of particular NK cell subsets, as well as gene regulation within each subset. The most striking difference was an increased frequency of adaptive NK cells in CMV⁺ donors. A small proportion of adaptive NK cells were still found in CMV⁻ individuals. A lower proportion of the adaptive cells within CMV⁻ individuals were NKG2C⁺ and, on average, expressed less NKG2C compared with CMV⁺ individuals, consistent with adaptive NK cells arising from infections other than CMV, but with CMV being a particularly strong driver for NK cell adaptive immunity.

In summary, in this study NK cells within healthy individuals displayed an exquisite division of labor (summarized in Figure 7F). The expansion, contraction, or alteration, of these subsets in pathological settings will give insight into NK cell responses during disease progression.

The authors thank the Genomic Technologies Core Facility, Faculty of Biology, Medicine, and Health (FMBH), The University of Manchester for performing scRNA-seq; Gareth Howell and the Manchester Collaborative Centre for Inflammation Research (MCCIR) Flow Cytometry Facility; and Clair Gardiner and Matt Hepworth for advice.  The authors also acknowledge help from IT Services and the use of the Computational Shared Facility at The University of Manchester.

This work was supported by the Medical Research Council (MC_PC_15072) a Wellcome Trust investigator Award (110091), Versus Arthritis fellowship (21745), the NIHR Manchester Biomedical Research Centre, Versus Arthritis grant 20385, the Manchester Academic Health Sciences Centre, and the MCCIR (funded by a precompetitive open-innovation award from GSK, AstraZeneca, and The University of Manchester), National Institutes of Health, National Institute on Aging grant AI130760, and an American Society for Hematology Junior Scholar Award (E.M.M.).

Contributions: All authors conceived the project and designed the experiments; S.L.S., P.R.K., J.D.W., K.B.S., E.M.M., B.M., S.S.C., A.H., R.S., M.L.S., and K.A. performed the experiments; S.L.S., P.R.K., P.M., E.M.M., S.S., E.H.S., A.H., R.S., M.L.S., K.A., and Y.-C.W. analyzed the data; S.L.S., P.R.K., and D.M.D. wrote the manuscript; and all authors critically appraised and edited the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

The current affiliation for S.S.C. is Qiagen, Houston, TX.

Correspondence: Daniel M. Davis, Manchester Collaborative Centre for Inflammation Research, University of Manchester, 46 Grafton St, Manchester M13 9NT, United Kingdom; e-mail: daniel.davis@manchester.ac.uk.