γδ T Cells in Merkel Cell Carcinomas Have a Proinflammatory Profile Prognostic of Patient Survival

Merkel cell carcinoma (MCC) is a rare neuroendocrine skin cancer caused by either oncoviral infection or UV carcinogenesis. In viral-positive MCC (VP-MCC), integration of the Merkel cell polyoma virus (MCPyV) genome into the host cell causes cellular transformation via the early transforming genes, which encode the small and large T-antigens (1, 2). In contrast, viral-negative MCC (VN-MCC) is associated with sun-induced somatic hypermutation (3–5). For both subtypes of MCC, immune evasion is likely important for tumor initiation and progression. Immunosuppression is a known risk factor for MCC, and age-related immunosenescence is also likely responsible for increased incidence primarily in the elderly. Spontaneous tumor regression and MCC of unknown primary origin are also well documented, suggesting immune clearance of the primary tumor has occurred (6, 7). Importantly, almost half of patients with advanced MCC respond to anti–PD-1/PD-L1 immune checkpoint inhibitor (ICI) therapies with durable objective responses (8–10).

There is strong evidence to support the role of host-adaptive immune responses in MCC. Heavy T-cell infiltrates are observed in approximately 20% of cases and have prognostic and predictive significance (11–13). T cells reactive to MCPyV peptides have been isolated from VP-MCC patients' blood and tumors, showing likely tumor specificity (14–17), while VN-MCC also exhibit T-cell responses, presumably due to recognition of neoantigens (5). A loss of antigen-presenting machinery in tumor cells, including the downregulation of the MHC class I complex, which is frequent in MCC (18), may facilitate immune escape and has been implicated in cell-based immunotherapy resistance (19).

Conventional MHC-restricted αβ T cells are important for cancer immune surveillance; however, there is also growing interest in “unconventional” T cells that lack MHC-restriction (20). γδ T cells are a major subset of such cells and are defined by expression of a T-cell receptor (TCR) encoded by the TCR-γ and TCR-δ genes (20). During thymic development, expression of a γδ TCR drives an effector phenotype distinct from the αβ T–cell lineage. γδ T cells circulate in a preprimed effector state, poised to elicit a rapid effector response (20). Upon activation, γδ T cells rapidly proliferate and are proinflammatory, producing cytokines such as IFNγ and IL17, and releasing cytotoxic granules of granzymes and perforin. γδ T cells are emerging as key players in tumor immune surveillance, and have been identified in tumors of the skin, colon, breast, and pancreas (21–27). However, their role in MCC is currently unknown.

Here, we studied tumor-infiltrating T cells in MCC using multiplex IHC/immunofluorescence (m-IHC/IF), flow cytometry, and single-cell RNA sequencing (scRNA-seq). We observed abundant CD4/CD8 double-negative (DN) T cells in some tumors, indicating the presence of unconventional T cells. These cells were characterized as γδ T cells using flow cytometry, and scRNA-seq revealed a proinflammatory gene-expression signature that was prognostic in an expanded cohort of patients.

Patient cohorts were collected from three centers according to protocols approved by their respective institutional review boards (IRB) according to the Helsinki Declaration. These centers were the Peter MacCallum Cancer Centre (Protocols: 14–113 and 10–16; n = 51), the Sydney Local Health District Ethics Committee [Protocol No. X15–0454 (prev. X11–0289) and HREC/11/RPAH/444 and Protocol No. X17–0312 (prev. X11–0023) and HREC/11/RPAH/32 (n = 12)], and the Northern Sydney Local Health District (Protocol HREC1312–417M; n = 6). Either written patient consent was obtained or a waiver of patient consent was granted by respective IRBs for access to archive samples (see Supplementary Tables S1–S3 for patient and sample details).

Formalin-fixed, paraffin-embedded (FFPE) sections (4 μm) were baked for 1 hour at 60°C on SuperFrost Plus slides prior to dewaxing using histolene followed by graded ethanol washes and then antigen retrieval (AR). Sections were stained for 30 minutes in primary antibody, blocked in 0.3% H₂O₂ for 10 minutes, incubated for 10 minutes in the appropriate secondary antibody and then 10 minutes in TSA reagent (OPAL 7-color IHC kit, PerkinElmer, catalog no. NEL811001KT). After staining, microwave heat removal of the antibody conjugates was performed in AR buffer, allowing for serial staining. Primary antibodies, tyramide signal amplification (TSA) dyes, order of staining and AR conditions are detailed in Supplementary Table S4. Sections were counterstained with DAPI. Single antibody control slides were stained as imaging controls. m-IHC/IF stained sections were imaged on the Vectra 3.0 Automated Quantitative Pathology Imaging System (PerkinElmer). InForm Software v2.2 (PerkinElmer) was used to spectrally unmix images, perform segmentation of tissue regions and individual cells, and cell phenotyping.

FFPE (4-μm) tissue sections baked on SuperFrost Plus slides were stained for MCPyV (Santa Cruz Biotechnology, catalog no. sc-136172), CD3 (Abcam, catalog no. ab135372) and TCRD (Santa Cruz Biotechnology, catalog no. sc-100289) using traditional 3,3′-Diaminobenzidine (DAB) IHC methods. Details of antibodies used, conditions applied, and the staining platform can be found in Supplementary Table S5.

