A Gene Signature Predicting Natural Killer Cell Infiltration and Improved Survival in Melanoma Patients

Immunotherapies have improved clinical treatment for several cancer types, including renal cell carcinoma (1), non–small cell lung cancer (2), hematologic malignancies (3), and melanoma (4). Although these treatments show promise, their efficacy is restricted to a subset of tumor types and patients with “immune hot” tumors. Improved understanding of how the innate immune system contributes to antitumor responses can help to open new therapeutic avenues, to increase the potential and efficacy of clinical treatments, and to expand the range of targeted malignancies. Natural killer (NK) cells, a subset of innate lymphoid cells, are effectors of innate immunity. The cytotoxic capabilities of NK cells allow them to kill tumor cells even at a relatively low ratio (e.g., 1:1; ref. 5). NK cells are necessary for clearance of cells that carry a viral burden or have undergone oncogenic transformation, and several in vivo studies have demonstrated a role for NK cells in limiting the metastatic dissemination of melanoma (6–9). NK cells can be targeted by immunotherapeutics (10, 11). Furthermore, NK cells can initiate an antitumor T-cell response by recruiting conventional type-1 dendritic cells (cDC1) through chemokine signaling (via XCL-1 and CCL5/RANTES; ref. 12), and support stimulatory DCs (sDC) through the expression of the Flt3 ligand (Flt3l in mice and FLT3LG in humans; ref. 13).

Regulators of NK cell activity include the cytokine IL15 (14), chemokines such as CCL5 (RANTES; ref. 15), growth factors such as TGFβ (16, 17), and cytokine-inducible SH2 containing protein (CIS), which modulates intracellular JAK–STAT signaling (18). Evidence suggests that modulation of NK cell populations is feasible for cancer treatment (19). Treatments based upon systemic administration of IL15 constructs have shown promise in leukemic and solid tumors (20–22). Although IL15 alone can stimulate immune targeting of cancers, cytokine-mediated NK cell activation and expansion can be increased in combination with IL12 and IL18 (23) or further amplified through deletion of CIS (encoded by CISH), a negative regulator of cytokine signaling and effector function. Cish^−/− mice are resistant to a range of metastatic cancers (18). Thus, it is clear that we are yet to fully elucidate the full repertoire of molecular systems that regulate in vivo NK cell activity.

Advances in nucleic acid sequencing technology and associated methods for data analysis have allowed the application of transcriptomic profiling to complex tumor samples (24, 25). Resulting data have enabled the development of mathematical methods, such as CIBERSORT (26), that infer the relative abundance of immune cells that have infiltrated into solid tumor samples. Although these tools have improved our understanding of immune infiltration (27–30) and advanced cancer research, the complexity of model fitting procedures makes it difficult for researchers to modify or investigate these gene lists. In addition, public tumor transcriptomic data offer opportunities to identify how changes in tumor phenotype are associated with changes in the relative abundance of immune cell subpopulations.

We have developed a single-sample gene-set scoring method that uses a rank-based metric to quantify the relative enrichment of a gene set within a sample transcriptome (31). Here, we have combined NK cell signatures from the LM22/CIBERSORT (26) and LM7 (32) gene sets and curated this list to produce a gene set that reflects the relative abundance of NK cells within a tumor sample. As melanoma tumor cells are highly immunogenic, we have focused upon the analysis of TCGA skin cutaneous melanoma (SKCM) data (33). We show that the relative expression of NK cell genes within metastatic tumors is associated with a survival advantage. Our scoring approach can be used to explore putative modulators of NK cell activity by examining their association with NK score and survival effects associated with their expression.

