High-dimensional single-cell analysis of human natural killer cell heterogeneity

Abstract

Natural killer (NK) cells are innate lymphoid cells (ILCs) contributing to immune responses to microbes and tumors. Historically, their classification hinged on a limited array of surface protein markers. Here, we used single-cell RNA sequencing (scRNA-seq) and cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) to dissect the heterogeneity of NK cells. We identified three prominent NK cell subsets in healthy human blood: NK1, NK2 and NK3, further differentiated into six distinct subgroups. Our findings delineate the molecular characteristics, key transcription factors, biological functions, metabolic traits and cytokine responses of each subgroup. These data also suggest two separate ontogenetic origins for NK cells, leading to divergent transcriptional trajectories. Furthermore, we analyzed the distribution of NK cell subsets in the lung, tonsils and intraepithelial lymphocytes isolated from healthy individuals and in 22 tumor types. This standardized terminology aims at fostering clarity and consistency in future research, thereby improving cross-study comparisons.

Fig. 1. CITE-seq analysis reveals three prominent subsets of peripheral blood NK cells in healthy individuals.

Based on dataset 5. a, WNN and UMAP (WNN_UMAP) visualization of NK cells sorted from healthy human blood with clusters identified by unsupervised hierarchical clustering (based on scRNA-seq and expression of 228 surface proteins). b, Dot plot of the 20 most distinguishing genes expressed for the three major subsets of human blood NK cells. Gene expression was analyzed using the using the two-sided Wilcoxon rank-sum test with Bonferroni adjustment. Ribosomal genes and mitochondrial genes were removed for clarity. The color indicates the Z-score scaled gene expression levels. c, Dot plot of the most distinguishing proteins expressed for the three major subsets of human blood NK cells. Protein expression was analyzed using the two-sided Wilcoxon rank-sum test with Bonferroni adjustment. Alternative protein names are shown in parentheses. The color indicates the Z-score scaled protein expression levels. d, WNN_UMAP visualization of the surface expression of the major discriminating proteins expressed at the surface of NK1, NK2 and NK3 cells. a.u., arbitrary units.

Fig. 2. The three most important NK cell populations can be subdivided into six subgroups.

Based on datasets 1–4a. a, UMAP visualization of NK cells sorted from healthy human blood, with clusters identified by unsupervised hierarchical clustering. b, Bar graph showing the proportion of cells within each cluster in all donors. Blue and pink bars are shown under HCMV-positive and HCMV-negative individuals, respectively. c, Violin plot of the scoring of the six NK clusters with respect to established NK1, NK2 and NK3 signatures (n = 13 samples). In the violin plots, the point is the median value. The error bars present the median +/- standard deviation. d, Violin plot of the scoring of the NK3 clusters with the NK3 signature in HCMV-positive and HCMV-negative individuals (n = 13 samples). e, Dot plot of the 20 most distinguishing genes for each subset of NK cells. Gene expression was analyzed using the two-sided Wilcoxon rank-sum test with Bonferroni adjustment. Ribosomal genes and mitochondrial genes were removed for clarity. The color indicates the Z-score scaled gene expression levels.

Fig. 3. Markers of interest, functions and metabolism characterizing NK cell populations.

Based on datasets 1–4a. a, Heatmap showing the differential expression of markers of interest among NK cell subsets. The color scale is based on z-score-scaled gene expression. The z-score distribution ranges from −2 (blue) to 2 (red). b,c, Selected GO terms showing enrichment in the three major populations (and six major subsets) of healthy human blood NK cells. Benjamini–Hochberg-corrected −log₁₀(P) values were calculated by a hypergeometric test. The black dotted line indicates the significance threshold, which is −log₁₀(0.05). d, Heatmap showing the differential enrichment for selected metabolic pathways among NK cell subsets. The color scale is based on the z score of the normalized enrichment score for each metabolic pathway. e, Assessment (z scores) of the response to different cytokines and chemokines in each subgroup, quantified by Cytosig. IFN, interferon; PGE2, prostaglandin E2. P values were computed by comparing z scores in one NK subset with those in other subsets using two-sided Student’s t-tests, and –log₁₀(P) values exceeding 10 were capped at 10 to facilitate visualization.

Fig. 4. Putative transcriptional trajectories connecting NK cell subpopulations.

a–e, Based on dataset 4a. f–g, Based on datasets 1–4a. a,b, RNA-velocity analysis and pseudotime inference based on velocity analysis in a representative sample. c,d, Confirmation of trajectories and pseudotime analysis using a diffusion-map approach (Destiny analysis of dataset 4a). e, Plot of pseudotime derived from diffusion-map analysis across all NK cell subgroups. f, Monocle-derived trajectory performed on the NK1A, NK1B and NK1C subsets and projected onto the UMAP, colored by clusters. g, Monocle-derived trajectory performed on NK1 subpopulations and projected onto the UMAP colored by pseudotime (inferred by Monocle3).