Representative tumor pieces (∼1 cm³) were used for dissociation. Surrounding nontumor tissue was macroscopically removed. Tumors were finely diced into RPMI1640 (Thermo Fisher Scientific, catalog no. A1049101) containing 20% FBS (Thermo Fisher Scientific, catalog no. 10099133), 1 mg/mL collagenase type 4 (Worthington Biochemical), 30 U/mL DNase (Sigma-Aldrich), 10 μg/mL gentamicin sulfate (Pfizer Australia) and incubated for 30 minutes at 37°C on a rocker. The digested tumor was teased through a 100-μm sieve. A subset of tumors from the Melanoma Institute of Australia were dissociated into single-cell suspensions using the human Tumor Dissociation Kit (Miltenyi Biotec, catalog no. 130–095–929) and gentleMACS Dissociator (Miltenyi Biotec, catalog no. 130–093–235), according to the manufacturer's instructions. Pelleted cells were resuspended in RPMI1640 containing 10% FBS and used for FACS analysis or resuspended in 10% DMSO/90% FBS and frozen at −80°C.

Thawed samples were incubated in DNAse I (Sigma-Aldrich, catalog no. 4716728001) for 20 minutes at 37°C. Cells were filtered, washed in PBS and incubated in PBS + 50 nmol/L dasatinib (StemCell Technologies, catalog no. 73082) + DNAse at 37°C for 1 hour. Cells were stained with viability dye LIVE/DEAD Fixable Near Infrared (Thermo Fisher Scientific, catalog no. L34975) followed by blocking with Fc block (BD Biosciences, catalog no. 564219; 500 ng/test) and 5% normal mouse serum (generated in-house). Cells were then stained in FACS buffer (2% FBS in PBS) with surface mAbs and tetramers (see Supplementary Table S6 for the reagents used and the dilutions of the reagents). Intracellular staining was performed using the BD Cytofix/Cytoperm Kit (BD Biosciences, catalog no. 554717). Remaining samples were fixed in 2% paraformaldehyde (PFA), washed, filtered, and immediately analyzed using a BD LSR Fortessa. Analysis of case P007 was performed independently where cells were labeled with mAbs in FACS buffer for 30 minutes at 4°C, washed twice in PBS, fixed in 2% PFA, washed, and analyzed on a BD LSRII.

Cells were stained as above. After sorting into sort buffer (PBS + 1% BSA and 0.5 U/μL of RNAse inhibitor, Thermo Fisher Scientific, catalog no. AM2694) using a BD FACSAria III, cells were counted, cell density adjusted, and immediately used in single-cell sequencing experiments.

Data were analyzed using Flowjo software (Treestar Inc.). For sorting and analysis of the 11-sample cohort, T cells were defined as LIVE/DEAD⁻, CD19⁻, CD14⁻, CD45⁺ and CD3⁺ after doublet exclusion. γδ T cells were also TCRγδ⁺. Example gating is shown in Supplementary Fig. S1.

cDNA from single T cells sorted into 96-well PCR plates was generated using Superscript VILO (Thermo Fisher Scientific, catalog no. 11754050) and 0.1% Triton X-100. TCR-γ and -δ transcripts were amplified by two rounds of semi-nested multiplexed PCR as described previously (28). TRGC or TRDC internal reverse primers were added to PCR products and then sent for Sanger sequencing at the Australian Genome Research Facility (ARGF; primer sequences were as in ref. 28). Results were analyzed using the IMGT database.

FACS-sorted cells were loaded on a Chromium Single cell Controller (10x Genomics) using the Single-Cell Immune-Profiling Solution (Chromium Single Cell 5′ Library & Gel Bead Kit, catalog no. PN-1000002), following the manufacturer's recommendations for GEM Generation and Barcoding, Post GEM-RT clean up, cDNA Amplification and quality control. TCR enrichment was carried out using the Chromium Single Cell V(D)J Enrichment Kit (catalog no. PN-1000005), human T cell (αβ transcripts), or custom-designed primers for γδ TCR transcripts enrichment (Supplementary Fig. S2). Gene expression and TCR libraries (αβ and γδ) were prepared following manufacturer's recommendations. All libraries were sequenced using 150 bp paired-end sequencing on the Illumina NovaSeq 6000. Five thousand and 50,000 read-pairs/cell were generated for TCR and gene expression libraries, respectively (see Supplementary Table S7 for scRNA-seq run summary statistics and Supplementary Data for gene expression and TCR-seq annotations, respectively).

RNA extraction was done using the RNeasy FFPE RNA extraction kit according to manufacturer's protocols (Qiagen, catalog no. 35047). RNA-seq library preparation and hybridization capture enrichment were done using either Illumina RNA-Access (Illumina, catalog no. 15049525) or Agilent RNA-Direct (Agilent, catalog no. G7564A) kits according to manufacturers' protocols. Total RNA-seq libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (rRNA depletion workflow; NEB, catalog no. E7760L) as previously described (29). Sequencing was done using Illumina HiSeq 2000, NextSeq and NovaSeq 6000 (Illumina; see Supplementary Table S8 and Supplementary Data for bulk RNA-seq statistics and processed count data, respectively). Raw read-level sequencing data cannot be made publicly available owing to restrictions on patient consent preventing sharing of potentially identifiable genetic data.