Data used in this study are available from listed repositories (Table 1). For TCGA SKCM data, RSEM abundance data without normalization were downloaded from the genomic data commons. For sorted immune cell populations (GSE60424, ref. 34; GSE24759, ref. 35) and melanoma cell line data (E-MTAB-1496), processed transcript abundance data were downloaded and used directly. For GSE24759, only samples derived from peripheral blood were used: data from colony-forming samples were excluded to exclude culturing effects, and CD56⁻/CD16⁺/CD3⁻ mature NK cell data were excluded due to apparent batch effects. Although the relative log expression plot appeared normal for NK cell samples (Supplementary Fig. S1A), a principal component analysis (PCA) indicated that the CD56^LoCD16^Hi NK cell samples showed separation from other NK cell samples on PC1 and PC3 (Supplementary Fig. S1B) with several samples clustering toward granulocytes and monocytes on PC4 (Supplementary Fig. S1C). In cases with gene multimapping (multiple probes/probe sets per gene), median values were used. To better examine the relative expression of marker genes across immune cell subsets, we examined CD45⁺ (“nonmalignant”) single-cell RNA-seq data from dissociated melanoma samples with annotated cell types (GSE72056; ref. 36). Further, single-cell RNA-seq from (checkpoint inhibition naïve) dissociated melanoma samples (GSE120575; ref. 37) were used for subsequent comparison and visualization.

Using relative log expression (RLE) plots to explore unwanted variation (38), log-transformed transcript abundance data (downloaded directly from GEO) were median-centered for each gene, and then within each sample the difference between the observed and population median of each gene was calculated. For the PCA, genes with abundance above the 10th percentile (5.34) within at least 4 samples (corresponding to the smallest sample group) were retained. Data were normalized using sklearn StandardScaler, before calculating principal components using the sklearn PCA function.

Genes identified in Fig. 1 were annotated as immune-associated using the Gene Ontology category GO:0002376—“immune system process”—together with all descendant or child processes (39). The original SKCM manuscript (33) was used for the associated TCGA immune gene list and TCGA classifications of “Immune high” patients.

Preprocessed single-cell RNA data were downloaded from respective sources (Table 1) and were analyzed through the Seurat (v. 2.3.4) pipeline (40).

Single-cell RNA-seq data from Tirosh and colleagues (36) have annotated cell types and thus CD45⁺ “nonmalignant” cells were used for signature curation. Cells with abundance data for fewer than 1,500 genes or more than 10,000 genes were removed. To improve downstream clustering and visualization of immune subtypes used, endothelial cells were removed before performing a PCA across the 1,000 most-variable genes, taking the first 50 principal components for uniform manifold approximation and projection (UMAP; ref. 41) visualization, with a “minimum distance” of 0.2 and “number of neighbors” equal to 50.

Preprocessed CD45⁺ single-cell RNA-seq data from treatment-naïve samples were taken from Sade-Feldman and colleagues (37). Cells with abundance data for fewer than 1,000 genes or more than 5,000 detected genes were removed, as were cells with more than 5% of reads derived from mitochondrial genes. PCA was performed across the 3,550 most-variable genes, and the top 20 principal components were used for visualization with UMAP, using a “minimum distance” of 0.1 and “number of neighbors” equal to 30.

A schematic overview of the workflow for NK signature curation is given in Supplementary Fig. S2. With the aim of applying a unidirectional (expected upregulated) gene set with singscore (31), a preliminary NK cell gene set was created by combining all “expected upregulated” genes from the CIBERSORT (LM22) active and resting NK cell gene sets (26), the LM7 NK cell gene set (32), and human orthologs for established NK cell markers from a number of mouse studies (refs. 18, 42, 43; collectively referred to as the “Huntington gene list”; Supplementary Fig. S2).

To improve the specificity of this signature against other immune subsets, we examined bulk RNA-seq data from sorted cell populations and single-cell RNA-seq data from melanoma metastases. SRA files for immune cell bulk RNA-seq data from healthy individuals were downloaded in September 2016 and converted to FASTQ format using the SRA toolkit. Reads were mapped to human genome hg19 using Rsubread (v. 1.32.0; ref. 44) and counts were quantified using featureCounts. Count-per-million (CPM) and RPKM values were calculated using edgeR (v. 3.24.0). For differential expression analysis, genes were retained if their abundance exceeded a threshold (CPM > 4) within at least 4 samples (i.e., all samples for a sorted cell type). The voom-limma (v. 3.38.2) pipeline (45) was used by applying the TREAT criteria (ref. 46; log|fold change| > 1) to calculate the t statistics, log|fold change|, and adjusted P values for all genes when comparing NK cells against all other sorted cell populations (i.e., excluding whole blood). The Homo.sapiens (v. 1.3.1) package was used for annotation. Results for NK signature genes are given in Supplementary Table S1.

