A new natural killer cell-specific gene signature predicting recurrence in colorectal cancer patients

The protective role of Natural Killer (NK) cell tumour immunosurveillance has long been recognised in colorectal cancer (CRC). However, as most patients show limited intra-tumoral NK cell infiltration, improving our ability to identify those with high NK cell activity might aid in dissecting the molecular features which underlie NK cell sensitivity. Here, a novel CRC-specific NK cell gene signature that infers NK cell load in primary tissue samples was derived and validated in multiple patient CRC cohorts. In contrast with other NK cell gene signatures that have several overlapping genes across different immune cell types, our NK cell signature has been extensively refined to be specific for CRC-infiltrating NK cells. The specificity of the signature is substantiated in tumour-infiltrating NK cells from primary CRC tumours at the single cell level, and the signature includes genes representative of NK cells of different maturation states, activation status and anatomical origin. Our signature also accurately discriminates murine NK cells, demonstrating the applicability of this geneset when mining datasets generated from preclinical studies. Differential gene expression analysis revealed tumour-intrinsic features associated with NK cell inclusion versus exclusion in CRC patients, with those tumours with predicted high NK activity showing strong evidence of enhanced chemotactic and cytotoxic transcriptional programs. Furthermore, survival modelling indicated that NK signature expression is associated with improved survival outcomes in CRC patients. Thus, scoring CRC samples with this refined NK cell signature might aid in identifying patients with high NK cell activity who could be prime candidates for NK cell directed immunotherapies.

Introduction

Despite improvements in the surgical and medical management late-stage colorectal cancer (CRC), 5-year survival rates for patients with late-stage colorectal cancer (CRC) remain extremely poor and innovative treatment strategies are needed. Although T-cell directed immunotherapies are strikingly effective in several solid cancer types, durable responses are limited to colorectal tumours exhibiting defective mismatch repair processes (dMMR), thus benefiting approximately 5% of all CRC patients (1). Downregulation of major histocompatibility complex class I (MHC-I), whose expression is a pre-requisite for T-cell-mediated immune killing, has been reported as a major cause of secondary resistance to checkpoint blockade therapy (2, 3). Thus, there is an urgent need to explore alternate immunotherapeutic strategies which harness other antigen-independent cell types.

In recent years, natural killer (NK) cells have emerged as a promising candidate for immunotherapeutic development. Adoptive transfer of primary NK cells suppresses T-cell mediated graft-versus-host disease and exacerbates graft-versus-tumour responses (4, 5). Similarly, a myriad of alloreactive NK cell lines - primarily NK92, but also NKL, KHYG01 and YTS – have demonstrated safety and efficacy in HLA-mismatched recipients in early-stage clinical trials (reviewed in (6)). This suggests that NK cell therapeutics may be harnessed in an “off-the-shelf” manner in the future, circumventing the difficulty associated with the ex vivo expansion of antigen-specific T-cell clones. Likewise, chimeric antigen receptor (CAR) NK cells have shown promising results in xenograft models (7, 8), and monoclonal antibodies which neutralise NK cell inhibitory receptors such as the KIR family (9, 10) and NKG2A (11) have entered early-stage clinical trials.

NK cell activity is inversely correlated with cancer incidence (12) and there is a wealth of evidence supporting the role of NK cells in controlling both spontaneous and experimental metastasis (13, 14). In CRC, NK cell infiltration has been identified as a positive prognostic marker in both primary (15, 16) and metastatic (16, 17) disease. NK cells differ from conventional lymphocytes in that they function in an MHC-I unrestricted manner in accordance with the “missing-self” hypothesis (18, 19). In this manner, NK cell-directed immunotherapies may overcome the restricted benefit of antigen-specific T-cell responses in tumours with high mutational diversity. Moreover, harnessing NK cell cytotoxicity is a promising opportunity in the treatment of immunologically “cold” tumours such as CRC that undergo loss of MHC-I expression. Indeed, recent studies have reported that 40% of patient-derived CRC organoids exhibit MHC-I loss which could not be rescued by IFN stimulation (20), concordant with clinical reports of MHC-I loss in approximately 60% of MSI CRCs (21). Yet, clinical enthusiasm has been tempered due to the scant NK load reported in most CRC tumors, despite high levels of chemokines and cytokines (22). Thus, novel means of determining which patients show high NK cell infiltration and activity and might therefore benefit from NK cell directed immunotherapies is needed.

Transcriptomic signatures are sets of genes whose coordinated expression has a verified association with specific biological parameters – such as cell type and phenotypic state – or clinical measures, including disease subtype, survival outcome or therapeutic response. Deconvoluting the relative abundance of immune cell subsets using transcriptomic signatures offers many benefits over the comparatively low-throughput methods of immunophenotyping which are currently available. In CRC alone, a multitude of prognostic signatures have been reported (23, 24) alongside several signatures which predict response to 5-FU based chemotherapies in Stage II-III disease (25–27). Additionally, RNAseq data from a cohort of 40 CRC patients with metastatic or relapsed disease was used to derive a 27-gene signature able to discriminate responders and non-responders to the FOLFOX6 chemotherapy regimen with accuracy of 92.5%, demonstrating the powerful role which signature analysis may play in personalising patient care (28).

In their simultaneous assessment of multiple markers, transcriptome-wide approaches circumvent the limitations of single-marker or low-throughput phenotypic assessments. This may also allow for previously unidentified markers or associations to be identified which may provide important insights into immune cell biology. However, a major consideration when employing such computational approaches is the suitability of the reference profiles, as the gene signatures derived from sorted cell types in healthy individuals may not accurately reflect those of the potentially dysregulated cells in a diseased individual. This is particularly pertinent in the case of NK cells, as activated NK cells have unique transcriptional profiles as compared with resting NKs from healthy volunteers (29). Additionally, as signatures often rely on highly expressed genes and thus often incorporate genes that are not that specific of a given cell type, the accuracy of cell type detemination may be compromised. Although prognostic NK cell gene signatures exist for highly immunogenic tumour types such as melanoma (30) and renal cell carcinoma (31), it is currently unclear how appropriate these are for use in tumours such as CRC which traditionally show a poor immune cell infiltration. Other works which have looked at the transcriptomic profiles of lymphocytes in CRC have focussed on defining residency versus exhaustion states (32). where the genesets used to define CD4+ and CD8+ T-cells had several genes overlapping with the NK cell geneset.

Here, we present a novel NK cell transcriptomic signature which can be used to infer NK cell abundance from bulk RNA sequencing data of primary CRC samples. Candidate NK cell-related genes were pooled from previously published works and in-house differential expression (DE) analyses and sequentially filtered to ensure their fidelity as NK cell markers with minimal off-target expression in tumour, stromal and other immune cell types. Single-cell RNAseq (scRNAseq) data from primary CRC samples were then used to validate each signature gene in tumour-infiltrating NK cells. We then show that high NK cell score is associated with upregulation of cytolytic and chemotactic transcriptional processes, and survival analysis revealed that patients with higher evidence of NK cell activity demonstrate significantly longer recurrence and disease-free intervals. Collectively, the NK cell signature allows for the identification of CRC patients with high NK cell activity, which may aid in defining the molecular characteristics associated with strong response to NK cell targeting immunotherapies.

Methods

Preparation of publicly available RNAseq data