scRNA-seq binary base call (BCL) files were demultiplexed and converted into FASTQ files using BCLtoFastq prior to alignment to hg38 with CellRanger (v3.1.0) and analyzed in R studio with Seurat (v3.1.1). Cells with a mitochondrial proportion greater than 10% and/or a feature count <800 or >4,000 were discarded. Doublet finding was performed with Scrublet (v0.2.1). Sample batches were combined, normalized and scaled with SCTransform. Cell cycle–associated genes were regressed from the data using the Seurat workflow. Uniform Manifold Approximation and Projection (UMAP) clustering was performed selecting the first nine principal components. scTCR-seq αβ and γδ TCR FASTQs were aligned with CellRanger (v3.0.2) to the default CellRanger VDJ reference and a custom γδ reference constructed from the IMGT database, respectively. Output contig annotations were filtered and analyzed in R studio with Seurat.

RNA-Seq FASTQ files were aligned to hg19 and raw counts were quantified with featurecount via bcbio-nextgen (v1.1.6a-b). Count data were filtered for lowly expressed genes, TMM-normalized and logCPM converted with limma (v3.42.0) and edgeR (v3.28.1). A previously published MCC microarray dataset was selected because of the availability of clinical data kindly supplied by the authors of the original study (Supplementary Table S1; ref. 11). The microarray data was sourced from the NCBI Gene Expression Omnibus (GEO) database (GSE22396; ref. 11) and read into R with GEOquery package (v2.52.0), filtered, and RMA-normalized with affy (v1.62.0). MCC cell line (GSE124857; ref. 30) and innate immune cell and T cell (GSE124731; ref. 31) RNA-seq datasets were sourced from NCBI GEO database.

Gene-set variance analysis (GSVA) using RNA-seq or microarray data was performed using GSVA (v1.34.0). Prior to GSVA, RNA-seq datasets were filtered and converted to log₂-CPM, and the microarray data was filtered and RMA-normalized. RNA-seq cohorts were integrated and batch corrected with limma (v3.42.0). Signature genes in the RNA-seq set were scaled by subtracting the mean gene logCPM values and microarray expression values were converted to z-scores for visualization of gene expression compared with GSVA scores.

Analysis of the cell phenotype populations and comparison across tumors was performed in GraphPad PRISM and R. Survival analysis was performed using survival package (v3.1–8). For MCC-specific overall survival, non-MCC–related events were censored. Statistical significance was calculated with a log-rank test. Multiple Cox regression analysis was performed, and Kaplan–Meier and forest plots were produced with the survminer package (v0.4.8).

The Jurkat-76.MG2 cell line has previously been described (32). For MG1, MG3, and MG4, matched TCR-γ and -δ genes were cloned into the pMIG expression vector as per MG2 (32). Vectors were then transduced, along with a pMIG vector encoding full-length CD3ϵ, δ, γ and ζ subunits, into β 2m-deficient, αβTCR-deficient jurkat-76 cells as previously described (33). C1R.MR1 (34) C1R.CD1a (35), C1R.CD1b(36), C1R.CD1c (37), and C1R.CD1d (38) lines are previously described. Jurkat-76 cells provided by Zhenjun Chen (University of Melbourne) who received them in 2015 from Mirjam H.M. Heemskerk (ref. 39; Leiden University Medical Center, Leiden, the Netherlands). The cell line has not been authenticated. The cell line has returned negative Mycoplasma tests. Cells were used for up to approximately 10 passages. Cells were maintained in RF10 complete media consisting of RPMI-1640 base (Thermo Fisher Scientific, catalog no. A1049101) supplemented with 10% FBS (Thermo Fisher Scientific, catalog no. 10099133), penicillin (100 U/mL), streptomycin (100 μg/mL), Glutamax (2 mmol/L), sodium pyruvate (1 mmol/L), nonessential amino acids (0.1 mmol/L), HEPES buffer (15 mmol/L), pH 7.2 to 7.5 (all from Thermo Fisher Scientific, catalog numbers: 15140163, 35050061, 11360070, 11140050, and 15630080, respectively), and 2-mercaptoethanol (50 μmol/L, Sigma, catalog no. M-7522).

5×10⁴ Jurkat-76.TCR cells were cocultured with 5 × 10⁴ cell-trace violet (CTV; Thermo Fisher Scientific, catalog no. C34557)–labeled C1R cells for approximately 18 hours. To label C1R cells with CTV, cells were washed twice with PBS before being resuspended in 500 μL PBS + 1 μmol/L CTV and incubated at 37°C for 15 minutes. 500 μL RF10 complete media (as above) was then added before cells were pelleted and resuspended in RF10 ready for counting. For blocking experiments, anti-MR1 (clone 26.5) or anti-CD1c (clone L161; BioLegend, catalog no. 361102 and catalog no. 331502, respectively) were included at 10 μg/mL. After culture, cells were stained with LIVE/DEAD NIR and anti-human CD69, fixed in 2% PFA and immediately acquired on an LSRFortessa (BD Biosciences). Jurkat cells were gated on the basis of FSC-A and SSC-A, and defined as CTV⁻, GFP⁺ cells after dead cell and doublet removal.