For CD45⁺ single-cell RNA-seq data from Tirosh and colleagues, if genes had evidence of extensive dropout (median value below 0.5 for all cell types), they were retained if the 75th percentile value for NK cells was higher than all other cell types and greater than 2.67 (90th percentile of abundance across all genes and cells, for genes with logTPM > 1 in at least 30 cells; Supplementary Fig. S2). All other genes (which passed the above dropout criteria) were retained if the median expression within NK cells was greater than the 75th percentile of all other cell types. Genes were retained if they passed these criteria within both the sorted bulk RNA-seq and single-cell RNA-seq data.

Next, we improved the specificity of our NK signature against genes that are expressed across solid tumors. For all genes, we examined the expression within adherent cell lines (i.e., hematopoietic and lymphoid cell lines were excluded) using RNA-seq data from the Cancer Cell Line Encyclopedia (CCLE; ref. 47). Genes were retained if their median expression (logTPM) was below 5.23 (25th percentile of nonzero abundance across all genes and cell lines). Results of these tests are given in Supplementary Table S1.

Gene set scoring was performed using the singscore approach (31). Briefly, genes are ranked by increasing transcript abundance, and for a set of target genes, the mean rank is calculated and normalized against theoretical minimum and maximum values. If directional gene lists are provided (i.e., a set of genes expected to be upregulated and a set of genes to be downregulated), as with the TGFβ EMT signature (48), then the mean rank of expected upregulated genes is calculated from genes ranked by increasing abundance. The mean rank of expected downregulated genes is calculated from genes ranked by decreasing abundance. These values are then normalized and summed. Accordingly, a high gene set score indicates that the pattern of gene expression in a sample is concordant with the pattern captured by the gene-expression signature.

As noted in the results and shown in Supplementary Fig. S3A and S3B, there are large survival differences between patients with primary and metastatic tumors. To avoid confounding effects from this, unless otherwise stated, we have focused on patients with metastatic tumors only who also had valid age and survival data. One patient with both a metastatic and primary sample was excluded.

Cox proportional hazard models and Kaplan–Meier survival curves were generated using python lifelines (v. 0.14.6; DOI: 10.5281/zenodo.1303381) with standard parameters. For individual gene Cox hazard models, a Bonferroni correction was applied to correct for multiple-hypothesis testing. Unless otherwise specified, separation of patient samples with a single parameter (e.g., age, gene expression) used the 33rd and 66th percentile values to threshold, and separation of patient samples with two parameters used the 40th and 60th percentile values to threshold.

General analyses and visualization were performed using python v3.6 with: pandas (49) for data handling; scipy (v 1.1.0), scikit-learn (v. 0.19.2), and numpy (v. 1.14.5) for numerical calculations; and matplotlib (v. 2.2.3) for plotting/visualization. For analyses with R, dplyr (v. 0.7.8) and tidyr (v. 0.8.2) were used for formatting data, and ggplot2 (v 3.1.0), gridExtra (v. 2.3), and RColorBrewer (v. 1.1-2) for visualization.

Computational scripts used in this work are available from our GitHub repository: https://github.com/DavisLaboratory/NK_scoring

Cutaneous melanoma is an ideal target for immunotherapy as the high mutational burden of this malignancy is associated with the generation of neoantigens, which can induce an immune response (50). Several reports demonstrate that immune infiltration signatures provide a prognostic indicator in melanoma (51), including the TCGA SKCM study, which demonstrated that this effect was independent of the underlying genomic subtype of the melanoma (33).

Due to significant survival differences between patients with a primary or metastatic tumor (Supplementary Fig. S3A and B), this report focuses on the 365 patients with only metastatic tumor samples. There are significant survival effects associated with the patient's age at diagnosis (Supplementary Fig. S3C), whereas male and female rates were not significantly different (Supplementary Fig. S3D and Table 2; Cox proportional hazards model, P = 0.20). Similarly, no tumor sites showed significantly different hazard coefficients (Supplementary Table S2).

