Discovery of specialized NK cell populations infiltrating human melanoma metastases

NK cells represent an important arm of the innate immune system, and their name reflects their ability to recognize and kill tumor cells following activation of germline-encoded receptors (1). NK cells hold promise for cancer immunotherapy, and a variety of approaches are being developed to enhance NK cell–mediated tumor immunity. These include antibodies that block inhibitory receptors (e.g. NKG2A, TIGIT, and CD96), antibodies that inhibit MICA/B proteolytic shedding by tumor cells, recombinant IL-15 superagonist, and adoptive transfer of alloreactive or genetically engineered NK cells (2–7). In vivo models in immunocompetent mice bearing syngeneic tumors have significantly contributed to the development of such immunotherapies. However, human NK cells are substantially different from murine NK cells, including major differences in NK cell receptors and their ligands (8). In addition, most of the current knowledge of human NK cells is primarily based on analysis of blood samples, and it remains unclear how closely tumor-infiltrating NK cells resemble these circulating populations.

Investigation of NK cells has traditionally emphasized their cytotoxic function, in particular in the context of tumor immunity (9). NK cells have been viewed as innate immune cells that act at an earlier stage than adaptive lymphocytes (T cells and B cells) (1). However, two recent studies demonstrated a key role for NK cells in T cell–mediated immunity. Specifically, NK cells were found to recruit cross-presenting XCR1⁺ DCs to tumors that are critical for T cell–mediated immunity through secretion of the chemokines XCL1, XCL2, and CCL5. Infiltration of human tumors by NK cells and cross-presenting DCs was found to correlate with response to checkpoint blockade with a PD-1 antibody (10, 11). Thus, NK cells may play a broader role in tumor immunity beyond their direct cytotoxic function.

An in-depth characterization of the tumor-infiltrating NK cells may help to guide the development of NK cell–based immunotherapies for human cancers. Therefore, we performed a single-cell RNA-seq (scRNA-seq) analysis of NK cells infiltrating human melanoma metastases, including metastases from patients whose tumors had progressed following checkpoint blockade. This analysis identified potentially novel tumor-infiltrating NK cell subpopulations and highlighted striking differences between tumor-infiltrating and circulating NK cells. These data may open new avenues of investigation on NK cell immunity in human cancer.