To study T cells in MCC, we performed quantitative analysis using mIHC/IF. FFPE tissue sections were stained with antibodies targeting T-cell epitopes (CD3, CD4, CD8, FOXP3), PD-L1, synaptophysin (tumor), and DAPI to detect cell nuclei (Fig. 1A). Between 2 and 35 high power (20× magnification) regions of interest (ROI = 1.338 mm²) were selected depending on tumor area encompassing the tumor–stroma interface and central tumor regions. Tumor and stromal compartments were segmented on the basis of synaptophysin staining and the T cells phenotyped on the basis of their respective markers (cytotoxic T cells were defined as CD3⁺CD8⁺; helper T cells were defined as CD3⁺CD4⁺; regulatory T cells were defined as CD3⁺FOXP3⁺; DN T cells were defined as CD3⁺CD4⁻CD8⁻). T-cell populations were enumerated within intratumoral regions or within complete ROIs (including the stroma), and cell counts were averaged across ROIs within individual tissue sections.

We used m-IHC/IF to analyze samples from a cohort of 39 patients treated at three institutions. Of the 39 patients, 28 had VN-MCC and 11 VP-MCC by MCPyV large T-antigen IHC staining. A total of 58 FFPE tumor samples were stained, including 34 primary and 24 metastatic lesions. Matched primary and metastatic lesions were available for 17 cases (Supplementary Table S2). Primary VP-MCC had significantly higher intratumoral CD8⁺ (P < 0.001) and CD4⁺ (P = 0.0456) cell infiltrates compared with primary VN-MCC. Among metastatic lesions, only CD4⁺ T cells were significantly higher in VP-MCC tumors in total ROIs (P = 0.03; Fig. 1B). No significant differences in PD-L1 staining were observed between VP-MCC and VN-MCC (Supplementary Fig. S3).

Tumors were classified into four groups based on the dominant CD4/CD8 staining profile (Fig. 1C). Among 34 primary tumors, the single largest group had a predominance of CD4⁺ T cells (50%) followed by CD8⁺ and DN T cells (with 12% each), while 26% had no dominant subset. Among the matched primary–metastatic pairs, the dominant T-cell subset was consistent in 10 of 17 (59%) cases (Supplementary Fig. S4). In 40% of MCC at least 20% of all infiltrating T cells were DN T cells, and a subset of MCC showed a dominant DN T-cell infiltration pattern.

Flow cytometry was used to characterize immune cells in 11 MCC tumors (3 primary and 8 metastatic) from 10 patients (Supplementary Table S3). Consistent with the m-IHC/IF data, CD4⁺ and CD8⁺ cells were the most abundant T cells, ranging from 14% to 82% and 16% to 77%, with medians of 57.1% and 39.4%, respectively. DN T cells ranged from 1% to 25% of CD3⁺ cells, with a median of 4.0% (Fig. 2A). DN T cells are a heterogeneous collection of unconventional T cells including CD1d-restricted Type I natural killer T (NKT) cells, mucosal-associated invariant T (MAIT) cells, and γδ T cells. Analysis of MCC samples revealed no appreciable NKT cells and very few MAIT cells (median 0.3%) as determined by CD1d–α-GalCer and MR1–5-OP-RU tetramer staining. In one tumor, MAIT cells represented 16.9% of T cells; however, this was a liver metastasis, and the liver is known to be enriched for MAIT cells (Fig. 2A; Supplementary Fig. S5; ref. 39). In contrast, γδ T cells accounted for a median of 2.8% of total T cells, ranging from 0.3% to 27.5%.

γδ T cells can be broadly subdivided on the basis of TCR-δ usage, with Vδ2⁺ and Vδ2⁻ γδ T cells representing functionally distinct populations. Vδ2 mAb staining revealed that most γδ T cells were Vδ2⁻, with Vδ2⁻ and Vδ2⁺ cells accounting for medians of 1.3% and 0.2% of total T cells, with ranges of 0.1% to 27.3% and 0% to 2.4%, respectively (Fig. 2A). Overall, γδ T cells accounted for a median of 73.1% of tumor-infiltrating DN T cells.

In one tumor (P055), >25% of total infiltrating T cells were DN T cells, and these were almost exclusively Vδ2⁻ γδ T cells (Fig. 2B). Phenotypic analysis of total Vδ2⁻ γδ T cells in this tumor demonstrated that although 90% were DN, approximately 10% were CD8⁺ (Fig. 2C). Moreover, the DN cells displayed an enrichment of terminally differentiated cells (CD45RA⁺, CD27⁻) relative to CD4⁺ and CD8⁺ T cells in the same tumor, and they expressed markers of activation and tissue retention such as CD69 and CD103, including co-expression of these markers by 27% of cells, which is indicative of a tissue-resident memory (TRM) phenotype (Fig. 2D). Similar to their CD8⁺ T-cell counterparts, these cells had high expression of immune checkpoint markers PD-1 and LAG3, suggesting an exhaustion phenotype (Fig. 2D).

To gauge the clonality of γδ T cells from P055, we performed TCR analysis on sorted single cells. Sequencing of 59 cells revealed an oligoclonal repertoire dominated by four clonotypes. The Vδ1/Vγ4 TCR (clone MG1), Vδ1/Vγ3 TCR (clone MG2), Vδ1/Vγ8 TCR (clone MG3) and Vδ1/Vγ4 TCR (clone MG4) clones represented 21%, 18%, 18%, and 12% of total γδ T cells, respectively (Supplementary Table S9). The oligoclonal repertoire indicated the possibility of antigen-driven clonal expansion, which was consistent with the cells expressing checkpoint markers and having an activated, resident phenotype.