To analyze survival effects associated with individual genes, we built a series of Cox proportional hazard models for each gene where patient age at diagnosis was included as the only covariate (together with transcript abundance for that gene). As shown in Fig. 1A, when focusing upon highly significant genes with a negative hazard coefficient, many are annotated as modulators of immune function (genes in bold) and a number of remaining genes are well-known immune modulators that lack associated GO annotations (e.g., CLEC2B, CD72, SRGN, and MIR155HG). The associated patient survival curves are shown for patient age (Fig. 1B), and a selection of genes (Fig. 1C–E). Higher expression of the hallmark inflammatory cytokine encoded by IFNG corresponds to improved survival outcomes (Fig. 1C), and several interferon-induced genes are also associated with a hazard reduction (e.g., IRF1, IFITM1; Fig. 1A). High tumor transcript abundances of the NK cell marker gene KLRD1 (Fig. 1D; also known as CD94) or the cytokine IL15 (Fig. 1E), which regulates NK cell (14, 42, 52) and T-cell activity (53), are also associated with improved long-term survival outcomes. Transcript abundance for the B2M gene encoding β2 microglobulin (Fig. 1F) has one of the most negative hazard coefficients, likely reflecting its role in MHC class I antigen presentation of neoantigens to CD8⁺ T cells and consistent with reports of this process for immune control of tumors (54). Further, a truncation mutant of B2M can confer resistance to PD-1 blockade in melanoma (55), and mutations in B2M have been shown to disrupt immune surveillance in lung cancer (56). The large negative hazard coefficient associated with HAPLN3, encoding a hyaluronan and proteoglycan link protein, suggests that this gene may warrant further investigation in the context of immune recognition and targeting.

Several transcriptomic data deconvolution methods have generated gene signatures that are predictive of tumor infiltration by specific immune cell subpopulations. For this work, we examined transcriptomic data from sorted immune cell populations and NK cell signatures from the LM22 (26) and LM7 (32) gene sets. A common critique of immune deconvolution methods is the high colinearity/cross-correlation between different signatures (32). Although this can be attributed to the similar transcriptional profiles of some lymphocyte subsets (demonstrated by similar positions of sorted NK and T-cell populations in PCA plots; Supplementary Fig. S1B and S1C), to an extent it also represents the cascading series of intercellular interactions that mediate immune activation within complex tissue samples. Accordingly, several immune-associated gene subsets are cross-correlated to a varying extent (Supplementary Fig. S4). To address this, we combined and curated the NK cell signatures from the LM7 and LM22 gene lists, together with human orthologs for well-known mouse NK marker genes (Supplementary Fig. S2; details given in Materials and Methods). Genes were retained if they had higher expression in NK cells relative to other immune populations as well as relatively low expression across adherent cell lines from the CCLE (Supplementary Table S1).

The relative expression of these 20 marker genes is shown across immune subsets from melanoma single-cell RNA-seq data (Fig. 2A), as well bulk RNA-seq (Fig. 2B) and microarray (Fig. 2C) data from sorted immune cell populations. Many NK marker genes have some expression within CD4⁺ and CD8⁺ T-cell populations (Fig. 2B and C). In particular, the cytotoxic mechanisms used by NK cells share many similarities to those used by CD8⁺ T cells, including secretion of granzymes (e.g., GZMA, GZMB, GZMK, and GZMM) and perforin (PRF1; ref. 57). As demonstrated by genes without extensive dropout in the single-cell RNA-seq data (e.g., CTSW, GZMB, IL2RB, KLRD1, and NKG7), however, the minimum abundance within NK cells tends to be greater than the 75th percentile of abundance within T-cell populations. For genes with minor dropout issues (e.g., NCR1, IL18RAP, XCL1, and XCL2; median abundance for all cell types around 0), the 75th percentile of abundance within NK cells corresponds to relatively high expression (logTPM of approximately 5 or above). There does appear to be a subset of unresolved cells (Fig. 2A, gray) with relatively high expression of genes including CTSW, GZMB, and NKG7. Examining the cell classifications by Tirosh and colleagues, many of these unresolved cells show expression of markers from both B-cell and T-cell populations, and thus these may be unresolved cell doublets.