Fig. 5. Putative ontogeny of the main NK populations.

Based on datasets 1–4a. a, UMAP visualization of the module score of individual cells scored with signatures derived from the main NK cell progenitor (ENKP) identified in mice. b, Violin plots of the module scores of individual cells scored with ENKP signatures. Data are shown as median ± s.d. (n = 13 samples). b,d, In the violin plots the point is the median value. The error bars present the median +/- standard deviation. c, UMAP visualization of the module score of individual cells scored with human blood ILCP signatures. d, Violin plots of the module scores of individual cells scored with signatures of human blood ILCPs Data are shown as median ± s.d. (n = 13 samples).

Fig. 6. Distribution of NK1, NK2 and NK3 cell subsets in tissues.

Based on dataset 7. a, UMAP visualization of the main populations of group 1 ILCs present in PBMCs, tonsil, lung and IEL, colored by their main populations and by their cluster and tissues (as defined in dataset 7). b, UMAP visualization of the module score of individual cells scored with signatures of NK1, NK2 and NK3 of ILC populations present in tonsil, lung and IELs. iel_prdm1, Intestinal intraepithelial lymphocytes; PRDM1+ NK cells; lung_bright, lung CD56bright NK cells; lung_dim, lung CD56bright NK cells; lung_znf NK, lung ZNF683+ NK cells; pbmc_bright, blood CD56bright NK cells; pbmc_dim, blood CD56dim NK cells; tonsil_bright, tonsil CD56bright NK cells; tonsil_dim, tonsil CD56dim NK cells; tonsil_znfNK, tonsil ZNF683+ NK cells.

Fig. 7. Distribution of NK1, NK2 and NK3 cell subsets in the blood of people with cancer and at the tumor bed.

Based on datasets 1–4a and 6. a, Bar graph showing the proportion of the three main NK populations in the blood and at the tumor bed in 22 cancer types (n = 676 samples). b, PCA on tumor-infiltrating and blood NK cells, grouped by NK population, cancer conditions and tissue. The PCA is based on the mean expression levels of the 2,000 genes most differentially expressed across tissue and conditions. Groups are colored on the basis of their tissue of origin. PC1 and PC2 explained 13.7% and 11% of the variance, respectively (n = 676 samples). MELA, melanoma; MM, multiple myeloma; RC, renal carcinoma; FTC, fallopian tube carcinoma; CLL, chronic lymphocytic leukemia; ALL, acute lymphocytic leukemia; PACA, pancreatic carcinoma; THCA,: thyroid carcinoma; LC, lung cancer; HCC, hepatocellular carcinoma; PRAD, prostate cancer; GC, gastric cancer; CRC, colorectal cancer; ICC, intrahepatic cholangiocarcinoma; HNSCC, head and neck squamous cell carcinoma; OV, ovarian cancer; BRCA, breast cancer; UCEC, uterine corpus endometrial carcinoma; ESCA, esophageal cancer; NB, neuroblastoma; NPC, nasopharyngeal carcinoma; BCC, basal cell carcinoma.

Fig. 8. Distinct transcriptional phenotypes of NK1, NK2 and NK3 cell subsets in the blood of people with cancer and at the tumor bed.

Based on datasets 1–4a and 6. a, PCA of blood NK cells, grouped by NK population and cancer conditions. The PCA is based on the mean expression levels of the 2,000 genes most differentially expressed across tissue and conditions. Groups are colored on the basis of NK cell subsets. PC2 and PC3 explained 8.1% and 7% of the variance, respectively (n = 676 samples). b, PCA of tumor-infiltrating NK cells, grouped by NK population and cancer conditions. The PCA is based on the mean expression levels of the 2,000 genes most differentially expressed across tissues and conditions. Groups are colored on the basis of NK cell subsets. PC2 and PC3 explained 13.9% and 12.7% of the variance, respectively (n = 676 samples).

Extended Data Fig. 1. The classification of NK cells into 3 main families is robust in other blood NK cell atlas.

Based on Dataset 6. a, Uniform Manifold Approximation and Projection (UMAP) of blood NK cells from Tang et al. pan-cancer NK cells atlas. Subsets constituting less than 1% of circulating NK cells were excluded which resulted in a total of 84,343 human blood NK cells for analysis. b, UMAP of blood NK cells from Dataset 5 scored with NK1, NK2, and NK3 signatures. c, RidgePlot visualization of the scoring of the clusters defined by Tang et al. (n= 676 samples).

Extended Data Fig. 2. The classification of NK cells into 3 main families is robust in other blood NK cell samples.