We previously showed that gene transfer of the TCR expressed by the second most abundant clone (MG2) into a Jurkat-76 reporter T-cell line imbued TCR-mediated reactivity to the monomorphic antigen-presenting molecule MHC-related protein-1 (MR1), suggesting that the self-protein MR1 may be the molecular target for this tumor-infiltrating γδ T–cell clone (32). To extend this analysis, we generated three more Jurkat-76 lines expressing the TCRs from the MG1, MG3, and MG4 clones. These lines were challenged with C1R antigen-presenting cells (APC) engineered to express distinct candidate γδ TCR ligands: the MHC-like molecules MR1, CD1a, CD1b, CD1c and CD1d (Supplementary Fig. S6). As expected, the MG2 line reacted strongly to MR1. MG3 and MG4 did not reproducibly respond to any of the C1R APCs, albeit with some nonreproducible low-level basal responsiveness to C1R cells that varied between repeated assays. In contrast, MG1 responded strongly to CD1c-expressing APCs, and this response was blocked by an anti-CD1c mAb, confirming the reactivity. Thus, the two most frequent γδ T–cell clones observed in P055 were CD1c- and MR1-reactive, respectively, and these two clones collectively represented 39% of the γδ T cells and 11% of the total T cells in that tumor.

We next applied scRNA-seq to sorted CD3⁺ (2,651) and γδ⁺ (1,357) cells from tumor P055. All cells were combined for dimensionality reduction and UMAP clustering (Fig. 3A). UMAP cell clusters corresponded to CD8A-expressing T cells (C1), CD4-expressing T cells (C2), FOXP3-expressing T cells (C4), and two distinct populations of γδ T cells (C0 and C3). Differential gene-expression analysis was performed contrasting each cell cluster to the remaining clusters (P < 0.05, fold change > 0.25; Fig. 3B). The majority of γδ T cells clustered within C0 and were characterized by high expression of TCRD and TCRG genes (e.g., TRDV1, TRDC, and TRGC2); markers of innate immunity (NKG7, GLNY, KLRC3, and FCGR3A, which encodes CD16); markers of cytotoxic cell function, including granzymes (GZMA, GZMH, and GZMB), perforin (PRF1); and markers of T-cell exhaustion (PDCD1, which encodes PD-1, LAG3, and HAVCR2, which encodes TIM3; Fig. 3C). Cells in cluster C3 had high expression of CD69 and innate-like T-cell markers, but lacked the exhausted and tissue-resident characteristics of cluster C0.

To complement Sanger sequencing–based TCR analysis, we adapted the same primer reagents for γδ TCR-seq to be compatible with 10x Genomics. Consistent with the Sanger sequencing–based data, the γδ T cells were mostly Vδ1 and Vδ3 clonotypes (Fig. 3D). Furthermore, the same top four clonotypes represented 80% of all γδ T cells (Fig. 3E). In contrast, scTCR-seq of αβ cells revealed a diverse repertoire of αβ TCRs in the CD4⁺ and CD8⁺ populations and no dominant clones (Supplementary Fig. S7).

We next sought to identify a gene signature from the scRNA-seq data what would reflect the proinflammatory and innate profile of γδ T cells in MCC (Fig. 4A). Because γδ T cells and CD8⁺ T cells shared the most common gene-expression patterns, we contrasted γδ clusters (C0, C3) and CD8A-expressing cells (C1) after excluding other T-cell clusters (log fold change > 0.25, minimum percentage of cells expressing >30%). Genes highly expressed (logCPM > 0.5) in MCC cell lines (30) were removed to exclude genes potentially expressed by MCC tumor cells, which would confound downstream bulk tumor analysis. Genes not correlated (r_s < 0.6) in bulk RNA gene-expression data across 34 MCC tumor samples were also removed. γδ TCR constant regions (e.g., TRDC and TRGC2) were removed as we found they are expressed in other T- and NK-cell subsets based on published bulk RNA-seq data, but TRDV1 was retained as it is more specific and reflects the predominant Vδ1 infiltrate in MCC (Supplementary Fig. S8; ref. 31). Correlation between mRNA expression and the intratumoral DN T-cell count (detected by mIHC/IF) showed that 7 of 13 genes were ranked above the 90th percentile (Fig. 4B; Supplementary Table S10). The expression of the individual genes comprising the 13-gene γδ T–cell signature in the scRNA-seq data is shown in Supplementary Fig. S9.

GSVA was next used to calculate enrichment scores in MCC tumors subject to bulk RNA-seq (Fig. 4C). As expected, the GSVA 13-gene γδ T–cell signature score correlated with intratumoral DN T cells (mean cells/ROI) detected by m-IHC/IF (r_s = 0.6). Lower positive correlations were observed for CD8⁺ (r_s = 0.27) and FOXP3⁺ cells (r_s = 0.27), but no correlation to CD4⁺ cells (r_s = 0.02; Fig. 4D). Similar results were observed for total T-cell infiltrates with respect to T-cell subset. The cell-type specificity of the GSVA γδ T–cell enrichment was further tested using published bulk RNA-seq data for sorted adaptive and innate immune cell types, including CD8, CD4, NK, Vδ1, Vδ2, NKT, and MAIT cells (Supplementary Fig. S10; ref. 31). High enrichment scores were observed not only for Vδ1⁺ and Vδ2⁺ γδ T cells, but also for NK cells. NKTs, MAITs, CD4⁺ and CD8⁺ T cells all had negative median GSVA scores.