Given our use of annotated cell subsets from Tirosh and colleagues for NK cell marker curation, we next examined the expression of selected lymphocyte markers across an independent, larger set of single-cell RNA-seq data from Sade-Feldman and colleagues. As shown (Fig. 3, top), annotated NK cells within the Tirosh data have little or no expression of T-cell markers such as CD3D and CD4, and lower expression of CD3E, whereas there is strong expression of remaining NK cell marker genes. Note that FCGR3A and NCAM1 were excluded from our curated signature (Fig. 2) due to overlap with other immune subsets and solid cancer cell lines, respectively; however they are included here as established markers to delineate lymphocyte populations. Examining the Sade-Feldman data, we find a subset of cells with lower expression of CD3E and no expression of CD4, yet relatively high expression of other indicated NK marker genes. KLRF1 (NKp80) appears to have particularly high expression within this subset that we have annotated as NK cells.

Gene set enrichment analyses are commonly used after differential expression to assess whether genes with the greatest changes are enriched for classifications of specific pathways or processes. An alternative “relative approach” (25) is to analyze the gene-expression patterns (transcript abundances) of individual samples and calculate the relative concordance of each one against specific gene set that has been defined to capture a particular molecular phenotype.

We have developed a gene-set scoring method, singscore (31), which uses the normalized mean rank of genes that are associated with a molecular phenotype or cellular behavior (48, 58). With this approach, a difference in score between two samples can be related to the percentile change in mean rank of the gene set, providing a metric that summarizes the concordance between the gene-expression profile of an individual sample and the specified gene sets. Using this scoring method with “Immune cluster” genes from the original TCGA SKCM publication (33) and introducing a single threshold, we can largely recapitulate their classification of “Immune High” samples (Supplementary Fig. S5A and S5B). Furthermore, we can extend this analysis to samples that were later added to the TCGA SCKM cohort (Supplementary Fig. S5C) and show a similar survival effect (Supplementary Fig. S5D).

We used our curated 20 gene NK cell signature to score metastatic tumors from the TCGA SKCM cohort. When samples are sorted by increasing NK score, the concordant expression pattern of these genes across metastatic melanoma samples becomes apparent: as expected, all these genes carry a positive correlation with the NK score (Fig. 4A). Using our NK signature to partition patients, there are strong survival effects for patients with either a high or moderate NK score (Fig. 4B). These survival effects are largely recapitulated across a selection of marker genes (Fig. 4B) associated with: T-cell infiltration (CD3D), cytokine signaling, which promotes NK cell and T-cell infiltration (IL15), a shared component of the IL2 and IL15 receptors (IL2RB), interferon-gamma responsive checkpoint markers (CD274), NK secreted chemokines/cytokines associated with dendritic cell recruitment and stimulation (CCL5, XCL1) and cytotoxic effector molecules (GZMB, FASLG).

Next, we investigated the relative survival effects of key cytokine and cytotoxic components linked to NK cell function (IL15, XCL1, XCL2, CCL5, FLT3LG, GZMA, GZMB, and FASLG) within patient samples from the low NK score and high NK score groups (Supplementary Fig. S6). For metastatic tumors with a low NK score, IL15 and FLT3LG showed some survival benefit; however, none of the genes examined (including these) showed a significant (P < 0.05) association with patient survival. Conversely, when examining metastatic tumors with a high NK score, there appears to be improved patient survival linked to relatively high expression of the chemokines XCL2 (also shown in Fig. 4C; P = 0.035) and CCL5 (P = 0.031). There are no survival differences in this group for the cytotoxic genes investigated (Supplementary Fig. S6), exemplified by GZMB (Fig. 4C).

It is difficult to disentangle the survival effects of other cytotoxic/effector immune cells, as improved survival effects are also seen for the T-cell marker CD3E (Fig. 4B). This is further reflected by the strong association between our NK score and scores derived from the TCGA immune cluster genes (Supplementary Fig. S7A) or the T-cell signature (Supplementary Fig. S7B), which share similar survival effects (Supplementary Fig. S7C and D). However, NK cells play a role in initiating the intercellular signaling cascade necessary for immune cell recruitment. We find good concordance between the Böttcher 5-gene NK cell signature (12) and our NK signature score (Supplementary Fig. S7E), as well as between our NK score and a score calculated using the Böttcher 4-gene DC cell signature (Supplementary Fig. S7F), again with similar survival effects (Supplementary Fig. S7G and S7H). The relative abundances of T-cell signature genes within TCGA SKCM samples and sorted immune populations are shown for reference (Supplementary Fig. S8). A number of genes associated with cytotoxic CD8⁺ T cells show relatively high expression within NK cell samples.