Datasets used in the present study are listed in Table S4. Raw counts files from publicly available datasets were downloaded rom Gene Expression Omnibus (GEO) using the NCBI portal (http://www.ncbi.nlm.nih.gov/geo/). CCLE data was downloaded as a PharmacoSet (PSet) through the PharmacoGx R/Bioconductor package (version 1.6.1). For in-house data, 3’ RNA-seq reads were aligned using HISAT2 against human genome GRCh37 (release 75) and the featureCounts tool from the RSubread package was used to quantify the number of reads for each gene per sample. The filterByExpr function from the edgeR package (v 3.28.1) was used to filter lowly expressed genes and calculated count- or transcript-per-million (CPM/TPM) values. The MyGeneset application of the ImmGen online databrower (http://www.immgen.org/Databrowser19/DatabrowserPage.html) was used for analysis of GSE15907.

The GSE107011 dataset is composed of 29 FACS-isolated cell types from the blood of healthy individuals. For this project, The NK cell, basophils, and neutrophil samples were annotated as per the authors’ annotation file. However, to circumvent the known issues associated with multiple comparisons, and to facilitate comparisons with data sets offering less cellular subtype resolution, the various maturation and functional states for other cell types were merged to form larger groups: naïve, switched, non-switched and exhausted B-cell samples were merged and annotated as “B-cells” (n_samples = 16); naïve, central memory, effector memory and terminal effector CD8+ T-cell samples were merged and annotated as “CD8+ T-cells (n_samples = 16); classical, intermediate and non-classical monocytes were merged and annotated as “monocytes” (n_samples = 12); myeloid and plasmacytoic dendritic cell samples were merged and annotated as “dendritic cells” (n_samples = 8). Plasma and progenitor cells were excluded from downstream analyses. Given the transcriptional overlap between NK- and T-cells, T- cell subsets (namely CD8+, CD4+, MAIT cells, Tfh, TH1, TH17, Th2, Treg, γδ+ and γδ -) were interrogated both individually and as part of a larger, merged group annotated as “pan-T-cells”. Resultingly, 17 discrete immune cell subtypes from this dataset were interrogated. For GSE60424, only blood samples from healthy individuals were included.

Differential expression analysis, gene set testing and survival modelling

The singscore (v1.6.0) R package was used for Single-sample gene set enrichment analysis against various molecular signatures. Depending on the direction of gene sets, we used different settings of simpleScore function as specified in the documentation. NK-high and NK-low groupings were defined as the top and bottom 10% of samples, respectively, when ranked by NK scores. DE analysis between NK-high and NK-low samples was performed using the voom-limma (v 3.42.2) pipeline (Law et al., 2014). After running eBayes, we considered genes with absolute log2FC > 1 (for ovexpression) or log2FC < −1 (for repression) and adjusted p-value < 0.05 as DEGs. The goana function from the limma package was used to perform gene ontology (GO) analysis and camera gene set testing was performed using MSigDB signatures retrieved from the WEHI bioinformatics portal (http://bioinf.wehi.edu.au/software/MSigDB/). Survival modelling was performed using the survminer (v0.4.8) and survival (v3.2-7) packages using the clinical annotation files provided from the sources listed in Table S4.

Data wrangling and visualization

All computational analyses were performed using R (version 3.6.1). For data wrangling and visualization, base R functions were using alongside several core packages from the tidyverse (v 1.3.0) R package. tidyr (v 1.1.2) and dplyr (v 1.0.2) were used for reading and manipulating the data, as well as ggplot2 (v 3.2.1) and cowplot (v 1.0.0) for data visualization. Heatmaps were generated using the complexHeatmap (v 2.2.0) or pheatmap (v 1.0.12) R packages.

Statistical analysis

All statistical analyses were performed using R. Data is expressed as mean ± standard deviation (SD) unless otherwise indicated. The minimum threshold for rejecting the null hypothesis was p<0.05. For results where statistics are shown, significance is denoted as: * = p<0.05; ** = p< 0.01; *** = p<0.001.

Code availability

The code used throughout this study is available on Github (https://github.com/cshembrey/NK_Signature_CRC).

Results

Collation of candidate NK cell signature genes

To identify a CRC-specific gene expression signature associated with NK cell abundance, we implemented a novel pipeline that involved curation and sequential refinement of putative NK cell genes against immune, tumour and stromal cells in multiple bulk and single cell RNAseq datasets (workflow outlined in Figure S1). Firstly, 605 unique genes were collated from eight partially overlapping sources (Figure 1). Four previously curated NK cell signatures were compiled: “CIBERSORT Active” (n_genes = 56) and “CIBERSORT Resting” (ngenes = 56) refer to the gene sets corresponding to activated and resting NK cells, respectively, as reported in the LM22 signature matrix used for the CIBERSORTx algorithm (33); The “Cursons Extended NK cell Geneset” gene set (n_genes = 112 genes) was derived from the supplementary data table of a melanoma-specific NK cell signature previously reported by Cursons and colleagues (30); The “Wang NK cell marker” gene set (ngenes = 13; 34) is composed of markers used to guide the immunophenotyping of different cell subsets in bulk RNAseq data from CRC cell lines and primary samples, and the “receptors” gene set (n_genes = 43) was compiled by mining the literature for various receptor subsets (eg. activating, inhibitory, chemokine or cytokine receptors) with documented expression on NK cells.

FIGURE 1

Figure 1 Collation of putative NK cell signature genes. UpSet plot of gene intersects from eight partially overlapping sources (listed at right). Gene set (blue bars) and intersection (red bars) sizes are indicated. Inset: Venn diagram of DEGs identified from pairwise comparisons of NK cells versus other immune cell types in GSE60424 and GSE107011. Only the “Union DEGs” (n_genes = 194) were retained as candidate signature genes.

Subsequently, three novel gene sets were derived from the results of differential expression (DE) analysis, where putative genes enriched in NK cells were identified by pairwise comparison of NK cells and at least one other immune cell type. DE analysis of the GSE60424 bulk RNAseq dataset (35), composed of six immune FACS-isolated cell types (NK cells, B-cells, CD4+ T-cells, CD8+ T-cells, monocytes and neutrophils), identified 280 DEGs (“GSE60424 DEGs”). Similarly, DE analysis of GSE107011 (36), an RNA-seq dataset composed of 29 FACS-isolated immune cell types, yielded 427 DEGs (“GSE107011 DEGs”). Finally, to enhance our resolution when discriminating between NK and T-cells, a subsequent “pan-T DEGs” (ngenes = 23; from GSE107011) was constructed from the DE genes when only NK versus T-cell comparisons were considered; here, the selected genes were those upregulated in NK cells relative to all T-cell subsets (eg. upregulated in NKs relative to CD8+ and CD4+ and MAIT etc.), irrespective of their expression level in other innate immune cell types.

Following compilation, 605 unique candidate NK genes were identified. To examine the representation of traditional NK cell markers versus potentially novel NK cell-related genes, the interconnectedness of this gene set was compared across the eight sources (Figure 1). Most genes identified (451/608; 74%) were uniquely found from our DE analyses. Surprisingly, there were zero genes which were conserved across all eight sources and only two genes, KLRD1 and KIR3DL2, were conserved across 7/8 sources. We hypothesize that the lack of consensus between these datasets may reflect the multiple discrepancies between studies in terms of criteria chosen for candidate gene selection and filtering, combined with variations in the NK subtypes present in each dataset.

To further refine the comparatively large GSE60424 and GSE107011 DEG sets, a “Union DEGs” gene set was created (Figure 1; inset). For inclusion in this merged gene set, a given candidate gene needed to be differentially upregulated in NK cells relative to all other immune cell types in one or both data sets (eg. upregulated in NK cells versus T-cells and B-cells and neutrophils) differing from the previous DEG analyses where the gene in question need only be upregulated in NK cells on a pairwise basis (eg. upregulated in NK cells versus T-cells, but not necessarily B-cells nor neutrophils). Of the 525 unique DEGs identified across GSE60424 and GSE107011, 194 genes fit this criterion. After reintegration of these “Union DEGs” with the genes derived from other sources, 295 candidate genes remained.

Although DE analysis identifies genes which are preferentially expressed by NK cells as compared with other cell types, the expression of individual DEGs by NK cells may still be very low. This is problematic when aiming to identify NK cell signals from tumour sequencing data, particularly given that NK cells represent a small proportion of the total cell number. Thus, we further refined out geneset by retaining only those candidates with higher median in NK cells relative to other immune cell subsets in GSE107011 (Figure 2A; Figure S2A; purple boxes) and GSE60424 (Figure S2B; purple boxes). Sets of passing genes for each dataset were derived by taking the intersect of the specific genes (n_genes = 49) from GSE107011 (Figure 2A) and GSE60424 (Figure S2B). Of these, 10 genes (ADGRB2, B4GALT6, LDB2, LIM2, LINC01451, LRRC43, PCDH1, PRSS57, RAMP1 and RNF165) were identified in both data sets.

FIGURE 2

Figure 2 NK cell specificity filtration against immune and CRC-supportive cell types. (A) Biplots depicting median expression for each candidate gene (rings; coloured by sum of sources) in NK cells from GSE107011 (logTPM) versus other immune cell types. The intersect of passing genes for each pairwise comparison (purple boxes) were retained as candidate genes. (B) Biplots depicting median expression for each candidate gene (rings; coloured by sum of sources) in NK cells from GSE107011 (logTPM) versus CCLE CRC cell lines (left column; logRPKM) and GSE90830 (CRC cell lines; right column; logRPKM). The union of failing genes for each pairwise comparison (blue boxes; n_genes = 12) were flagged for removal from the candidate geneset. (C) Lineplot of candidate gene expression in the stromal cell (CD31+; gold), fibroblast (FAP+; green) and epithelial cell (EpCAM+; light blue) compartments of six primary CRC samples (GSE39396), normalised to leukocyte (CD45+; purple) expression. Genes with normalised expression >1 were flagged for removal. CCLE, Cancer Cell Line Encyclopedia; CoT, CRC primary tumour; LT, CRC liver metastasis.

Interestingly, this analysis highlighted that many of the “classical” NK cell markers that are in multiple sources (Figure 2A; yellow/green rings) are also expressed at very high levels in “unconventional” T-cell subsets such as MAIT, γδ+ and γδ- T-cells. As illustrated by GSE60424 (Figure S2B), where T-cells are exclusively grouped as CD8+ or CD4+, many FACS-based studies do not include these relatively niche T-cell subsets, possibly leading to the identification of putative NK cell-specific genes which are in fact highly expressed by unconventional T-cell subsets.

The two sorted NK cell-containing data sets interrogated thus far (GSE60424 and GSE107011) were derived from the peripheral blood of healthy individuals; this approach has limitations, as the gene expression profiles of such cells may not necessarily reflect those of tissue-infiltrating NK cells nor NK cells in the context of cancer. To address this concern, we performed an independent DE analysis (see methods for details) on an NK-cell containing scRNAseq data set (GSE146771 (37), composed of SMART-Seq2 and 10X subseries) generated from a cohort of CRC patients (n_samples = 20; from 18 unique patients) with Stage II-III disease of varying pathological grade, MSI status, and extent of nodal involvement. This analysis identified several additional NK cell marker genes, and allowed us to cross-validate the 49 genes derived from bulk RNA-seq analysis (mentioned above) at the single cell level. Following these analysis steps, 82 candidate genes were prioritised.

Signature gene specificity filtration against tumour and stromal cells

As the NK cell signature is designed to resolve the NK cell fraction from bulk RNAseq data of tumour tissue, it is imperative that the genes in the signature should not be expressed by the tumour cells themselves. To determine whether any of the 82 identified genes were expressed by tumour cells, additional specificity thresholding was performed against two sources of CRC cells (Figure 2B). As bulk RNAseq from CRC tumour tissue would be expected to contain immune cell fractions, two cohorts of CRC cell lines were used in this analysis: The Cancer Cell Line Encyclopedia (38) (n_samples = 57) and GSE90830 [nsamples = 44 (39)]. These datasets revealed that 12 genes had relatively higher expression in CRC cells compared with NK cells (Figure 2B; blue boxes), and they were therefore excluded from further analyses.

Despite the immune- and tumour-specific filtration steps previously performed, the possibility of contaminant expression of our candidate genes by other non-NK cell types such as including fibroblasts, stromal cells and non-transformed epithelial cells has not been accounted for. To address this, the relative expression of our candidate genes was interrogated in the leukocyte (CD45+), stromal cell (CD31+) and epithelial cell (EpCAM+) fractions isolated from the tumours of 6 patients with CRC (GSE39397 (39); as well as against cultured, normal colon mucosa-derived fibroblasts (CCD-Co-18). Failing genes (ngenes = 14) were defined as those with significantly higher expression relative to leukocytes in a non-leukocyte subset (Figure S3; red headers, summarised in Figure 2C).

Interrogating expression of signature genes in tumour-infiltrating NK cells

Having validated the specificity of our candidate genes at the immune, tumour and stromal cell level in bulk RNAseq, we returned to GSE146771 to confirm this specificity in scRNAseq data. By interrogating the gene expression distribution across Uniform Manifold Approximation and Projection (UMAP) plots (Figure 3), the specificity of each gene for particular NK cell subsets was elucidated. For example, in the SMART-Seq dataset three NK cell subgroups could be discerned based upon marker expression, providing greater cellular resolution: CD16+, a marker expressed on virtually all CD56^dim NK cells; GZMK+, a cytotoxicity marker frequently associated with the CD56^bright subpopulation; and CD103+, a marker of tissue residency. Certain “NK cell specific” genes such as SH2D1B and KIR2DL3 showed equivalent expression in each of the three NK subsets, whereas other genes were subtype selective. For example, LINGO2, PRSS57 and RNF165 were preferentially expressed by the CD16+ subgroup, whereas WIPF3 expression was high across both the CD103+ and GZMK+ subgroups but sparse in the CD16+ population.

FIGURE 3

Figure 3 Candidate gene expression in CRC-infiltrating immune cells. UMAP plots of dissociated primary CRC samples from (A) GSE146771 (SMART-Seq2 scRNAseq; n_samples = 10) and (B) GSE146771 (10X scRNAseq; n_samples = 10) coloured by cell type (at left) and candidate gene expression. Maximum expression (LogTPM) is indicated in parentheses above each plot.

This approach also allowed us to identify genes which, although differentially upregulated in NK cells, showed moderate basal expression across multiple cell types (Figures S4, S5; “High DEGs”; pale blue headers) as well as genes whose expression was no longer specific to NK cells once evaluated at the single-cell level (Figures S4, S5; “Non-specific”; red headers). For example, many candidate genes - including PLAC8, SLC15A4, SLFN13 and ST8SIA6 – have high expression in NK cells although they are promiscuously expressed by multiple cell types, with particularly high expression in the myeloid subset, warranting their exclusion from the final gene list.

Validating the subtype specificity of the curated NK cell signature

Having refined the NK cell signature using CRC-infiltrating NK cells, we next interrogated whether the signature exhibited any biases towards NK cells from particular subsets or sources. GSE133383 (40) contains transcriptomic data for both the immature CD56^bright and mature CD56^dim subsets of NK cells isolated from the blood, lymphoid organs (spleen, lymph nodes and bone marrow) and lungs of four healthy donors. As is the case with most transcriptional analyses of NK cells, when interrogated in GSE133383 our signature genes cluster according to NK cell maturation state rather than tissue source (Figure 4A). Whilst a substantial number of genes which appear to have relatively uniform expression across both the CD56 bright and dim subsets (Figure 4A), multiple subtype-enriched markers have also been retained. Enriched markers for the populous CD56^dim class include the KIR family members, GZMB (40–42) and the chemokine receptors CXCR1 and CXCR2 (43). Two recently identified markers of terminally mature NK cells, HAVCR2 and CX3CR1 (44) show preferential expression in the CD56^dim subset as expected. Analogously, enriched markers for the CD56^bright subset include CD56 (NCAM1) itself, XCL1 (44) and KLRC1 (45), whose inclusion indicate that the signature has sufficient resolution to detect this relatively minor population. Collectively, this demonstrates that our signature captures both general and subset-enriched NK cell markers, suggesting that it is representative of all types of NK cell.

FIGURE 4

Figure 4 Profiling refined gene set expression in NK cell subsets. (A) Heatmap of refined gene set expression in CD56^bright (orange) and CD56^dim (purple) NK cells isolated from the blood (red), spleen (yellow), bone marrow (blue), lung (pink) and lymph node (green) of four healthy donors (GSE133383). (B) Correlation matrix of refined gene set. Genes are annotated (at left) according to whether they show preferential expression in the CD56^dim (orange), CD56^bright (purple) or neither (green) NK cell population (C) Boxplot expression of selected genes in CD56^dim (orange boxplots) vs CD56^bright (purple boxplots) NK cell populations. Headers are coloured based on whether the gene is preferentially expressed in CD56^dim (orange headers) vs CD56^bright (purple headers) NK cells or shows equivalent expression across both subsets (green headers). Significance was assessed with Student’s T-test; ***p-value < 0.001. ; ****p<0.0001; ns: non-significant.

Many of the genes which are lowly-expressed in GSE133383 (Figure 4A; LIM2, SPTSSB, LGALS9B, PRSS57, ADGRG3, and, to a lesser extent, DRAXIN), which lacks intestinal NK cells samples, are highly expressed in particular subsets of CRC-infiltrating NK cells (Figures 4A, B), possibly reflecting a role for these genes as CRC-specific NK cell markers induced by the tumour microenvironment. It is noteworthy that the peripheral blood CD56^bright subset – the most abundant NK cell subset in humans - do not cluster together, supporting the idea that this signature is more CRC-specific than previously reported NK cell signatures derived from blood-based analyses.

Next, as it is expected that genes which are considered robust markers of a given cell type should have coordinated expression patterns, the correlation between the individual genes in the signature was assessed. Correlation analysis demonstrated that the signature genes cluster into two major blocks where expression is highly cross-correlated, representing the genes selective for CD56^dim or CD56^bright NK cells (Figure 4B; orange and purple annotations at left). Subtype selectivity for each gene was defined based on its relative expression in each of these subsets (Figure 4C and Figure S6). That the two major blocks are anti-correlated may reflect the progressive loss of the CD56^bright markers as NK cells mature and reinforces the idea that our signature allows for pan-NK cell rather than subset-specific detection. Situated between the two major clusters are a subset of genes (including GNLY, TXK and PTGDR) which appear to have equivalent expression in CD56^bright and CD56^dim cells (Figure 4B; green annotations), likely corresponding to a “transitional” NK cell phenotype reported by several groups (42–44, 46).

As a final visualisation measure, the expression of the NK cell signature in aggregate (ie. profiling expression of the signature as a whole, rather than the expression of each individual gene) was performed in an independent human dataset, GSE22886 (47). This dataset is composed of twelve different leukocyte subsets isolated from the PBMCs of healthy donors and was used to avoid the biases of testing our signature on the GSE60424 and GSE107011 data from which it was partially derived. In GSE22886, our NK cell signature was clearly enriched in the NK cell samples relative to the other cell types (Figure 5A). Due to the high sequence and transcriptomic homology between human and murine NK cells (48), we next interrogated the performance of the human NK cell signature in GSE15907 (Immgen) (49); a microarray dataset generated following the ex vivo isolation of multiple immune lineages from adult B6 mice (Figure 5B). Notably, in the murine context our NK cell signature clearly discriminates NK cells from CD8+ T-cells, γδ-T-cells and NKT cells, highlighting the flexibility of this novel NK cell signature to be used in in vivo studies to detect either endogenous murine or xeno-transplanted human NK cells.

FIGURE 5

Figure 5 NK cell signature efficiently discriminates NK cells from other haematopoietic compartments. Boxplots of NK cell signature expression in sorted peripheral blood cells from (A) GSE22886 (Log2) independent human PBMC data and (B) GSE15907 (Immgen) murine data. Each boxplot represents one sample (coloured by cell type, at right) and individual points are single signature genes.

Differential gene expression analysis confirms that cytotoxic and migratory programs are associated with high NK score

Having finalised the NK cell signature (Table S1, n_genes = 43), we next sought to determine which genes were concomitantly up- or downregulated in the samples whose transcriptomic profiles showed strong evidence of this signature. To identify these samples, we used a single-sample, rank-based gene-set scoring method termed singscore (50) to score samples against our NK signature. Here, a sample with a high NK score is interpreted as having high evidence of NK cell activity, whereas a sample with low NK score exhibits limited NK cell signature expression. Using samples sourced from two large, publicly available repositories of primary CRCs, the TCGA colorectal adenocarcinomas (TCGA-COAD; n_samples = 454), and GSE39582 (n_samples = 566 (51), we defined the NK-high and NK-low groupings based on the top 10% and bottom 10% of scored samples, respectively. We then performed differential expression analysis to compare the gene expression profiles (GEPs) of samples with high scores to those with low scores in each of the two data sets.

A strong positive correlation was observed between the logFC of genes in TCGA and GSE39582 datasets (Figures 6A, B; Spearman correlation R = 0.72, p-value < 2.2 x 10^-16), suggesting similar transcriptional programs of NK inclusion/exclusion in both data sets.

FIGURE 6

Figure 6 NK score is associated with high chemokine and cytolytic activity. (A) Scatterplot of the LogFCs of genes comparing NK-high and NK-low groups in TCGA and GSE39582 datasets; Genes that were differentially upregulated (red points) or repressed (blue points) in both datasets are highlighted. DEGs were defined as those with adjusted p-value < 0.05 and absolute LogFC > 1 or LogFC < -1. (B) Venn diagram of upregulated (upper) and down-regulated (lower) DEGs identified in (A). (C) Barcode plots showing enrichment of genes in the “NK cell mediated immunity” gene set (GO:0002228) in NK-high samples in both data sets (D) PCA plot of top 50 most significant GO terms in the NK-high group, distributed by semantic similarity of genes within each term (E) STRING network analysis of protein-protein interactions between the overlapping upregulated DEGs from (A), coloured by biological process. Thickness of the connecting line indicates the strength of evidence for the predicted interaction. GO, gene ontology.

48 genes were differentially upregulated in NK-high samples from both datasets (Figure 6B, upper panel, Table S2). Of these, there was a strong enrichment of genes encoding cytotoxic effectors which are critical for NK cell killing, including the granule proteins NKG7 and granulysin (GNLY), as well as multiple members of the granzyme family including GZMB, GZMA and GZMH. Additionally, there was an overrepresentation of genes encoding ligands for chemokines implicated in NK cell trafficking, including CXCL9, CXCL10 and CCL5 (RANTES). Importantly, these ligands are expressed by tumour cells rather than the NK cells, suggesting that the high NK cell density in these tumours is at least partly driven by a tumour-intrinsic factor.

Conversely, 45 genes conserved between the two data sets were differentially repressed in the NK-high group (Table S3). The parallel downregulation of multiple genes which have been associated with increased metastasis and poor prognosis in CRC was observed, including Dachshund family transcription factor 1 (DACH1) (52) and the metalloprotease meprin-α (MEP1A) (53) and CXCL14 (54). Interestingly, recent evidence indicates that CXCL14 is essential for MHC-I upregulation (55) and as such, it stands to reason that CRC samples with high NK cell activity show decreased CXCL14 expression.

To identify biological processes associated with the NK-high phenotype, gene ontology (GO) enrichment analysis was conducted. Competitive gene set testing against the whole transcriptome of NK-high vs NK-low samples confirmed significant enrichment of the GO term “NK cell mediated immunity” in the NK-high group (Figure 6C). Moreover, the top 50 most significantly enriched GO terms in the NK-high group, when clustered by semantic similarity, converged on umbrella terms including “immune response”, “leukocyte activation”, “leukocyte cell-cell adhesion” and “cytokine production” (Figure 6D).

To visualise the functional synergy of the 45 overlapping upregulated DEGs (from Figure 6B) which define the NK-high group, the Search Tool for Retrieval of Interacting Genes/Proteins (STRING) platform was used to construct a protein-protein interaction (PPI) network (Figure 6E). Consistent with DE results, the PPI network was primarily centred around CC- and CXC-family chemokine ligands (CCL5, CCL4, CCL8, CXCL9, CXCL10, CXCL11, CXCL13, CXCL18) as well as cytolytic effectors (GNLY, GZMH, GZMK, GZMB, GZMK, NKG7, PRF1). A minor node related to antigen processing and presentation (HLA-DMA, HLA-DMB, HLA-DPA1, TAP1, CD74) may point towards a concomitant activation of effector T-cells in the NK-high group; however, chronic antigen stimulation is also known to induce lymphocyte exhaustion, resulting in reduced cytotoxic function (56).

Collectively, these data indicate that the NK-high group is defined by strong chemotactic signalling and high cytolytic activity. As these genes are signposts of an immune-active microenvironment, their upregulation in the NK-high group suggests that the NK cell signature selects for functionally active NK cells, rather than merely NK cell presence.

High NK score is associated with improved survival outcomes and other clinical parameters

Given that NK cell load has been associated with better patient prognosis in CRC (15, 17), we next performed Kaplan-Meier survival analysis to determine whether signature expression was associated with survival probability. Due to significant survival differences between patients with a primary or metastatic disease, we focussed only on those patients with Stage I-III CRC. For all survival analyses, patients were stratified by NK score where “NK-High” and “NK-Low” are defined as samples above and below the median NK score, respectively.

For the TCGA-COAD cohort, the NK-High group had significantly increased disease-free interval (DFI; defined as the period from date of diagnosis until a tumour progression event e.g. locoregional recurrence or distant metastasis) as compared with the NK-Low group (Figure 7A, log rank p-value = 0.0054, p-value from multivariate Cox regression model [accounting for age, stage and MMR status] = 0.02). Similarly, the GSE39582 cohort exhibited a trend towards significantly prolonged recurrence-free survival (RFS), defined as the period between surgical resection and a tumour progression event, at the univariate level (Figure 7B, log rank p-value = 0.053). There were no significant differences in overall survival (OS; Figures S7A, C) nor progression-free interval (PFI; Figure S7B) between the NK-High and NK-Low groups in either dataset.

FIGURE 7

Figure 7 NK score is associated with survival outcome and other clinical parameters. Kaplan-Meier survival curves for patients stratified by NK score (where “NK-High” and “NK-Low” are defined as samples above and below the median NK score, respectively) for Stage I-III CRC patients in the (A) TCGA-COAD (DFI) and (B) GSE39582 (RFS) cohorts. Survival differences were tested using both log-rank and multivariate Cox proportional hazards models (adjusted for age, tumour stage and MMR status) with corresponding p-values indicated. (C) Boxplots of association of NK score with CMS status and (D) mutational load in TCGA-COAD patients. (E) Boxplots of association of NK score with MSI (F) CIMP (G) BRAF and (H) KRAS status in GSE39582 patients. (Student’s T-test; **p-value < 0.01; ***p-value < 0.001; ****p < 0.0001). DFI, disease-free interval; RFS, recurrence-free survival, MMR, mismatch repair; CIMP, CpG island methylator phenotype.

We next interrogated the relationship between NK score and various clinical and molecular parameters such as patient history, molecular subtype, and driver mutation status (see additional analyses in Supplementary Data). In both the GSE29582 (Figure S8A) and TCGA (Figures S9B, C) cohorts, NK scores were increased in early-stage disease (Stage I-II) as compared with late-stage disease (Stages III-IV). In the TCGA data, NK scores were considerably enriched in the immunogenic CMS1 subtype (Figure 7C) and in those tumours with higher mutational load (Figure 7D). In GSE29582, high NK score was significantly increased in those patients with MSI (Figure 7E) and CpG Island Methylator Phenotype (CIMP+; Figure 7F) disease. This was corroborated in the TCGA data (Figures S9A, B), where high NK score was also significantly enriched in those with MSI-associated clinical parameters such as tumour hypermutation (Figure S9D) and MLH1 silencing (Figure S9E).

With respect to CRC driver mutations, NK score was significantly enriched in those patients with BRAF mutant (Figure 7G) or KRAS wildtype (Figure 7H) genotypes, although there was no association with TP53 status (Figure S8D). There was no association between NK score and tumour characteristics such as tumour location (Figure S8B; distal versus proximal colon) and histological subtype (Figure S9F; mucinous vs non-mucinous). Similarly, no association was found between NK score and factors related to clinical history including adjuvant chemotherapy status (Figure S8C), prior CRC diagnosis (Figure S9G) or evidence of synchronous disease (Figure S9H). In sum, these results support the clinical utility of using this newly derived NK cell signature in the context of CRC.

Discussion

Bioinformatic approaches which allow for the deconvolution of immune cell subsets from bulk sequencing data have revolutionised our ability to assess the immune landscape of individual tumours. Transcriptomic signatures which predict intra-tumoral NK cell infiltration have been shown to indicate improved patient survival in various cancers (30, 31). However, despite accumulating evidence that NK cell load is prognostic in CRC, there are currently no extensively curated NK cell signatures explicitly designed for use in the context of CRC.

Here, through extensive computational curation of established NK cell-related genes with putative markers discovered through DE analysis of various bulk RNAseq and scRNAseq datasets, we define a comprehensive NK cell signature specific for CRC samples. Survival modelling in two large cohorts of primary CRC patients indicated that NK signature expression is associated with prolonged progression- and recurrence-free intervals, consistent with previous reports on the beneficial prognostic impact of NK cell infiltration in other solid cancers (57). That these metrics are associated with tumour progression (rather than overall survival) may relate to the fact that NK cells are believed to play a greater role in the suppression of metastases rather than the prevention of primary tumours (58, 59).

Although the NK cell signature presented herein is not the first of its kind, we believe that it is the most specific and the most extensively curated NK cell signature for the precise identification of the wide range of NK cell subsets in CRC samples due to its derivation in the context of this tumour type. Practical limitations have meant that most NK cell gene signatures, including that employed in the widely used CIBERSORT immune cell deconvolution tool, were curated using NK cells isolated from the peripheral blood. As NK cells exhibit tissue-specific phenotypes (40, 42, 60), and the tumour microenvironment is known to rewire NK cell transcriptional programs (60–62), it is difficult to assess how accurately these pre-existing NK cell signatures perform when extrapolated for use in the context of cancer. Moreover, striking differences in the transcriptomic profiles of circulating NK cells versus tumour-infiltrating NK cells from the same patients have also been reported (62), emphasising the importance of using scRNAseq data from tumour-infiltrating NK cells to faithfully derive a reference transcriptomic profile for CRC-associated NK cells.

Another issue facing pre-existing NK cell signatures has been a potential lack of specificity. For example, although a four-gene signature (composed of NCR1, PRF1, CX3CL1 and CX3CR1) was shown to accurately distinguish NK-high vs NK-low subgroups in clear cell renal cell carcinoma (31), these genes are ubiquitously expressed by many immune cell types. More broadly, the transcriptional programs of NK cells are highly overlapping with those of multiple T-cell subsets, particularly γδ-T-cells (63). Numerous studies have demonstrated that several archetypical NK cell markers including NKG2A (64), NCRs (65, 66) and various KIRs (67) are expressed by both traditional αβ- and unconventional γδ-T-cells. Our strategy of prioritising genes based on their NK-cell specificity rather than only those with the highest expression may therefore increase precision when teasing apart the NK cell contribution from complex mixtures of other immune cells.

In both the TCGA and GSE39582 cohorts investigated, NK scores were higher in patients with early-stage disease, aligning with results of a large, pan-cancer meta-analysis reporting that NK cell infiltration is lower in advanced-stage tumours (57). Similarly, high NK score was associated with several clinically useful molecular parameters such as MSI disease, CIMP positive status and BRAF mutation which have each been previously linked with high TIL and NK cell infiltration (68, 69). These features are all defining characteristics of the highly immunogenic CMS1 molecular subtype of CRC (70). Moreover, DE analysis between the NK high- and low-scoring patients identified transcriptional and biological processes associated with high cytotoxic and chemokine activity in samples with high NK scores. As the signature efficiently discriminates murine NK cells from other immune cell types, this signature may also prove useful in prospective in vivo studies.

It remains unclear how optimally transcriptional signatures validated in a particular cancer will perform in alternate tumour types. Thus, future work may focus on cross-validating NK cell signatures in contexts outside their tumour-of-derivation, and/or on defining a pan-cancer NK cell signature which is less sensitive to subtype-specific influences. Prospective studies may also aim to dissect the transcriptional profile of intra-metastatic NK cells and determine whether, or in what ways, this diverges from that of primary CRC-infiltrating NK cells.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

CS, FH, and MF contributed to conception and design of the study. FH provided funding to support all of the study. CS sourced data and performed all computational analyses, with guidance from FH and MF. MF wrote sections of the codes used in the analysis. FH and MF supervised CS throughout the study. CS wrote the first draft of the manuscript. MF wrote sections of the manuscript. FH contributed to manuscript revision. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by The National Health and Medical Research Council of Australia (project grant 1164081, FH) and by the Tour de Cure Foundation (Senior research grant, FH).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.1011247/full#supplementary-material

References

1. Birenda K, Hwang JJ. Advances in immunotherapy in the treatment of colorectal cancer. Am J Hematol/Oncol (2017) 13:4–8.

Google Scholar

2. Rodig SJ, Gusenleitner D, Jackson DG, Gjini E, Giobbie-Hurder A, Jin C, et al. MHC proteins confer differential sensitivity to CTLA-4 and PD-1 blockade in untreated metastatic melanoma. Sci Trans Med (2018) 10:3342. doi: 10.1126/scitranslmed.aar3342

CrossRef Full Text | Google Scholar

3. Lee JH, Shklovskaya E, Lim SY, Carlino MS, Menzies AM, Stewart A, et al. Transcriptional downregulation of MHC class I and melanoma de- differentiation in resistance to PD-1 inhibition. Nat Commun (2020) 11:1897. doi: 10.1038/s41467-020-15726-7

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Olson JA, Leveson-Gower DB, Gill S, Baker J, Beilhack A, Negrin RS. NK cells mediate reduction of GVHD by inhibiting activated, alloreactive T cells while retaining GVT effects. Blood (2010) 115:4293–301. doi: 10.1182/blood-2009-05-222190

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lundqvist A, McCoy JP, Samsel L, Childs R. Reduction of GVHD and enhanced antitumor effects after adoptive infusion of alloreactive Ly49-mismatched NK cells from MHC-matched donors. Blood (2007) 109:3603–6. doi: 10.1182/blood-2006-05-024315

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Shembrey C, ND; H, Hollande F. Impact of tumor and immunological heterogeneity on the anti-cancer immune response. Cancers (Basel) (2019) 11:1217. doi: 10.3390/cancers11091217

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu E, Tong Y, Dotti G, Shaim H, Savoldo B, Mukherjee M, et al. Cord blood NK cells engineered to express IL-15 and a CD19-targeted CAR show long-term persistence and potent antitumor activity. Leukemia (2018) 32:520–31. doi: 10.1038/leu.2017.226

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Han J, Chu J, Keung Chan W, Zhang J, Wang Y, Cohen JB, et al. CAR-engineered NK cells targeting wild-type EGFR and EGFRvIII enhance killing of glioblastoma and patient-derived glioblastoma stem cells. Sci Rep (2015) 5:11483. doi: 10.1038/srep11483

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Bagot M, Porcu P, Marie-Cardine A, Battistella M, William BM, Vermeer M, et al. IPH4102, a first-in-class anti-KIR3DL2 monoclonal antibody, in patients with relapsed or refractory cutaneous T-cell lymphoma: An international, first-in-human, open-label, phase 1 trial. Lancet Oncol (2019) 20:1160–70. doi: 10.1016/S1470-2045(19)30320-1

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Vey N, Dumas P, Recher C, Gastaud L, Lioure B, Bulabois C, et al. Randomized phase 2 trial of lirilumab (anti-KIR monoclonal antibody, mAb) as maintenance treatment in elderly patients (pts) with acute myeloid leukemia (AML): Results of the effikir trial. Blood (2017) 130:889. doi: 10.1182/blood.V130.Suppl_1.889.889

CrossRef Full Text | Google Scholar

11. Cohen R, Fayette J, Posner M, Surgery G--. Phase II study of monalizumab, a first-in-class NKG2A monoclonal antibody, in combination with cetuximab in previously treated recurrent or metastatic. Innate-PharmaCom (2018). doi: 10.1158/1538-7445.AM2018-CT158

CrossRef Full Text | Google Scholar

12. Imai K, Matsuyama S, Miyake S, Suga K, Nakachi K. Natural cytotoxic activity of peripheral-blood lymphocytes and cancer incidence: An 11-year follow-up study of a general population. Lancet (2000) 356:1795–9. doi: 10.1016/S0140-6736(00)03231-1

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Hanna N, Burton RC. Definitive evidence that natural killer (NK) cells inhibit experimental tumor metastases in vivo. J Immunol (1981) 127:1754–8.

PubMed Abstract | Google Scholar

14. Sathe P, Delconte RB, Souza-Fonseca-Guimaraes F, Seillet C, Chopin M, Vandenberg CJ, et al. Innate immunodeficiency following genetic ablation of Mcl1 in natural killer cells. Nat Commun (2014) 5:4539. doi: 10.1038/ncomms5539

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Coca S, Perez-Piqueras J, Martinez D, Colmenarejo A, Saez MA, Vallejo C, et al. The prognostic significance of intratumoral natural killer cells in patients with colorectal carcinoma. Cancer (1997) 79:2320–8. doi: 10.1002/(SICI)1097-0142(19970615)79:12<2320::AID-CNCR5>3.0.CO;2-P

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kondo E, Koda K, Takiguchi N, Oda K, Seike K, Ishizuka M, et al. Preoperative natural killer cell activity as a prognostic factor for distant metastasis following surgery for colon cancer. Digestive Surgery. (2003) 20:445–51. doi: 10.1159/000072714

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Donadon M, Hudspeth K, Cimino M, Di Tommaso L, Preti M, Tentorio P, et al. Increased infiltration of natural killer and T cells in colorectal liver metastases improves patient overall survival. J Gastrointestinal Surg (2017) 21:1–11. doi: 10.1007/s11605-017-3446-6

CrossRef Full Text | Google Scholar

18. Kärre K, Ljunggren HG, Piontek G, Kiessling R. Selective rejection of h-2-deficient lymphoma variants suggests alternative immune defence strategy. Nature (1986) 319:675–8. doi: 10.1038/319675a0

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Karlhofer FM, Ribaudo RK, Yokoyama WM. MHC class I alloantigen specificity of ly-49+ IL-2-activated natural killer cells. Nature (1992) 358:66–70. doi: 10.1038/358066a0

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Dijkstra KK, Cattaneo CM, Weeber F, Chalabi M, van de Haar J, Fanchi LF, et al. Generation of tumor-reactive T cells by Co-culture of peripheral blood lymphocytes and tumor organoids. Cell (2018) 174:1586–1598.e12. doi: 10.1016/j.cell.2018.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Dierssen JWF, de Miranda NFCC, Ferrone S, van Puijenbroek M, Cornelisse CJ, Fleuren GJ, et al. HNPCC versus sporadic microsatellite-unstable colon cancers follow different routes toward loss of HLA class I expression. BMC Cancer (2007) 7:33. doi: 10.1186/1471-2407-7-33

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Halama N, Braun M, Kahlert C, Spille A, Quack C, Rahbari N, et al. Natural killer cells are scarce in colorectal carcinoma tissue despite high levels of chemokines and cytokines. Clin Cancer Res (2011) 17:678–89. doi: 10.1158/1078-0432.CCR-10-2173

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Tan IB, Tan P. Genetics: An 18-gene signature (ColoPrint®) for colon cancer prognosis. Nat Rev Clin Oncol (2011) 8:131–3. doi: 10.1038/nrclinonc.2010.229

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Chen H, Sun X, Ge W, Qian Y, Bai R, Zheng S. A seven-gene signature predicts overall survival of patients with colorectal cancer. Oncotarget (2017) 8:95054–65. doi: 10.18632/oncotarget.10982

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gao S, Tibiche C, Zou J, Zaman N, Trifiro M, O’Connor-McCourt M, et al. Identification and construction of combinatory cancer hallmark-based gene signature sets to predict recurrence and chemotherapy benefit in stage II colorectal cancer. JAMA Oncol (2016) 2:37–45. doi: 10.1001/jamaoncol.2015.3413

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Tong M, Zheng W, Li H, Li X, Ao L, Shen Y, et al. Multi-omics landscapes of colorectal cancer subtypes discriminated by an individualized prognostic signature for 5-fluorouracil-based chemotherapy. Oncogenesis (2016) 5:e242–2. doi: 10.1038/oncsis.2016.51

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Cao B, Luo L, Feng L, Ma S, Chen T, Ren Y, et al. A network-based predictive gene-expression signature for adjuvant chemotherapy benefit in stage II colorectal cancer. BMC Cancer (2017) 17:844. doi: 10.1186/s12885-017-3821-4

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Watanabe T, Kobunai T, Yamamoto Y, Matsuda K, Ishihara S, Nozawa K, et al. Gene expression signature and response to the use of leucovorin, fluorouracil and oxaliplatin in colorectal cancer patients. Clin Trans Oncol (2011) 13:419–25. doi: 10.1007/s12094-011-0676-z

CrossRef Full Text | Google Scholar

29. Costanzo MC, Kim D, Creegan M, Lal KG, Ake JA, Currier JR, et al. Transcriptomic signatures of NK cells suggest impaired responsiveness in HIV-1 infection and increased activity post-vaccination. Nat Commun (2018) 9:1–16. doi: 10.1038/s41467-018-03618-w

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Cursons J, Souza-Fonseca Guimaraes F, Foroutan M, Anderson A, Hollande F, Hediyeh-Zadeh S, et al. A gene signature predicting natural killer cell infiltration and improved survival in melanoma patients. Cancer Immunol Res (2019) 7(7):1162–74. doi: 10.1158/2326-6066.CIR-18-0500

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Eckl J, Buchner A, Prinz PU, Riesenberg R, Siegert SI, Kammerer R, et al. Transcript signature predicts tissue NK cell content and defines renal cell carcinoma subgroups independent of TNM staging. J Mol Med (2012) 90:55–66. doi: 10.1007/s00109-011-0806-7

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Foroutan M, Molania R, Pfefferle A, Behrenbruch C, Scheer S, Kallies A, et al. The ratio of exhausted to resident infiltrating lymphocytes is prognostic for colorectal cancer patient outcome. Cancer Immunol Res (2021) 9:1125–40. doi: 10.1158/2326-6066.CIR-21-0137

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wang J, Mouradov D, Wang X, Jorissen RN, Chambers MC, Zimmerman LJ, et al. Colorectal cancer cell line proteomes are representative of primary tumors and predict drug sensitivity. Gastroenterology (2017) 153:1082–95. doi: 10.1053/j.gastro.2017.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Linsley PS, Speake C, Whalen E, Chaussabel D. Copy number loss of the interferon gene cluster in melanomas is linked to reduced T cell infiltrate and poor patient prognosis. PloS One (2014) 9(10):e109760. doi: 10.1371/journal.pone.0109760

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Monaco G, Lee B, Xu W, Mustafah S, Hwang YY, Carré C, et al. RNA-Seq signatures normalized by mRNA abundance allow absolute deconvolution of human immune cell types. Cell Rep (2019) 26:1627–1640.e7. doi: 10.1016/j.celrep.2019.01.041

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O’Brien SA, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell (2020) 181:442–59.e29. doi: 10.1016/j.cell.2020.03.048

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin AA, Kim S, et al. The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature (2012) 483:603–7. doi: 10.1038/nature11003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Calon A, Espinet E, Palomo-Ponce S, Tauriello DVF, Iglesias M, Céspedes MV, et al. Dependency of colorectal cancer on a TGF-β-Driven program in stromal cells for metastasis initiation. Cancer Cell (2012) 22:571–84. doi: 10.1016/j.ccr.2012.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Dogra P, Rancan C, Ma W, Toth M, Senda T, Carpenter DJ, et al. Tissue determinants of human NK cell development, function, and residence. Cell (2020) 180:749–63.e13. doi: 10.1016/j.cell.2020.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Fehniger TA, Cai SF, Cao X, Bredemeyer AJ, Presti RM, French ARR, et al. Acquisition of murine NK cell cytotoxicity requires the translation of a pre-existing pool of granzyme b and perforin mRNAs. Immunity (2007) 26:798–811. doi: 10.1016/j.immuni.2007.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Pfefferle A, Netskar H, Ask EH, Lorenz S, Goodridge JP, Sohlberg E, et al. A temporal transcriptional map of human natural killer cell differentiation. bioRxiv (2019), 630657. doi: 10.1101/630657

CrossRef Full Text | Google Scholar

43. Lima M, Leander M, Santos M, Santos AH, Lau C, Queirós ML, et al. Chemokine receptor expression on normal blood CD56+ NK-cells elucidates cell partners that comigrate during the innate and adaptive immune responses and identifies a transitional NK-cell population. J Immunol Res (2015) 2015:839684. doi: 10.1155/2015/839684

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Yang C, Siebert JR, Burns R, Gerbec ZJ, Bonacci B, Rymaszewski A, et al. Heterogeneity of human bone marrow and blood natural killer cells defined by single-cell transcriptome. Nat Commun (2019) 10(1):3931. doi: 10.1038/s41467-019-11947-7

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Smith SL, Kennedy PR, Stacey KB, Worboys JD, Yarwood A, Seo S, et al. Diversity of peripheral blood human NK cells identified by single cell RNA sequencing. Blood Adv (2019) 4:1388–406. doi: 10.1182/bloodadvances.2019000699

CrossRef Full Text | Google Scholar

46. Yu J, Mao HC, Wei M, Hughes T, Zhang J, Park IK, et al. CD94 surface density identifies a functional intermediary between the CD56bright and CD56dim human NK-cell subsets. Blood (2010) 115:274–81. doi: 10.1182/blood-2009-04-215491

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Abbas AR, Baldwin D, Ma Y, Ouyang W, Gurney A, Martin F, et al. Immune response in silico (IRIS): Immune-specific genes identified from a compendium of microarray expression data. Genes Immunity (2005) 6:319–31. doi: 10.1038/sj.gene.6364173

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Crinier A, Milpied P, Escalière B, Piperoglou C, Galluso J, Balsamo A, et al. High-dimensional single-cell analysis identifies organ-specific signatures and conserved NK cell subsets in humans and mice. Immunity (2018) 49:971–986.e5. doi: 10.1016/j.immuni.2018.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Bezman NA, Kim CC, Sun JC, Min-Oo G, Hendricks DW, Kamimura Y, et al. Molecular definition of the identity and activation of natural killer cells. Nat Immunol (2012) 13:1000–8. doi: 10.1038/ni.2395

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinf (2018) 19:1–10. doi: 10.1186/s12859-018-2435-4

CrossRef Full Text | Google Scholar

51. Marisa L, de Reyniès A, Duval A, Selves J, Gaub MP, Vescovo L, et al. Gene expression classification of colon cancer into molecular subtypes: Characterization, validation, and prognostic value. PloS Med (2013) 10:e1001453. doi: 10.1371/journal.pmed.1001453

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Hu X, Zhang L, Li Y, Ma X, Dai W, Gao X, et al. Organoid modelling identifies that DACH1 functions as a tumour promoter in colorectal cancer by modulating BMP signalling. EBioMedicine (2020) 56:102800. doi: 10.1016/j.ebiom.2020.102800

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wang X, Chen J, Wang J, Yu F, Zhao S, Zhang Y, et al. Metalloproteases meprin-α (MEP1A) is a prognostic biomarker and promotes proliferation and invasion of colorectal cancer. BMC Cancer (2016) 16:383. doi: 10.1186/s12885-016-2460-5

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zeng J, Yang X, Cheng L, Liu R, Lei Y, Dong D, et al. Chemokine CXCL14 is associated with prognosis in patients with colorectal carcinoma after curative resection. J Trans Med (2013) 11:6. doi: 10.1186/1479-5876-11-6

CrossRef Full Text | Google Scholar

55. Westrich JA, Vermeer DW, Silva A, Bonney S, Berger JN, Cicchini L, et al. CXCL14 suppresses human papillomavirus-associated head and neck cancer through antigen-specific CD8+ T-cell responses by upregulating MHC-I expression. Oncogene (2019) 38:7166–80. doi: 10.1038/s41388-019-0911-6

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Bucks CM, Norton JA, Boesteanu AC, Mueller YM, Katsikis PD. Chronic antigen stimulation alone is sufficient to drive CD8 + T cell exhaustion. J Immunol (2009) 182:6697–708. doi: 10.4049/jimmunol.0800997

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Nersesian S, Schwartz SL, Grantham SR, MacLean LK, Lee SN, Pugh-Toole M, et al. NK cell infiltration is associated with improved overall survival in solid cancers: A systematic review and meta-analysis. Trans Oncol (2021) 14(1):100930. doi: 10.1016/j.tranon.2020.100930

CrossRef Full Text | Google Scholar

58. Malaise M, Rovira J, Renner P, Eggenhofer E, Sabet-Baktach M, Lantow M, et al. KLRG1+ NK cells protect T-bet-Deficient mice from pulmonary metastatic colorectal carcinoma. J Immunol (2014) 192:1954–61. doi: 10.4049/jimmunol.1300876

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Takeda K, Hayakawa Y, Smyth MJ, Kayagaki N, Yamaguchi N, Kakuta S, et al. Involvement of tumor necrosis factor-related apoptosis-inducing ligand in surveillance of tumor metastasis by liver natural killer cells. Nat Med (2001) 7:94–100. doi: 10.1038/83416

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Platonova S, Cherfils-Vicini J, Damotte D, Crozet L, Vieillard V, Validire P, et al. Profound coordinated alterations of intratumoral NK cell phenotype and function in lung carcinoma. Cancer Res (2011) 71:5412–22. doi: 10.1158/0008-5472.CAN-10-4179

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Gillard-Bocquet M, Caer C, Cagnard N, Crozet L, Perez M, Fridman WH, et al. Lung tumor microenvironment induces specific gene expression signature in intratumoral NK cells. Front Immunol (2013) 4:19. doi: 10.3389/fimmu.2013.00019

PubMed Abstract | CrossRef Full Text | Google Scholar

62. De Andrade LF, Lu Y, Luoma A, Ito Y, Pan D, Pyrdol JW, et al. Discovery of specialized NK cell populations infiltrating human melanoma metastases. JCI Insight (2019) 4(23):e133103. doi: 10.1172/jci.insight.133103

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Pizzolato G, Kaminski H, Tosolini M, Franchini DM, Pont F, Martins F, et al. Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVδ1 and TCRVδ2 γδ T lymphocytes. Proc Natl Acad Sci USA (2019) 116:11906–15. doi: 10.1073/pnas.1818488116

PubMed Abstract | CrossRef Full Text | Google Scholar

64. D’Ombrain MC, Hansen DS, Simpson KM, Schofield L. γδ-T cells expressing NK receptors predominate over NK cells and conventional T cells in the innate IFN-γ response to plasmodium falciparum malaria. Eur J Immunol (2007) 37:1864–73. doi: 10.1002/eji.200636889

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Von Lilienfeld-Toal M, Nattermann J, Feldmann G, Sievers E, Frank S, Strehl J, et al. Activated γδ T cells express the natural cytotoxicity receptor natural killer p44 and show cytotoxic activity against myeloma cells. Clin Exp Immunol (2006) 144:528–33. doi: 10.1111/j.1365-2249.2006.03078.x

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Correia DV, Fogli M, Hudspeth K, Gomes Da Silva M, Mavilio D, Silva-Santos B. Differentiation of human peripheral blood Vδ1+ T cells expressing the natural cytotoxicity receptor NKp30 for recognition of lymphoid leukemia cells. Blood (2011) 118:992–1001. doi: 10.1182/blood-2011-02-339135

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Miyazaki K, Yamaguchi M, Imai H, Kobayashi T, Tamaru S, Nishii K, et al. Gene expression profiling of peripheral T-cell lymphoma including {gamma}{delta} T-cell lymphoma. Blood (2009) 113:1071–4. doi: 10.1182/blood-2008-07-166363

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Angelova M, Charoentong P, Hackl H, Fischer ML, Snajder R, Krogsdam AM, et al. Characterization of the immunophenotypes and antigenomes of colorectal cancers reveals distinct tumor escape mechanisms and novel targets for immunotherapy. Genome Biol (2015) 16:64. doi: 10.1186/s13059-015-0620-6

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Pagès F, Mlecnik B, Marliot F, Bindea G, Ou FS, Bifulco C, et al. International validation of the consensus immunoscore for the classification of colon cancer: A prognostic and accuracy study. Lancet (2018) 391:2128–39. doi: 10.1016/S0140-6736(18)30789-X

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Guinney J, Dienstmann R, Wang X, De Reyniès A, Schlicker A, Soneson C, et al. The consensus molecular subtypes of colorectal cancer. Nat Med (2015) 21:1350–6. doi: 10.1038/nm.3967

PubMed Abstract | CrossRef Full Text | Google Scholar