GSVA using the 13-gene γδ T–cell score was next applied to available gene-expression datasets with clinical annotation (Supplementary Table S11). Two datasets represented patients treated before the introduction of ICI therapies, including an unpublished RNA-seq cohort from Melbourne (RNA-Access/RNA-Direct; n = 39) and a published dataset representing patients treated in the United States and Europe (GSE22396; n = 29; ref. 11). A third series included RNA-seq data for patients treated with anti–PD-1/PD-L1 antibodies (NEB-Next; n = 21; ref. 29). GSVA was applied independently to each dataset, classifying tumors and patients into three groups based on the Q1/Q3 percentile rank thresholds (low <25th, medium 25th–74th, and high ≥75th; Supplementary Fig. S11). To maximize statistical power for survival analysis we combined all data representing 89 MCC patients. A Kaplan–Meier analysis showed a significant disease-specific survival (DSS) difference between groups, with patients harboring tumors with high γδ T–cell enrichment scores having a longer DSS compared with patients with tumors with low and medium enrichment scores (P = 0.034; Fig. 4E). Cox regression analysis was performed using the GSVA score intervals, sex, age at diagnosis, MCPyV status, American Joint Committee on Cancer (AJCC) stage at diagnosis, and immunotherapy treatment (n = 83; Fig. 4F). High enrichment for the γδ T–cell gene signature showed a low HR (HR = 0.26, P = 0.021) suggesting a 74% decrease in risk of death.

Four patients treated with ICI therapies had high tumor γδ T cell–enrichment scores, of whom three had a complete response to treatment. Case P007 represented VN-MCC of unknown primary presenting initially in a submental lymph node and then serially as metastatic lesions in the axilla, near the lip and finally in a retroperitoneal lymph node (Fig. 5A and B). All of these tumors were surgically resected. Following the removal of a retroperitoneal tumor mass, the patient received adjuvant therapy with the anti–PD-1 therapy pembrolizumab; the patient remained disease free following 3 years of continuous treatment. Archival FFPE tissue was available for three MCC lesions, including the lip, axillary node, and retroperitoneal mass. m-IHC/IF showed a prominent T-cell infiltration and PD-L1 staining at the invasive tumor margin in the retroperitoneal mass (Fig. 5C). Despite the discordance in total T-cell infiltrate between P007 tumors, the percentage of DN T cells as a fraction of total CD3⁺ T cells was similar (49%–70% DN T cells). Flow cytometry analysis of the retroperitoneal metastasis confirmed the presence of DN T cells (Fig. 5D). Moreover, IHC targeting TCRD confirmed heavy infiltration of γδ T cells (Fig. 5E). Similar to our observation in tumor P055, flow cytometry showed DN T cells had high PD-1 expression (Fig. 5F).

T-cell infiltration is a common feature of MCC tumors and associated with improved survival, reflecting the immunogenic properties of this disease (11–13). Here, we confirmed the abundance of infiltrating CD4⁺ and CD8⁺ T cells in MCC and identified a large fraction of T cells without CD4 or CD8 cell-surface expression. This feature would have been missed by previous IHC-based studies as they did not use multiparameter methods. Flow cytometry revealed that γδ T cells represented the majority of DN T cells in MCC tumors. In the context of T-cell inflammation, γδ T cells can express markers of T-cell activation, exhaustion and tissue residency, suggestive of an antigen-driven response. We showed that γδ T cells in MCC had a proinflammatory and innate gene-expression profile predictive of better patient survival. These results indicate for what we believe to be the first time that γδ T cells are likely an important effector cell type in MCC tumors.

Tumor-infiltrating γδ T cells have been described in several other cancer types. γδ T–cell infiltrates in melanomas and cutaneous squamous cell carcinomas (cSCC) are higher in tumors than healthy skin or tonsil and associated with a better prognosis (21, 23, 40). There is conflicting data relating to γδ T–cell infiltrates in breast, colorectal, and pancreatic cancers, with variable prognostic significance reported (24–27). Like MCC, other skin cancers have a predominance of Vδ1⁺ cells, although in melanoma, an increased proportion of Vδ2⁺ over Vδ1⁺ cells are associated with improved survival (21). In cSCC, Vδ2⁺ cells were found to be proinflammatory and IFNγ-expressing, whereas Vδ1⁺ cells have a predominantly Th17 phenotype (23). In MCC, we detected a dominance of Vδ2⁻ γδ T–cell infiltrates, but these cell harbored an innate gene-expression pattern, similar to that reported in triple-negative breast cancers (26). Although we did not stimulate γδ T cells from MCC tumors to measure cytokine expression, these cells had detectable IFNG gene expression, and it is known that in the resting state, IFNG mRNA expression in T cells is predictive of IFNγ secretion upon stimulus (31). These data suggest that most γδ T cells in MCC have a proinflammatory function.