Based on Dataset 4b. a, UMAP based on 2 independent samples of NK cells sorted from healthy human blood with clusters identified by unsupervised hierarchical clustering and their scoring with NK1, NK2, and NK3 signatures. b, c, Dot plot and UMAP visualization of some of the most discriminatory markers expressed at the transcriptional level by the three major subsets of human blood NK cells.

Extended Data Fig. 3. The NK cells in human blood can be divided into six subgroups.

Based on Dataset 1,2,3 and 4a a, Clustree plot of sc3 stability of clusters at different clustering resolution (from k=0.5 to k=1.4). b, Mean sc3 stability as a function of the granularity of clustering resolution. c, UMAP visualization of the plot of NK cells sorted from healthy human blood with clusters identified by unsupervised hierarchical clustering at a granularity of 0.7 (optimal resolution according to sc3 stability). d, UMAP visualization of the subpopulations of NK cells from the blood of healthy individuals with clusters identified by unsupervised hierarchical clustering after removing proliferating cells and populations representing less than 3% of total NK cells. e, UMAP visualization of NKG2C protein expression. f, Pie chart showing the proportion of each subgroup in the NK cell population in blood.

Extended Data Fig. 4. The NK cells in human blood can be divided into six subgroups.

a,c and d-f: Based on Based on Dataset 1,2,3 and 4a b: Based on Dataset 5. a, Bar graph showing the proportion of cells within each cluster in the datasets. (n= 13 samples) b, Violin plot of the scoring of all CD45^pos cells from Dataset 5 with the 13 genes characteristic of human NK cells as defined by Crinier et al. The error bars present the median +/- standard deviation. (n= 8 samples) c, Violin plot of the scoring with the 13 genes characteristic of human NK cells as defined by Crinier et al. The error bars present the median +/- standard deviation. (n= 13 samples) d, e, f, UMAP visualization of the expression of some key markers of NK1, NK2 and NK3 populations.

Extended Data Fig. 5. Markers of interest, functions and metabolism characterizing NK cell populations.

Based on Dataset 1,2,3 and 4a. a, Heatmap showing the differential expression of the genes composing three metabolic pathways of interest among NK cell subsets. The color scale is based on z-score-scaled gene expression. The z-score distribution ranges from −2 (blue) to 2 (red).

Extended Data Fig. 6. Dissection of the trajectory leading from NKint to NK1C.

Based on Dataset 1,2,3 and 4a. a, Dynamic heatmap of the evolution of the top 150 markers that evolve most along the pseudotime of the trajectory leading from the NKint subset to the NK1C subset.

Extended Data Fig. 7. Master regulators genes characteristic for each subset of NK cells in the blood and putative ontogeny of the main NK populations.

Based on Dataset 1,2,3 and 4a. a, Heatmap showing the differential expression of true transcription factors detected in NK cell subsets. The color scale is based on the z-score of the regulon activity. The z-score distribution ranges from −2 (blue) to 2 (red).

Extended Data Fig. 8. Distribution of NK1, NK2 and NK3 cell subsets in tissues.

Based on Dataset 7. a, ViolinPlot visualization of the module score of individual cells scored with signatures of NK1, NK2 and NK3 of ILC populations present in tonsil, lung and IELs and grouped by clusters (as defined in Dataset 7). The error bars present the median +/- standard deviation. (Lung: n = 4 samples, Tonsil: n = 6 samples, IEL: n = 4 samples).

Extended Data Fig. 9. Distribution of NK1, NK2 and NK3 cell subsets in the blood of cancer patients and at the tumor bed.

a: Based on dataset 1,2,3,4a. b-e: Based on dataset 6. a, Heatmap depicting accuracy of the label transfer for subset annotation tested on 20 % of the cells heldout to train the reference. b, ViolinPlot visualization of the maximum prediction score per cell in the blood NK samples of Dataset 6. Cells are grouped by their predicted identity. The error bars present the median +/- standard deviation. c, ViolinPlot visualization of the maximum prediction score per cell in the tumor-infiltrated NK samples of Dataset 6. Cells are grouped by their predicted identity. The error bars present the median +/- standard deviation. d, ViolinPlot visualization of the module score of individual blood NK cells of Dataset 6 scored with signatures of NK1, NK2 and NK3. Cells are grouped by their predicted identity. The error bars present the median +/- standard deviation. e, ViolinPlot visualization of the module score of individual tumor-infiltrated NK cells of Dataset 6 scored with signatures of NK1, NK2 and NK3. Cells are grouped by their predicted identity. The error bars present the median +/- standard deviation.

Extended Data Fig. 10. Distinct transcriptionnal phenotypes of NK1, NK2 and NK3 cell subsets in the blood of cancer patients and at the tumor bed.

Based on dataset 6. a, Heatmap showing the Spearman correlation between NK1, NK2 and NK3 populations in healthy individuals and across 22 different cancer types in both blood and tumor. The error bars present the median +/- standard deviation. (n= 676 samples).