An advantage of the singscore approach is that it allows gene signatures associated with cell phenotype to be combined for further investigation. Work on innate anti–PD-1 resistance in melanoma found similarities with markers of MAPK inhibitor resistance (59), and melanoma phenotype switching has been linked to general drug/MAPK inhibitor resistance (60). To investigate this further, we examined the relative association between our NK score and several phenotype-switching–associated gene set scores, including: (i) a proliferative, epithelial phenotype; (ii) an invasive, mesenchymal phenotype; and (iii) a mesenchymal phenotype where EMT has been induced by TGFβ (48). We found no association between NK score and epithelial score (Fig. 5A); however, there was an association with mesenchymal score (Fig. 5B), such that tumors with less mesenchymal characteristics have lower NK scores, whereas highly mesenchymal tumors had a range of NK scores, suggesting that a subset of these tumors had higher NK cell infiltration. There was no association with NK score relative to the TGFβ-specific EMT gene score (Fig. 5C), although there was a positive association between TGFβ EMT score and mesenchymal score (Fig. 5D). As we had previously observed, whereas samples with a high TGFβ EMT score also had high mesenchymal gene-expression scores, a subset of highly mesenchymal samples showed no evidence of TGFβ-driven EMT. Neither the TGFβ EMT score (Fig. 5E) or NK score (Fig. 5F) had any association with age; however, when we partitioned patients by NK score and TGFβ EMT score, younger patients with evidence of good NK cell infiltration and low TGFβ activity had significantly favorable survival outcomes (Fig. 5G); this effect was not present for older patients (Fig. 5H).

Although most in vivo experiments indicate a primary role for NK cells in limiting metastatic colonization (7), these data suggest that not only are NK cells associated with established metastatic tumors, but the presence of NK cell infiltrate is associated with an improved prognostic outcome.

A role for NK cells in driving a robust immune response has been demonstrated through the effects of the NK–DC cell axis in modulating immunotherapy response (13). To further investigate potential modulators of NK cell infiltration, we next examined transcriptomic data from the LM-MEL panel, which contains representative cell lines for both the proliferative and invasive phenotype. Gene sets were filtered to retain only those present in both the TCGA and LM-MEL data, and gene set scoring was repeated for both data sets to facilitate comparison between tumor samples and the corresponding cell line models (Supplementary Fig. S9A–S9C). In the absence of a high NK score, patients with a high mesenchymal score show no survival effects associated with the TGFβ-EMT score (Supplementary Fig. S9D).

A number of melanoma cell lines from the LM-MEL panel appear to be associated with various subsets of high/low mesenchymal score and TGFβ-EMT score (Supplementary Fig. S9A–S9C; colored scatter markers). By contrasting genes correlated or anticorrelated with NK score across the TCGA data against these cell line data, we can identify markers that may be derived from the melanoma tumor and exert an immunomodulatory effect (Supplementary Fig. SE). To demonstrate the association with phenotype switching, the markers CDH1 and MITF are included (61, 62). Further, TGFβ activity in these cell lines is associated with THBS1 (63). Consistent with our observation that NK score tends to be higher in more mesenchymal tumor samples, many of the positively correlated genes tend to have higher expression in the MITF-low cell lines. Similarly, many of the anticorrelated genes tend to have lower expression in the MITF-low cell lines.

Several notable genes are present within these lists (Supplementary Fig. S9E). Again, B2M is identified, together with several HLA genes, suggesting that more mesenchymal-like melanomas may be more immunogenic in part because of increased antigen presentation associated with this phenotype. The presence of IL18 in this list may be partially explained by the known links between IL18 and NK cell activity (64). From the genes inversely correlated with NK score, we note that CMTM4 is a positive regulator of PD-L1 (together with its paralog CMTM6 which is not present; ref. 65) and appears to have lower expression within more mesenchymal cell lines.

Immune “checkpoint” inhibitor antibodies, which drive tumor-resident cytotoxic lymphocyte hyperactivation, have revolutionized cancer therapy. Although much clinical research is directed toward programs that underlie immunotherapy resistance and immune-related adverse events, we lack an in-depth understanding of the mechanisms dictating response, and we do not have markers to identify patients who are the most likely to respond in the context of metastatic disease. Checkpoint inhibitors primarily block inhibitory pathways in tumor-resident T cells; however, interest in other effector populations, such as NK cells, is growing (10), with studies showing that NK cells have a role in immunotherapy success (13).