Several prior studies have described methods for in silico prediction of γδ T–cell tumor infiltrates. Indeed, the CIBERSORT method showed that γδ T–cell enrichment was the best predictor of survival across 18,000 tumors of 39 cancer types (41). However, the original LM22 reference used for CIBERSORT was criticized for lack of specificity owing to colinearity of γδ T cells with other immune cell types, leading to efforts to improve the accuracy of reference datasets (42). A reliance on reference datasets representing cells from blood is another potential pitfall as functionality of these cells may be different to their counterparts in tumors. Furthermore, focusing on specific Vδ1⁺ or Vδ2⁺ subsets may limit the signature specificity given their different gene-expression profiles (43). In this study, we profiled γδ T cells isolated directly from tumor rather than peripheral blood. The 13-gene γδ T–cell GSVA score correlated with DN T cells by mIHC/IF, consistent with our assumption that the majority of DN T cells in MCC are γδ T cells. However, it is clear that individual genes comprising the 13-gene signature may also be expressed by other innate cell types, such as NK cells. Previous quantification of NK infiltrates in MCC relied upon CD16 IHC staining (12, 44). That said, CD16 is not a specific marker of NK cells (45), and indeed we showed that the gene encoding CD16 (FCGR3A) was expressed by a substantial fraction of γδ T cells and correlated with DN T cells (Supplementary Table S9). However, given that we cannot exclude the contribution of NK cells to the 13-gene γδ T–cell score, the interpretation of the scoring in individual cases should be treated with the appropriate level of caution.

Outstanding questions remain regarding the antigens detected by γδ T cells. γδ T–cell infiltration does not correlate with MCC viral status, suggesting that antigens common to both VP-MCC and VN-MCC are recognized by γδ T cells. We observed clonal expansion of Vδ2⁻ (Vδ1⁺ or Vδ3⁺) clonotypes by scTCR-seq, which, combined with the expression of activation and tissue retention markers and immune checkpoint molecules, is consistent with antigen-driven expansion. Subsets of Vδ2⁻ cells are known to detect a range of monomorphic MHC class I–like molecules including CD1d and CD1c. Here, we demonstrate functional recognition of CD1c by the most abundant clone within the tumor P055, MG1. We previously reported a novel class of γδ T cells that recognize the MHC-like molecule MR1 (32). In that study, we briefly reported the second most abundant TCR clonotype (MG2) derived from MCC lesion P055 was MR1 reactive. The role of MR1 and CD1c in the context of MCC is unclear. Both are antigen-presenting molecules that present small metabolites and lipids, respectively. CD1c is typically expressed by a range of myeloid and lymphoid cell types and can present leukemia-associated lipid antigens (46). MR1 is expressed by all nucleated cells, and appears to present as-yet-undefined tumor-associated antigens from a wide range of tumor types (47, 48). It is thus conceivable that MR1 expression on MCC cells or MR1/CD1c-expressing tumor-infiltrating cells are presenting tumor-associated antigens for surveillance by γδ T cells. Accordingly, MR1 and CD1c may represent novel immunotherapeutic targets, and further studies on this axis are clearly warranted.

The prospect that γδ T cells can be a target of existing ICI therapies is highly relevant to contemporary treatment of patients with MCC, although our understanding is currently limited. We reported complete response to ICI treatment in three patients with MCC where tumors had high γδ T–cell enrichment scores. There is evidence that γδ T cells may influence response to ICI in other diseases. In melanoma patients treated with anti-CTLA4 (ipilumimab), total γδ T cells in peripheral blood had no prognostic relevance but better patient survival was observed on the basis of increased frequency of Vδ2⁺ over Vδ1⁺ subsets (22). This is consistent with Vδ2⁺ cells being cytotoxic effectors in melanoma. Similar to our observations in MCC, γδ T cells found in melanomas and patient blood commonly express immune-checkpoint molecules including LAG3, TIM3, and PD-1 (21). Functional experiments to show the direct effects of ICI on γδ T cells is currently limited but PD-1 blockade has been reported to upregulate IFNγ secretion by γδ T cells in a leukemia model (49). Accordingly, γδ T cells are likely targeted by ICI but the functional outcomes require further investigation. The detection of γδ T cells or a gene-set enrichment score for γδ T cells, as described in this study, should be explored in future biomarker analysis for ICI-responsive MCC.

In summary, we have described γδ T cells as an abundant effector cell in MCC, thereby expanding our understanding of the immune response to this highly immunogenic disease. Future work will require analysis of a larger independent series of MCC and deeper characterization of the function of γδ T cells to fully understand their prognostic and therapeutic predictive value. Given the cytotoxic and proinflammatory effector function of γδ T cells, together with their lack of MHC class I restriction, which makes them impervious to tumor MHC class I downregulation, targeting these cells may represent an alternative approach to cellular therapeutics through adoptive or allogeneic transfer into patients with MCC cancer, and it is an approach being explored in other cancers (50).