Single-cell characterization of NK cells infiltrating human melanoma metastases. We performed a scRNA-seq analysis of NK cells from human melanoma metastases that required surgical resection. In 4 of 5 patients (CY129, CY155, CY158, and CY164), tumors had progressed following treatment with checkpoint blockade (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.133103DS1). Two of these patients had an initial response to checkpoint blockade (adaptive resistance, CY158 and CY164), while two patients failed to respond (intrinsic resistance, CY129 and CY155). NK cells were sorted by flow cytometry from tumors and matching blood samples as CD45⁺CD56⁺CD3^–CD14^–CD15^–CD163^– viable lymphocyte-size cells (Figure 1A). NK cell frequencies were lower in tumors compared with matching blood samples, as previously reported in lung cancer (12). The frequencies of NK cells in blood and tumors were also lower than the frequencies of T cells (CD45⁺CD3⁺CD14^–CD15^–CD163^– (Figure 1B and Supplemental Figure 1, A and B). The NK cell frequency among CD45⁺ lymphocytes (after exclusion of CD14⁺, CD15⁺, and CD163⁺ myeloid cells) ranged from 2.47% to 46.10% in the blood and 0.46% to 6.48% in tumors (Figure 1B). For example, in patient CY158 NK cells represented 4.04% of lymphocytes in the tumor and 30.4% of lymphocytes in the blood (Figure 1, A and B).

Figure 1

Isolation and scRNA-seq characterization of NK cells from human melanoma metastases and matching blood samples. (A) Gating strategy for the isolation of NK cells from blood and melanoma metastasis from patient CY158 by flow cytometry. (B) Percentage of blood (left) and tumor-infiltrating lymphocytes (right) from multiple patients that were NK cells (after exclusion of myeloid cells). (C) Comparison of NK cell clusters in blood versus melanoma metastasis (patient CY158) by scRNA-seq. UMAP plots were used to visualize blood and tumor-infiltrating NK cell clusters, and the percentage of NK cells in each cluster is indicated for blood and tumor NK cells (left). NK cell clusters are color coded, and key differentially expressed genes are shown for each cluster (right).

We examined the transcriptome of NK cells isolated from melanoma metastases and matching blood samples from 5 patients (CY129.2, CY155, CY158, CY160, and CY164) by scRNA-seq using the 10X Genomics platform. A nonlinear dimensionality reduction technique, uniform manifold approximation and projection (UMAP) was used to visualize NK cell clusters (13). In each patient, we observed major differences in NK cell clusters in tumors versus blood samples (Figure 1C and Supplemental Figure 2, A–C). NK cell clusters with similar gene expression profiles were identified in tumors from all 5 patients; the only exception was a NK cell population characterized by differential expression of interferon-inducible genes (including ISG15 and IFI6) that only was present in the tumor from patient CY158 (Figure 1C).

Integrated analysis of circulating NK cells from all 5 patients led to the identification of 4 blood NK cell clusters (bNK.0–bNK.3) (Figure 2A and Supplemental Figure 3A). A similar analysis of tumor-infiltrating NK cells highlighted 8 clusters (tNK.0–tNK.7) (Figure 2A, Supplemental Figure 3B, and Supplemental Excel File). Two NK cell clusters identified in tumors shared transcriptional features with blood NK cells but were present at strikingly different frequencies at these sites. The predominant NK cell cluster in tumors (tNK.0, 38.9%) shared part of the transcriptional signature with a minor blood NK cell population (bNK.2, 6.1%) (SELL, XCL1, and XCL2). In addition, the predominant blood NK cell cluster (bNK.0) shared expression of a number of key genes with tumor cluster tNK.2 (including FGFBP2, FCGR3A, PRF1, GZMB) but was smaller in tumor-infiltrating (13.5%) compared with blood (80.4%) NK cells. A small cluster of proliferating NK cells was also present in both blood (bNK.3, 1.5%) and tumors (tNK.6, 3.6%), and these cells shared expression of cell cycle genes (including PCNA and MKI67). The other 5 NK cell clusters identified in melanoma metastases (tNK.1, tNK.3, tNK.4, tNK.5, and tNK.7) were quite distinct in their transcriptional programs from blood NK cells (Figure 2). We also attempted to perform pathway analysis using the genes that were differentially expressed by these NK cell populations; however, the number of genes was not sufficiently large and did not generate enrichment scores with statistical significance for most of the gene sets. The single-cell data demonstrate that NK cells are indeed present in human melanoma metastases, including patients who progressed following treatment with PD-1 or CTLA-4 mAbs. However, a larger patient panel would be required to analyze the correlation of these NK cell populations with resistance. Tumor and blood NK cell populations from the same patients demonstrated striking transcriptional differences.

Figure 2

Single-cell RNA-seq analysis of NK cells using integrated data from 5 patients. (A) NK cell populations in blood (left) and metastases (right) analyzed by scRNA-seq using merged data from 5 patients (CD129.2, CY155, CY158, CY160, and CY164). NK cell clusters in blood and metastases were visualized using UMAP plots. The percentage of NK cells in each cluster and key differentially expressed genes for each cluster are indicated. (B) mRNA transcripts for selected genes in blood (top) and tumor-infiltrating (bottom) NK cell populations were visualized using UMAP plots. The intensity of the blue color indicates the level of expression for indicated genes in individual cells and is scaled separately between blood and tumor-infiltrating NK cells within the integrated data set from 5 patients.

NK cells and ILC3-like cells in human melanoma metastases. The markers used for isolation of NK cells (Figure 1A) are also expressed by innate lymphoid cells (ILCs), which include the ILC1, ILC2, and ILC3 subpopulations (14). ILCs are tissue-resident cells, and it was therefore likely that ILCs were not present in blood samples (15). However, it was important to assess whether some of the NK cell clusters identified in melanoma metastases could represent ILCs. A previous scRNA-seq study identified genes differentially expressed by NK cells, ILC1, ILC2, and ILC3, and we used these gene sets to assemble transcriptional signatures for each of these innate lymphocyte populations (16). As expected, all NK cell clusters in the blood had a strong NK cell gene expression signature but a low score for all 3 ILC signatures (Figure 3, A and B, and Supplemental Figure 4, A and C). Furthermore, 7 of 8 NK cell clusters identified in tumors also had a high score for the NK cell gene expression signature (Figure 3A). The exception was cluster tNK.4, which had a moderately high score for the ILC3 signature but a low score for the NK cell, ILC1, and ILC2 signatures (Figure 3, A and B, and Supplemental Figure 4, B and D). Therefore, NK cells (7 of 8 clusters, total of 92.3% cells) and ILC3-like cells (1 cluster, 7.7% of cells) were identified in human melanoma metastases.

Figure 3

Identification of NK cells and ILC3-like cells in melanoma metastases. (A and B) Identification of NK cell and ILC3 cell populations. Published single-cell data from innate lymphocytes isolated from human tonsil tissue were used to define gene expression signatures for NK cells and ILC3 (16). UMAP plots show the degree of similarity between NK cell (A) and ILC3 (B) gene expression signatures for blood and tumor NK cells. (C) Cytotoxicity gene expression signature (GZMA, GZMB, GZMH, GZMK, GZMM, PRF1, GNLY, and NKG7) for NK cells isolated from blood (left) and melanoma metastases (right). UMAP plots indicate the scores for this signature across NK cell clusters and are scaled separately between blood and tumor-infiltrating NK cells for the integrated data set from 5 patients.

Gene expression programs related to antitumor function of NK cells. The investigation of NK cells in tumor immunity has primarily focused on their cytotoxic function, but recent studies have highlighted an important role of NK cells in recruitment of cross-presenting DCs that are critical for induction of T cell–mediated tumor immunity (10, 11). We used a panel of genes encoding key cytotoxicity proteins to assemble a cytotoxicity gene expression signature (GZMA, GZMB, GZMH, GZMK, GZMM, PRF1, GNLY, and NKG7). This cytotoxicity signature was high in most blood NK cells, in particular clusters bNK.0, bNK.1, and bNK.3. In tumor NK cells, we observed a gradient for this cytotoxicity signature that was the highest in the tNK.2 and tNK.5 clusters and at an intermediate level in 5 other clusters. Only the tNK.4 cluster was negative for this cytotoxicity signature, consistent with its designation as a cluster of ILC3-like cells (Figure 3C). We also observed different patterns of expression for granzyme genes across NK cell clusters: GZMA and GZMB expression was high in most blood NK cells and showed a gradient similar to the cytotoxicity signature in tumor NK cells (Figure 2B and see below). Interestingly, GZMK expression showed a distinct pattern: it was low in most blood NK cells but high in many tumor NK cell clusters (see below).

We also observed striking differences in the expression of chemokine genes between tumor and blood NK cells. The chemokines XCL1 and XCL2 (that bind to the XCR1 chemokine receptor) were recently shown to play a critical role in recruiting cross-presenting DCs to tumors (11). Expression of these 2 chemokine genes was substantially higher in tumor NK cells (clusters tNK.0, tNK.3, tNK.6, tNK.7) compared with blood NK cells (Figure 4, A and C). In addition, we observed high expression of another set of chemokine genes (CCL3, CCL4, CCL4L2, and CCL5) in a subset of tumor-infiltrating NK cells (clusters tNK.1, tNK.2, tNK.5, and tNK.6) (Figure 4, B and C). These chemokines bind to CCR5 and other chemokine receptors and play a critical role in recruitment of T cells and other immune cells. CCL5 is also known to contribute to the recruitment of cross-presenting DCs (11, 17). Thus, tumor-resident NK cells express many chemokine genes that are important for recruitment of DCs, T cells, and other immune cell populations. It is unclear how the expression of these chemokine genes is regulated in tumors, but we note that tumor-infiltrating compared with blood NK cells expressed substantially higher levels of FOS and JUN, which encode the subunits of the AP-1 transcription factor (Figure 2B).

Figure 4

Expression of chemokine genes by blood and tumor-infiltrating NK cells. (A and B) UMAP plots showing expression of chemokine genes, XCL1 and XCL2 (A) as well as CCL3, CCL4, CCL4L2, and CCL5 (B), in blood and tumor-infiltrating NK cells. (C) Expression of each one of the chemokine genes by NK cells isolated from blood (top) and melanoma metastases (bottom). The intensity of the blue color indicates the level of expression for indicated genes in individual cells and is scaled separately between blood and tumor-infiltrating NK cells for the integrated data set from 5 patients.

The single-cell data also demonstrated functional specialization among tumor-infiltrating NK cell populations: 4 clusters of tumor NK cells showed high expression of XCL1/XCL2, while a distinct set of tumor NK cell clusters showed high expression of CCL3, CCL4, CCL4L2, and CCL5. These tumor NK cell populations may thus create distinct microenvironments. This pattern of chemokine gene expression was related to the cytotoxicity signature: NK cells with a lower cytotoxicity signature (clusters tNK.0, tNK.3, tNK.7) expressed higher levels of XCL1 and XCL2 than clusters with a higher cytotoxicity signature (tNK.1, tNK.2, and tNK.5). In contrast, CCL3, CCL4, CCL4L2, and CCL5 were expressed at a higher level by tumor-infiltrating NK cells with a higher cytotoxicity signature (Figure 3C and Figure 4, B and C). These data and recent publications (10, 11) demonstrate that the role of NK cells in tumor immunity needs to be reconsidered in a broader context: NK cells not only kill tumor cells but also recruit key immune cell populations required for protective tumor immunity.

Expression of activating and inhibitory receptors by tumor-infiltrating NK cells. NK cells integrate signals from the extracellular environment through a series of activating and inhibitory receptors (8). Among the genes encoding activating receptors, a high level of expression was observed for KLRF1 (NKp80 protein) in a large fraction of blood and tumor NK cells (Figure 5A). The AICL gene, which encodes the ligand for NKp80, is expressed in both hematological malignancies and solid tumors (18). Signals for other well-established activating NK cell receptors were lower (NCR1, NCR3, CD226, and KLRK1), but is important to note that some mRNAs tend to yield rather weak signals by scRNA-seq, even though both mRNA and protein are quite abundant in the relevant cell populations. For example, although KLRK1 mRNA (which encodes the NKG2D protein) was low in all NK cell populations, including blood NK cells, HCST mRNA was high (encodes DAP10, the adaptor molecule for NKG2D). Consistent with that explanation, published reports demonstrated that NKG2D protein can be detected on blood NK cells from melanoma patients, although at lower levels compared with healthy donors (19, 20).

Figure 5

Expression of genes encoding activating and inhibitory surface receptors on NK cells. (A and B) Expression of activating (A) and inhibitory (B) receptors in blood (top) and tumor (bottom) specimens. The intensity of the blue color indicates the level of expression of selected genes in individual cells and is scaled separately between blood and tumor-infiltrating NK cells within the integrated data set from 5 patients.

We also observed interesting expression patterns for receptors with established inhibitory function in NK cells. Tumor-infiltrating NK cells expressed higher levels of the KLRC1 gene (encodes NKG2A protein) than blood NK cells, and the KLRD1 gene (CD94 protein) was highly expressed by most tumor and blood NK cells (Figure 5B). This suggests that a large fraction of melanoma-infiltrating NK cells express the inhibitory NKG2A-CD94 receptor, which recognizes HLA-E. We also observed a strong signal for the KLRB1 gene (CD161 protein) in both tumor and blood NK cells (Figure 5B). CD161 is known to inhibit NK cell–mediated cytotoxicity following binding to the CLEC2D ligand on tumor cells and APCs (21, 22). The signals for most other inhibitory receptors were weaker, but distinct expression patterns emerged: CD96 was expressed across tumor NK cell clusters, while expression of other receptors was limited to one or a small subset of tumor NK cell clusters (such as CD160 and KIR2DL3) (Figure 5B and Supplemental Figure 5). These data identify activating (e.g., KLRF1) and inhibitory (e.g., KLRB1, KLRC1) receptor genes with high levels of expression in tumor-infiltrating NK cells.

Validation of scRNA-seq data by flow cytometry. We used flow cytometry to validate key findings from the scRNA-seq data and extend the analysis to a larger population of melanoma patients. We identified NK cells using well-established markers (CD45⁺CD56⁺CD3^–CD4^–CD8a^–CD14^–CD15^–CD163^–CD19^– viable lymphocyte-size cells), and then used the CD16a (FCGR3A gene) and FGFPB2 markers based on the scRNA-seq data to identify key NK cell subpopulations. This analysis identified 3 cell populations: (a) FGFBP2⁺CD16a⁺ NK cells that corresponded to blood (bNK.0) and tumor (tNK.2) clusters, which expressed key cytotoxicity genes (PRF1 and GZMB); (b) FGFBP2^–CD16a⁺ NK cells that could be resolved from the CD16a⁺ population using FGFPB2 as a marker; these cells corresponded to bNK.1 and potentially tNK.1, tNK.5, and tNK.6, which expressed FCGR3A but not FGFBP2; and (c) FGFBP2^–CD16a^– NK cells that corresponded to bNK.2 and abundant tumor NK cell populations, which did not express FGFBP2 or FCGR3A (primarily but not exclusively tNK.0 and tNK.3). ILC3-like cells (tNK.4) also did not express FGFBP2 or FCGR3A (Figure 2A and Figure 6, A and B).

Figure 6

Validation of NK cell subpopulations by flow cytometry. (A) Expression of genes used for identification of NK cell populations (FGFBP2 and FCGR3A) and effector molecules (GZMA and GZMK). The intensity of the blue color indicates the level of expression for selected genes in individual cells and is scaled separately between blood and tumor-infiltrating NK cells within the integrated data set from 5 patients. (B) Validation of 3 NK cell populations identified by scRNA-seq in blood samples by flow cytometry using FGFBP2 and CD16a as markers. NK cells were identified by gating on CD45⁺ and CD56⁺ cells that were negative for CD3ε, CD19, CD14, CD15, CD163, and a dead cell marker. A representative analysis is shown for patient CY165. (C) Quantification of 3 NK cell populations in blood and tumor samples based on FGFBP2 and CD16a markers. Labeling for granzymes A and K is also shown for each of the 3 populations. The genes named FCGR3A, FGFBP2, GZMA, and GZMK encode the proteins CD16a, FGFBP2, granzyme A, and granzyme K, respectively. Each dot represents an individual patient. MFI, mean fluorescence intensity. Statistical analysis was performed by 2-way ANOVA, Bonferroni’s post hoc tests; *P < 0.05, **P < 0.01, ***P < 0.001.

We then used the FGFBP2 and CD16a markers to examine these 3 NK cell populations in tumor and matching blood samples from a total of 7 melanoma patients. The predominant NK cell population in blood samples was positive for both FGFBP2 and CD16a, while a large fraction of tumor-infiltrating NK cells was negative for both markers (Figure 6C), consistent with the scRNA-seq data (Figure 6A). Therefore, CD16a and FGFBP2 markers identify 3 NK cell populations that are present with different frequencies in blood versus tumors. We further examined the expression of granzymes A and K in these NK cell populations. The expression of granzymes A and K was higher in blood compared with tumor NK cells. In addition, granzyme A levels were higher in CD16a⁺ blood NK cells compared with FGFBP2^–CD16a^– blood NK cells (Figure 6C), consistent with the scRNA-seq data (Figure 6C). Of note, flow cytometry analysis showed that granzyme K protein levels were lower in tumor-infiltrating compared with blood NK cells, while the scRNA-seq showed that GZMK mRNA was higher in tumor-infiltrating NK cells. This difference between granzyme K protein and GZMK mRNA levels could be explained by release of granzyme proteins by tumor-infiltrating NK cells following activation by tumor cells.

We also used flow cytometry to characterize melanoma cells from these metastases. HLA-A/B/C proteins were detected on the surface of melanoma cells in 8 patients but were low or undetectable in 3 other patients (Supplemental Figure 6, A and C). In particular, 1 of the tumors studied by scRNA-seq (CY155) had undetectable surface HLA-A/B/C proteins, and this patient had progressed following treatment with a PD-1 mAb (intrinsic resistance). MICA/B proteins are ligands for the activating NKG2D receptor (5). MICA/B protein was detected on the surface of melanoma cells in the majority of cases, although at a relatively low level (Supplemental Figure 6, B and D). Shed MICA was detected in serum samples from 7 of 9 melanoma cases (Supplemental Figure 6E), consistent with previous studies of proteolytic shedding of MICA/B proteins (5, 23). As expected, shed MICA was not detected in sera from healthy subjects.

Analysis of NK cells from distinct regions of a melanoma metastasis. One of the metastases (CY129.2 sample) was of sufficient size to enable analysis of NK cells from 3 distinct regions, the cortex (margin) and center (core) of the tumor as well as a distinct nodule in the cortex (Figure 7A). NK cells could be isolated from all 3 regions (Figure 7B), and 6 NK cell populations were identified in the combined data from the 3 regions (Figure 7C), which had gene expression programs similar to those described above (Figure 1C and Figure 2A). All 6 NK cell populations were identified in the 3 distinct regions of this metastasis, with some regional variation. NK cells expressing AREG, XLC1, XCL2, and FOS appeared to be enriched in the nodule (38.7%) compared with the cortex (21.6%) or center (14.5%) (Figure 7D). This analysis demonstrates that the identified NK cell populations are present in different regions of a metastasis, although systematic analysis of a large number of metastases may identify preferential localization of NK cell populations in particular tumor regions.

Figure 7

Regional analysis of NK cells in a melanoma metastasis. (A) Illustration of metastasis regions from which NK cells were isolated (patient CY129.2). (B) Flow cytometry analysis of NK cells isolated from each of the tumor regions. (C) UMAP representation of single-cell data for identification of NK cell clusters. NK cells from all 3 regions of the metastasis were combined for this UMAP plot. (D) Percentage of NK cell populations assigned to clusters (defined in C) in the 3 regions of this metastasis.

These studies identified striking differences in the transcriptome of NK cells isolated from human melanoma metastases compared with circulating NK cells from the same patients. While the majority of circulating NK cells (~80%) were represented by a cluster of cells with a well-defined cytotoxicity gene expression program, tumor-infiltrating NK cells exhibited greater diversity in their transcriptional states. We identified 7 clusters of NK cells in human melanoma metastases plus an additional cluster with an ILC3-like gene expression signature. These 7 tumor-infiltrating NK cell clusters appear to reflect functional specialization, both in terms of expression of cytotoxicity as well as chemokine/cytokine genes. In particular, we detected higher expression levels of the XCL1 and XCL2 chemokine genes in 3 clusters with lower expression of cytotoxicity genes, while expression of the CCL3, CCL4, CCL4L2, and CCL5 chemokine genes was higher in 3 other clusters with a higher cytotoxicity signature. NK cells were recently shown to recruit cross-presenting DCs to tumors that are critical for CD8 T cell–mediated tumor immunity. These DCs express the XCR1 chemokine receptor that is activated by the XCL1 and XCL2 chemokines (11, 24). The distinct clusters of tumor-resident NK cells may therefore play different functional roles in tumor immunity, including direct cytotoxicity against malignant cells as well as recruitment of other important immune cell populations, including cross-presenting DCs.

A recent scRNA-seq study compared circulating human NK cells with tissue-resident NK cells in the spleen; the findings of this study identified 2 blood and 4 splenic NK cell clusters. Similar to our observations, blood and splenic NK cells showed large differences in their transcriptomes (25). NK cells from the major splenic cluster differentially expressed the FGFBP2, FCGR3A, PRF1, and GZMB genes, similar to the predominant NK cell cluster in the blood (bNK.0 in our study) and 1 of the melanoma NK cell clusters (tNK.2). NK cells from another splenic cluster differentially expressed the XCL1 and XCL2 chemokine genes, similar to several of the melanoma NK cell clusters (including tNK.0). This study and our data are consistent with the hypothesis that the transcriptional state of NK cells is dynamically regulated by the tissue microenvironment.

Our data demonstrate that NK cells are present in human melanoma metastases, even in patients who failed to respond to checkpoint blockade. This study was not designed to address which of the identified NK cell clusters may be related to responsiveness/resistance to checkpoint blockade. This question could be addressed in a future prospective study that includes pretreatment and on-treatment tumor biopsies at defined time periods. It will also be of considerable interest to determine if the different NK cell populations are localized in distinct microenvironments within tumors. The differentially expressed genes that characterize each NK cell cluster represent markers for the identification of these NK cell populations in tissue sections by multicolor immunofluorescence.

These single-cell data provide a wealth of information on human NK cells from melanoma metastases and may inform the future development of strategies for enhancing NK cell–mediated immunity against tumors. For example, the identification of activating and inhibitory receptors expressed by NK cell subpopulations within tumors may enable therapeutic manipulation of these cell populations. The apparent functional specialization of NK cell populations within tumors may also provide opportunities to increase distinct aspects of NK cell function, such as cytotoxic activity or secretion of chemokines that recruit cross-presenting DCs.

Characterization of human NK cells infiltrating melanoma metastases. Freshly resected tumor tissue was dissociated using a Tumor Dissociation Kit (Miltenyi Biotec, catalog 130-095-929). Red blood cells in PBMC samples were lysed using ACK buffer. NK cells were isolated using a BD Aria flow cytometer (BD Biosciences) from melanoma lesions and corresponding blood samples as lymphocyte-size single viable cells that were CD45⁺ (Biolegend, catalog 304037 and 304022), CD56⁺ (Biolegend, catalog 355504 and 318312), CD3ε^– (Biolegend, catalog 300447, 300330, 317322), CD4^– (eBiosciences, catalog 4345826), CD8a^– (Biolegend, catalog 301038), CD14^– (Biolegend, catalog 301843), CD15^– (Biolegend, catalog 323028), and CD163^– (Biolegend, catalog 333628). These sorted NK cells were used for scRNA-seq analysis.

In follow-up experiments, tumor and blood NK cells were also stained with anti-FGFPB2 (Biolegend, catalog 346603), GZMA (Biolegend, catalog 507221), GZMK (Biolegend, catalog 370513), CD62L (Biolegend, catalog 304808), NKG2D (Biolegend, catalog 320828), and CD16a (Biolegend, catalog 302038) antibodies and analyzed by flow cytometry. All fluorochrome-labeled antibodies were purchased from BioLegend, BD Biosciences, or eBiosciences. For intracellular staining, NK cells were fixed and permeabilized with the True-Nuclear Transcription Factor Buffer Set (Biolegend, catalog 424401) according to the recommendations of the manufacturer. Samples were analyzed using a Cytoflex flow cytometer (Beckman Coulter).

scRNA-seq. Immediately after sorting of NK cells (~13,000 NK cells per sample), cell suspensions were washed in 0.05% RNase-free BSA in PBS. The 10X Genomics 3′ V2 single-cell assay (10X Genomics) was used for construction of scRNA-seq libraries. Reverse transcription, cDNA amplification, and library preparation were all performed according to the manufacturer’s instructions. Libraries were sequenced using an Illumina HiSeq 2500 on rapid-run mode or an Illumina HiSeq 4000, which yielded >25,000 reads per cell.

Computational analysis of scRNA-seq data. We used the 10X Genomics’ Cell Ranger software for the demultiplexing, alignment, filtering, barcode counting, and unique molecular identifier (UMI) counting steps. The analysis was performed using the Seurat 3.0 package (26). We first processed each individual data set separately prior to combining data from multiple samples. For each data set, we selected the 1,500 most variable genes. Subsequently, we ran principal component analysis (PCA) and used the first 15 principal components (PCs) to perform Louvain clustering and UMAP embedding (13, 27). We checked well-defined marker genes for each cluster to identify potential contaminating cell populations, such as T cells (CD3D, CD3E, and CD3G), B cells (IGHG1, IGHG2, and JCHAIN), macrophages (LYZ), and melanoma cells (MLANA); these cells were removed prior to subsequent analyses. We also removed outlier cells with very high total UMI counts (top 0.2 %) from each data set.

We then compared the paired blood and melanoma-infiltrating NK cell populations of each patient using Seurat’s integration algorithm (26) and then also separately integrated the 5 melanoma-infiltrating NK cell samples and 5 blood NK cell samples. We used the 3,000 most variable genes from each sample and the first 15 PCs to choose 1,000 anchor genes for integration. Afterward, we repeated PCA, clustering (resolution = 0.3 for melanoma-filtrating NK cells and 0.16 for blood NK cells), and UMAP embedding on the integrated data sets. Finally, we performed differential expression tests on the integrated data sets to identify the genes significantly upregulated in each cluster compared with all other cells (adjusted P < 0.05) as well as the genes differentially expressed between blood and melanoma-infiltrating NK cells within each major cluster.

For gene sets representing specific cellular functions or pathways, we also computed AUCell scores for each cell in order to evaluate the variation of activity of gene sets across cell population (28). We handpicked gene sets for cytotoxicity and chemokine activity based on analysis of the data. For NK cell, ILC1, ICL2, and ILC3 identities, we referred to the gene signatures defined by a previous study using the genes significantly upregulated in 1 of these innate cell types compared with the other 3 populations (adjusted P < 1 × 10^–3); specifically for ILC identities, the same genes were not significantly upregulated in NK cells (adjusted P > × 10^–2) (16).

We also performed a ranked gene set enrichment analysis (GSEA) with the GSEA 4.0.1 software (Broad Institute Inc., Massachusetts Institute of Technology, and Regents of the University of California) using the differentially expressed genes (Supplemental Excel File).

Data and materials availability. Processed data generated for this study are available through the Gene Expression Omnibus database (GSE139249).

Statistics. All statistical analyses were performed using GraphPad Prism 8 software. We used 2-way ANOVA and Bonferroni’s post hoc tests for statistical analyses.

Study approval. Melanoma tissue samples were obtained from patients who required surgery at Brigham and Women’s Hospital for the treatment of nonresponsive lesions. Blood samples were also collected at the time of surgery. The human subjects protocol was approved by the institutional review board of Brigham and Women’s Hospital, and written informed consent was received from all participants prior to inclusion in the study.