Clinically, peripheral blood NK cell activity is inversely correlated with cancer incidence (66). NK cell infiltration in human tumors is associated with better prognosis in squamous cell lung carcinoma, as well as gastric and colorectal carcinomas (7). In melanoma cells, researchers have found specific HLA-I allelic losses in up to 50% of patients analyzed, and even when expressed on melanoma cells, expression of specific HLA class I molecules is often insufficient to inhibit NK cell–mediated cytotoxicity (67). These data hint that metastatic melanoma may be susceptible to NK cell–mediated killing and therapies that enhance NK cell activity should be investigated further. Along these lines, cutaneous, subcutaneous, and lymph node melanoma biopsies from a cohort of stage III and IV unresectable metastatic melanoma patients being treated with anti–PD-1 (pembrolizumab) were analyzed for immune infiltrate (13). Metastatic melanoma biopsies from patients who responded to pembrolizumab had higher NK cell infiltration compared with nonresponders, and this NK infiltration was correlated with DC infiltration (13). However, Barry and colleagues did not observe an association between regulatory and effector T-cell populations in responders (13), even though CD8⁺ T-cell proliferation has been reported as a marker of pembrolizumab on-target effect and tumor regression (68).

Although the kinetics of immune infiltration into tumors are difficult to study directly, Barry and colleagues suggested that NK cells support DC persistence and survival by producing the DC growth factor FLT3 ligand, and demonstrated these effects using an Flt3l transgenic mouse model (13). This evidence together with another report (12) suggests that NK cell infiltration may precede the majority of DC infiltration because intratumor NK cells are a major source of XCL1, a chemoattractant for XCR1⁺ DCs. These data on melanoma immune infiltrate motivated our study, in which we have investigated the transcriptomic and matched clinical data available through the TCGA SKCM cohort (33).

We found that metastatic melanoma tumors that have a high NK score are associated with better patient survival, consistent with the results from Böttcher and colleagues and a range of in vivo animal survival studies (6, 17). This effect was recapitulated by the NK cell–secreted chemokines XCL1 and CCL5, as well as NK cell effectors GZMB and FASLG. Our comparison of cytokine and cytotoxic gene expression within the high NK score patient subset suggests that greater expression of cytokines such as CCL5 and XCL2 may have a greater effect on survival than higher expression of cytotoxic effectors such as GZMB and FASLG. This may represent a saturation of the effect of cytotoxic molecules within activated cytotoxic lymphocytes (already expressing these molecules), but this may also reflect the ability of NK cells to recruit and support DCs, driving amplification of the antitumor response through intercellular signaling programs. Although the gene set derived by Böttcher and colleagues performs well (12), they calculated NK scores using mean log₂ abundance data, which can be susceptible to outliers and places a greater weighting on genes with high transcript abundance. We demonstrated that although singscore was developed for larger gene sets, it performs well with the 5-gene Böttcher NK signature (NCR3, KLRB1, PRF1, CD160, and NCR1).

It is hard to compare the accuracy of these signatures without validation data for NK cell infiltration. Our results demonstrate the application of our computational method for estimating the abundance and heterogeneity of different immune subsets across different tumors and patients. Our methods allow variations in these relative immune scores to be compared against other phenotype-associated gene sets. Indeed, we extracted survival effects when our NK score was examined in the context of melanoma phenotype switching and TGFβ signaling. Phenotype switching is a regulatory program involved in the progression of melanoma (62, 69) that has been linked to vemurafenib resistance (70) and general drug resistance (60). It allows tumor cells to transition between proliferative (“epithelial-like”) and invasive (“mesenchymal-like”) behaviors. In melanoma, TGFβ drives EMT/phenotype-switching programs (48), which are mediated by signaling molecules such as thrombospondin 1 (63).