N.A. Gherardin reports grants from Australian Research Council and National Health and Medical Research Council during the conduct of the study. A. Caneborg reports grants from Medical Research Future Fund during the conduct of the study. K.F. Shannon reports personal fees from Merck Serono outside the submitted work. A. Hong reports personal fees from Bayer outside the submitted work. P.J. Neeson reports grants from Bristol-Myers Squibb, Roche/Genentech, Compugen, Allergan, and Merck Sharp & Dohme outside the submitted work. J.C. Becker reports personal fees from Merck Serono, Pfizer, and Sanofi, and grants from 4SC and IQVIA outside the submitted work. S. Sandhu reports grants from AstraZeneca, Merck Sharp & Dohme, Merck Serono, Genentech, Amgen, and Novartis and personal fees from Merck Sharp & Dohme, Bristol-Myers Squibb, and AstraZeneca outside the submitted work. G.A. McArthur reports other support from Genentech/Roche and Array BioPharma during the conduct of the study, as well as other support from Genentech/Roche and Array BioPharma outside the submitted work. R.A. Scolyer reports grants from National Health and Medical Research Council during the conduct of the study, as well as personal fees from Qbiotics, Novartis, Merck Sharp & Dohme, NeraCare, Amgen, Bristol-Myers Squibb, Myriad Genetics, and GlaxoSmithKline outside the submitted work. R.J. Hicks reports other support from Telix Pharamceuticals outside the submitted work. D.I. Godfrey reports grants from National Health and Medical Research Council and Australian Research Council during the conduct of the study, as well as other support from CSL, grants from Cancer Council Victoria, personal fees from Avalia Immunotherapies, and a patent for “Compositions and uses thereof” pending outside the submitted work. R.W. Tothill reports grants from National Health and Medical Research Council, Victorian Cancer Agency, and Medical Research Future Fund during the conduct of the study. No disclosures were reported by the other authors.

N.A. Gherardin: Formal analysis, investigation, methodology, writing–original draft, writing–review and editing. K. Waldeck: Data curation, formal analysis, investigation, writing–review and editing. A. Caneborg: Formal analysis, writing–review and editing. L.G. Martelotto: Formal analysis, methodology, writing–review and editing. S. Balachander: Formal analysis, writing–review and editing. M. Zethoven: Formal analysis, writing–review and editing. P.M. Petrone: Formal analysis, writing–review and editing. A. Pattison: Formal analysis, supervision, writing–review and editing. J.S. Wilmott: Formal analysis, writing–review and editing. S.M. Quiñones-Parra: Formal analysis, writing–review and editing. F. Rossello: Formal analysis, writing–review and editing. A. Posner: Formal analysis, writing–review and editing. A. Wong: Formal analysis, writing–review and editing. A.M. Weppler: Data curation, writing–review and editing. K.F. Shannon: Data curation, writing–review and editing. A. Hong: Data curation, writing–review and editing. P.M. Ferguson: Data curation, writing–review and editing. V. Jakrot: Data curation, writing–review and editing. J. Raleigh: Data curation, writing–review and editing. A. Hatzimihalis: Data curation, investigation, writing–review and editing. P.J. Neeson: Resources, supervision, writing–review and editing. P. Deleso: Data curation, writing–review and editing. M. Johnston: Data curation, writing–review and editing. M. Chua: Data curation, writing–review and editing. J.C. Becker: Resources, writing–review and editing. S. Sandhu: Resources, data curation, writing–review and editing. G.A. McArthur: Resources, writing–review and editing. A.J. Gill: Resources, formal analysis, writing–review and editing. R.A. Scolyer: Resources, formal analysis, writing–review and editing. R.J. Hicks: Resources, supervision, writing–review and editing. D.I. Godfrey: Resources, formal analysis, supervision, writing–review and editing. R.W. Tothill: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.

The authors are grateful to following people: Dr. Paul Savage (Brigham Young University, Provo, UT, USA) for providing α-galactosylceramide analogue PBS44 used for production of CD1d-α-GalCer tetramers; Prof. David Fairlie (University of Queensland, Brisbane, Queensland, Australia) for providing 5-ARU used for production of MR1-5-OP-RU tetramers; staff from the flow cytometry facilities at the Department of Microbiology and Immunology at the Peter Doherty Institute and the Melbourne Brain Centre, both at the University of Melbourne, as well as the Peter MacCallum Cancer Centre; University of Melbourne Centre for Cancer Research and Peter MacCallum Cancer Centre Genomics Core Facilities and the Genomics Platform Group at University of Melbourne for their support in generating primary data; Stephen Fox at the Peter MacCallum Cancer Centre for patient samples; Paul Nghiem, Kelly Paulson, and Kristina Lachance at the Fred Hutchinson Cancer Research Center for providing clinical annotation for their published MCC microarray data. This work was supported by the National Health and Medical Research Council of Australia (NHMRC; 1113293) and the Australian Research Council (ARC; CE140100011). D.I. Godfrey was supported by an NHMRC Senior Principal Research Fellowship (1020770, 1117766). R.A. Scolyer is supported by an NHMRC Program Grant, Practitioner Fellowship (1141295), the Ainsworth Foundation, the ClearBRIDGE Foundation, the Cameron Family, and colleagues at Melanoma Institute Australia and Royal Prince Alfred Hospital. R.J. Hicks is supported by an NHMRC Practitioner Fellowship (1108050). R.W. Tothill is supported by a Victorian Cancer Agency Mid-Career Fellowship and has received funding from the Melanoma and Skin Service group at Peter MacCallum Cancer Centre, the University of Melbourne (Research Support Grant Scheme), and the Medical Research Future Fund (1169955). A. Pattison is supported by the Joseph Herman Trust at the University of Melbourne.

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.