Although long-term survival effects are associated with patient age, we found no association between NK score and age (Fig. 5E). However, younger patients with low TGFβ-EMT and high NK scores had a significant survival advantage (Fig. 5G; P < 1 × 10^–3 compared with both NK^lo with TGFβ-EMT^lo, or NK^lo with TGFβ-EMT^hi). These results suggest that younger patients may receive a greater benefit from NK targeted immunotherapies, perhaps reflecting a higher capacity for a robust immune response in young patients.

Tumor inflammation influences patient response to immune-checkpoint blockade (68). Why infiltrate is absent or heterogeneous for certain tumor types is an outstanding question that drives immunotherapy drug R&D efforts in order to increase indications for class-leading drugs such as pembrolizumab and nivolumab. A large pan-cancer analysis of such immune-checkpoint inhibitors found that tumor mutational burden (TMB) was correlated with patient outcomes, with high TMB patients (top 10% TMB by histology) having better survival than low TMB patients (bottom 80%; ref. 71). Although immune infiltration data were not presented in this study, it would be of interest to examine links between TMB and immune infiltrate, especially because there appears to be a cutoff (>15%) for TMB to influence the hazard ratio (71). In non–small cell lung cancer, smokers have higher TMB (72) and more PD-1L compared with nonsmokers (73), which can account for superior PD-1 blockade response rates. Our analysis of the TCGA data also found PD-1L expression (CD274) to correlate with predicted immune infiltrate, NK cell score and melanoma patient survival (Fig. 4B), suggesting TMB may be linked to immune infiltration and highlighting the need for further work to examine this interplay.

The NK cell gene signature and NK cell gene score that we describe here can be applied to existing and future cancer data sets becoming available thanks to the efforts of cancer research consortia. The information from such gene signature analyses will allow researchers to stratify responders and nonresponders to conventional treatments, identify patients who are likely to benefit from NK cell–based immunotherapies, and facilitate the development of prognostic markers for personalized immunotherapeutics.

N.D. Huntington reports receiving a commercial research grant from, has received honoraria from speakers bureau of, has ownership interest (including stock, patents, etc.) in, and is a consultant/advisory board member for Biotech. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J. Cursons, F. Souza-Fonseca-Guimaraes, A. Behren, N.D. Huntington, M.J. Davis

Development of methodology: J. Cursons, F. Souza-Fonseca-Guimaraes, M. Foroutan, M.J. Davis

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Cursons, A. Anderson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Cursons, F. Souza-Fonseca-Guimaraes, M. Foroutan, A. Anderson, F. Hollande, S. Hediyeh-Zadeh, A. Behren, N.D. Huntington, M.J. Davis

Writing, review, and/or revision of the manuscript: J. Cursons, F. Souza-Fonseca-Guimaraes, A. Anderson, F. Hollande, A. Behren, N.D. Huntington, M.J. Davis

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J. Cursons, M. Foroutan, N.D. Huntington

Study supervision: J. Cursons, F. Souza-Fonseca-Guimaraes, N.D. Huntington, M.J. Davis

This work is supported in part by project grants from the National Health and Medical Research Council (NHMRC) of Australia (#1147528 to J. Cursons; #1128609 to M.J. Davis; #1124784, #1066770, #1057812, and #1124907 to N.D. Huntington; #1164081 to F. Hollande; and #1140406 to F. Souza-Fonseca-Guimaraes). F. Souza-Fonseca-Guimaraes was supported by NHMRC Early Career Fellowship (1088703), National Breast Cancer Foundation (NBCF) Fellowship (PF-15-008), and grant #1120725 awarded through the Priority-driven Collaborative Cancer Research Scheme and funded by Cure Cancer Australia with the assistance of Cancer Australia. A. Behren is the recipient of a Fellowship from the Victorian Government Department of Health and Human Services acting through the Victorian Cancer Agency. N.D. Huntington is an NHMRC CDF2 Fellow (1124788), a recipient of a Melanoma Research Grant from the Harry J Lloyd Charitable Trust, Melanoma Research Alliance Young Investigator Award, and a CLIP grant from Cancer Research Institute. M.J. Davis was supported by NBCF Career Development Fellowship ECF-14-043 and is the recipient of the Betty Smyth Centenary Fellowship in Bioinformatics. This study was made possible through Victorian State Government Operational Infrastructure Support and Australian Government NHMRC Independent Research Institute Infrastructure Support scheme. Results published here are based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.