Preexisting Skin-Resident CD8 and γδ T-cell Circuits Mediate Immune Response in Merkel Cell Carcinoma and Predict Immunotherapy Efficacy

Abstract

Merkel cell carcinoma (MCC) is an aggressive neuroendocrine skin cancer with a ∼50% response rate to immune checkpoint blockade (ICB) therapy. To identify predictive biomarkers, we integrated bulk and single-cell RNA sequencing (RNA-seq) with spatial transcriptomics from a cohort of 186 samples from 116 patients, including bulk RNA-seq from 14 matched pairs pre- and post-ICB. In nonresponders, tumors show evidence of increased tumor proliferation, neuronal stem cell markers, and IL1. Responders have increased type I/II interferons and preexisting tissue resident (Trm) CD8 or Vδ1 γδ T cells that functionally converge with overlapping antigen-specific transcriptional programs and clonal expansion of public T-cell receptors. Spatial transcriptomics demonstrated colocalization of T cells with B and dendritic cells, which supply chemokines and costimulation. Lastly, ICB significantly increased clonal expansion or recruitment of Trm and Vδ1 cells in tumors specifically in responders, underscoring their therapeutic importance. These data identify potential clinically actionable biomarkers and therapeutic targets for MCC. Significance: MCC serves as a model of ICB response. We utilized the largest-to-date, multimodal MCC dataset (n = 116 patients) to uncover unique tumor-intrinsic properties and immune circuits that predict response. We identified CD8 Trm and Vδ1 T cells as clinically actionable mediators of ICB response in major histocompatibility complex-high and -low MCCs, respectively.

Fig. 1:. Differentially expressed genes identify key immune and tumor signatures distinguishing immunotherapy responders and non-responders.

A, Patient sample heatmap of assays and clinical metadata for each sample in our dataset. Each column represents a single patient with the available assays and associated clinical characteristics. In total, 116 patients were included in this study. 101 patients received immunotherapy. 71 patients were assayed by bulk RNA-seq, 64 by single cell (sc)RNA-seq, 41 by multiplex immunofluorescence (mIF), and 20 by spatial transcriptomics (GeoMx and CosMx). VP, Virus-positive; VN, Virus-negative B, Schematic of the study overview. Patient characteristics, paraffin-embedded (FFPE) sample sequencing breakdown with sample counts, single-cell processing workflow with sample counts, and single-cell library breakdown with sample counts. Fifty-four tumor samples and 55 peripheral blood mononuclear cell (PBMC) samples were stained with 45 antibody-derived tags (ADT) and enriched for live, CD45+ cells (tissue), or CD3+ cells (PBMCs) via fluorescence-activated cell sorting (FACS). Cells were processed using 10x Chromium and four sequencing libraries were prepared for each sample: 5’ cDNA, ADT, ⍺β VDJ, and γδ VDJ. C, Volcano plot of differentially expressed (DE) genes between immunotherapy (IO) responders (R) and non-responders (NR). Bulk RNA-seq was conducted on 63 FFPE tumor samples from 62 immunotherapy naïve, systemic therapy-naïve patients, all samples prior to evaluable responses to ICB. Differentially expressed genes were identified using the DESeq2 package, with the following covariates: sex, viral status, and tissue. Selected protein-coding genes with a false discovery rate (FDR) < 0.1 were labeled. A total of 800 genes were upregulated in responders, and 436 genes were upregulated in non-responders. D, Pathway analysis of differentially expressed (DE) genes in responders. Only protein-coding genes with FDR <0.1 were included in the analysis. The top five enriched gene sets from the selected gene set libraries are shown for the responder DE genes. Gene sets were ranked based on the Enrichr combined score, which is a measure of both Fisher’s exact test p-value and the z-score of the deviation from the expected rank of the gene set. E, Three gene sets enriched in responders with associated genes are represented on the volcano plot of DE genes from responders and non-responders. Genes with FDR < 0.1 are labeled. F, Enrichment gene list analysis of DE genes found in non-responders. Only protein-coding genes with FDR <0.1 were included in the analysis. The top five enriched gene sets from the selected gene set libraries are shown for non-responder DE genes. Gene sets were ranked based on the Enrichr combined score. G, Two gene sets enriched in non-responders with associated genes are represented on the volcano plot of DE genes from responders and non-responders. Only genes with FDR < 0.1 are labeled.

Fig. 2:. Response to ICB is associated with tumor cell upregulation of an interferon-γ signature and negatively correlated with a tumor cycling signature.

A, Uniform Manifold Approximation and Projection (UMAP) projection of 139,845 cells among all 49 tissue samples in single-cell (sc) transcriptome datasets with major cell-type annotation. PCA inputs were batched-corrected using the Harmony package and cells were clustered using the standard Seurat pipeline. B, Dot plot showing the scaled expression of the top five marker genes for each cell type annotated in the tissue scRNA-seq dataset. C, Differential gene expression analysis of 60,902 tumor cells pseudobulked by tissue sample between ICB responders and non-responders was performed using the edgeR LRT test with sex and virus status as a covariate. 34 genes were upregulated in responders, and 25 genes were upregulated in non-responders. Protein-coding genes with FDR <0.25 are labeled. D, Gene set enrichment analysis (GSEA) of DE genes found in immunotherapy responders (R) and non-responders (NR). Hallmark gene sets were ranked by normalized enrichment score (NES), with the top 5 gene sets from responders and non-responders listed. *** p adjust<0.001, ** p adjust< 0.01, * p adjust< 0.05 E, SCENIC regulon analysis was run on 60,902 tumor cells. The top 5 regulons differentially expressed in responders (R) and non-responders (NR) are shown, ranked by false discovery rate (FDR), Wilcoxon test. F, Scatter plot of CD45+ infiltration (% CD45+ cells from FACS) with logit transformation versus Hallmark Interferon-γ signature module score in tumor cells from scRNA-seq, which depicts a positive Pearson correlation between interferon- γ response in the tumor and overall immune infiltration. Samples are colored according to ICB response. IFNG score was calculated on tumor cells by the ModuleScore function in Seurat using the MSIGDB “Hallmark Interferon Gamma Response” genes as input. G, A violin plot of IFNG expression of each cell type split by samples and ranked by mean expression. H, Scatter plots of immune proportion (% CD45+ cells from FACS) with logit transformation or tumor cell IFNG score compared to the proportion of proliferating tumor cells among the total cells with log transformation in each scRNA-seq tissue sample. Samples are colored by ICB response, which shows a negative Pearson correlation between the proportion of proliferating tumor cells and either the immune infiltration or IFNG response in the tumor cells, with non-responders clustered in the high proliferation, low IFNG or low CD45+ % quartile. I, A Kaplan-Meier curve of post-immune checkpoint blockade (ICB) survival shows that samples with higher levels of tumor cycling cell signature (derived from scRNA-seq) had lower survival. Samples from the bulk RNA-seq dataset were scored based on the “tumor cycling” cell type signature derived from our single cell data with gene set variance analysis (GSVA) and split into the top or bottom quartile (see Methods). Statistical test: Cox proportional hazards.

Fig. 3:. CD8 TILs form two distinct developmental trajectories composed of circulatory T cells and Tissue resident memory subsets.

A, Uniform Manifold Approximation and Projection (UMAP) projections of 24,448 CD8 T cells in 15 cell clusters. Outlines of the clusters show the formation of two meta-clusters: circulating T cells (Tcirc) and tissue-resident memory (Trm) T cells. Cells were subset from all tissue scRNA-seq samples and clustered using the standard Seurat pipeline. B, Heatmap of the antibody-derived tag (ADT) protein average scaled expression of each CD8 cluster with hierarchal clustering into two major meta-clusters distinguishing Tcirc and Trm cells. C, Bar chart of the proportion of co-expressing, co-inhibitory markers statistically enriched on T resident (Trm) compared to T circulatory (Tcirc) cells (χ², p <2.2e-16). D, Dot plot of selected average marker gene expression for each CD8 lineage with functionally grouped genes. logFC: log2 fold change E, Diffusion mapping projection of CD8 T cells with cell cycle genes regressed out using the Destiny package. The diffusion map shows two separate lineages, which are colored by Trm or Tcirc. Cells in the CD8 cluster C04 are depicted in gray. C04 show characteristics of both Trm and Tcirc cells. DC: diffusion component F, Heatmap of the log-likelihood ratio (LLR) of shared TCR clones between CD8 clusters with the heatmap clustered by LLR demonstrates the formation of two meta-clusters that match the Trm and Tcirc lineages. *, p <0.05, Fisher’s test G, Bar chart of all expanded CD8 clonotypes in the tissue dataset. Each line represents the proportion of cells in each lineage of the single clonotype. The clonotypes were ranked from the highest proportion of cells in the Trm lineage to the lowest proportion (highest in the Tcirc lineage). For 88% of the clonotypes, 100% of the cells were restricted to either the Trm or Tcirc lineage. H, Bar chart of the proportion of clones that were also found in the PBMC single-cell dataset split by developmental lineage. Only matched blood samples were included in the analysis. Chi-square test for significant differences between proportions. I, UMAP projection (same as panel A) colored by the TCR clone size. J, Violin plots of clonal size of all clonotypes found in the Trm or Tcirc lineage.

Fig. 4:. CD8 Trm cells are tumor specific and prognostic of ICB response.

A, Violin plot of cellular expression of the NeoTCR8(23) module score split by CD8 cluster. Clusters are colored by developmental lineage: T circulatory cells (Tcirc) or T resident cells (Trm). B, Cell density plot showing the distribution of cells with the MO_MCC_095c1 TCR clonotype, a highly expanded Trm TCR clone found in virus-positive MCC samples. Outline shows the distribution of Tcirc and Trm cells. C, Jurkat-76 cells modified with an NFAT-GFP reporter and CD8α/CD8β were transduced and selected to express the MO_MCC_095c1 TCR. These cells were co-cultured with an autologous lymphoblastoid cell line loaded with pooled peptides from the small and large T antigens of Merkel cell polyoma virus. After co-culture, GFP expression on the reporter cells was measured as a measure of antigen recognition, leading to T cell activation. CEF (CMV, EBV, Influenza) peptide pool, pool 1, and pool 2, pool 3 are representative examples of negative hits, whereas pools 8 and 12 are positive hits. See Supplementary Fig. S14 shows the complete results for all pools. *, p <0.05, One-way ANOVA. D, Structural model of the MO_MCC_095c1 TCR in complex with its target peptide, small T (sT) antigen residues 31–39, bound to its putative MHC allele, HLA B*07:02, reveals a canonical TCR-pMHC binding mode and multiple polar contacts that stabilize the interaction (dashed lines). TCRα is magenta, TCRβ cyan, peptide slate (stick representation), and MHC green. The high confidence of the model by AlphaFold score (0.82 out of 1) indicates a high likelihood of structural accuracy. E, Scatter plot of all significant GLIPH2(26) complementary-determining region 3 (CDR3) motifs found in the CD8 tissue TCR dataset. Axis has been adjusted to show motifs that are >80% found in virus-positive samples and >80% found in Trm cells. These virus-positive (VP) Trm motifs are shared between multiple patients and are highly expanded and represent likely shared anti-MCC TCR motifs. “%” indicates a variable amino acid found in the GLIPH motif. F, Representative multiplex immunofluorescence (mIF) images show skin resident CD8 T cells (CD3+CD8+CD103+) infiltrating the tumor parenchyma, marked by pan-cytokeratin (PanCK), and are cytotoxic, marked by granzyme B (GZMB). Scale bar:100 μm. G, A scatter plot of tissue-resident CD8 T cell densities (percentage of CD3+CD8+ expressing CD103 per mm²) from 29 mIF samples from 26 patients (17 responders, 9 non-responders). Wilcoxon test. H, A violin plot of the ratio of Trm to Tcirc cells split by response to immunotherapy (IO). Wilcoxon test. I, Kaplan-Meier curve of survival from the start of ICB from the bulk RNA-seq dataset. Samples were split based on the top and bottom quartiles of the Trm gene set variation analysis score derived from our single-cell dataset (see Methods section). Statistical test: Cox proportional hazards.

Fig. 5:. Landscape of γδ T cells and NK cells in MCC mirrors Trm and Tcirc trajectories found in CD8 T cells suggesting functional convergence in T-cell driven responses.

A, Scatter plot of the interferon-γ response score (module score from the hallmark gene set) of tumor cells versus the proportion of γδ T cells in a sample that shows the subset of γδ T cell-enriched samples that also have a low IFNG score. B, UMAP projection of 9722 γδ/NK cells among 16 clusters from tissue scRNA-seq data. Cells were clustered using the standard Seurat pipeline. C, Differentially expressed genes (pseudobulked, edgeR LRT test) in γδ/NK cells compared to samples with low tumor HLA I expression to samples with high tumor HLA I expression. D, A chord diagram of γ and δ V chain pairings for all cells in the tissue scRNA-seq dataset. E, Dot plot of average expression of selected significant, Wilcoxon gene markers of γδ and NK cells in both the tissue datasets separated by developmentally distinct δ-chain or NK lineage. δ1 and δ3 T cells show higher levels of exhaustion markers, activation markers, killer cell immunoglobulin-like receptor (KIR) family members, cytotoxic effectors, CXCR6, adhesion molecules, and exhaustion-related transcription factors (TF). δ2 T cells express high levels of TNF, LTB, naïve/T central memory (Tcm) markers, and AP-1 family members (FOS). F, Violin plot of selected significant Wilcoxon protein markers separated by tissue/peripheral blood mononuclear cell (PBMC) samples and grouped by δ-chain or NK lineage. G, A volcano plot of differentially expressed protein markers (ADT) of Vδ1 cells versus Vδ2 cells. Vδ1 cells upregulated exhaustion and activation markers. Vδ2 cells upregulated naïve and T central memory markers. All labeled proteins had P_adj. < 0.05. Wilcoxon test. H, Multiple alignment and consensus sequence of δ CDR3 of the top expanded Vδ1 clones shows clonal focusing with sequence similarity between patients. See Methods for details on TCR selection. I, UMAP projection of Vδ1 cells (same projection as panel b) colored by pseudotime as calculated by Monocle demonstrate that “root” cells are found mainly in C08_Vd1 whereas differentiated cells are found in C02_Vd1. Root cells were selected based on their highest expression of naïve-like markers, including TCF7, CD45RA, CD127, and CCR7. J, Heatmap of genes that are significantly and non-linearly corelated with Vδ1 pseudotime. Vδ1 T cells have high expression of naïve markers like TCF7 as well as AP-1 family members at the beginning of pseudotime. During differentiation, Vδ1 cells upregulate co-inhibitory, cytotoxic, activating, and Trm makers. Genes are colored by scaled expression across pseudotime. K, Representative image of multiplex immunofluorescence (mIF) showing cytotoxic granzyme B (GZMB+) γδ T cells δ TCR (dTCR+) infiltrating the tumor parenchyma, marked by pan-cytokeratin (PanCK), as indicated by arrows. Scale bar: 50 microns.

Fig. 6:. Spatial and temporal analysis of tumor infiltrating lymphocytes

A, Violin plots of the relative expression of GeoMx regions of interest (ROI) of Tissue resident memory (Trm)/circulatory T (Tcirc) signature genes (derived from our scRNA-seq dataset) by module score in the stroma or tumor parenchyma. B, UMAP projection of 614,209 cells segmented using the CosMx assay. The cells are colored according to the inferred cell type. Cells were clustered in Seurat and batch-corrected using SCTransform. C, Unsupervised hierarchical clustering heatmap of cell type co-localization as measured by the degree of co-localization (DoC) within CosMx fields of view (FOV). Neighboring cells were identified within a 54 micron radius of each cell, and then the (DoC) was determined using the SpatialTIME(34) method. A high degree of co-localization was observed among T cells, B cells, DCs, and pDCs, is outlined in the upper left. Ave: average; DoC: degree of co-localization. Mac: Macrophages, Endo: endothelial cells, Fibro: fibroblasts, Tumor cyc: tumor cycling cells, pDC: plasmacytoid dendritic cells, DC: dendritic cells, Epi: epithelial cells, Mono: monocytes D, A sample FOV showing co-localization of T cells, B cells, and pDC at the tumor-stroma interface using CosMx cell segmentation representation (left) and morphological markers showing para-nuclear dot-like expression of Pan cytokeratins (CK), with high CD298 expression across immune cells. Scale bar: 50 μm. E, Bar charts of the cell types found in T cell neighborhoods (54 μm radius) showing prevalence of pDC, DC, and B cells in T cell neighborhoods of responders (R) versus T cell neighborhoods of non-responders (NR). Comparisons were made using the chi-square test, ****, p value < 2.2e-16. F, A dot plot of receptor-ligand interactions from CellChat on the scRNA-seq dataset identifies multiple activating chemokines, cytokines, and surface receptor interactions between B cells, DCs, and pDCs to exhaustion, CD8 Trm cells (CD8 Tex.). G, Differential gene expression analysis (DE-seq2) comparing matched pre- versus post-ICB sample pairs across 7 non-responders and 7 responders shows upregulation of inflammatory genes after ICB treatment. Among responders, 1560 genes were upregulated and 995 genes were downregulated after therapy. (DE-Seq2, P_adj. <0.1,) Only two differentially expressed genes were identified among the non-responders. H, Enrichr gene list analysis of DE genes in pre- versus post-ICB matched pairs from responders shows inflammatory pathways upregulated and cell cycle pathways downregulated after treatment. Only protein-coding genes with FDR <0.1 were included in the analysis. The top five enriched gene sets from the Hallmark gene set libraries are shown for the DE genes. Gene sets were ranked based on the Enrichr combined score, which is a measure of both Fisher’s exact test p-value and the z-score of the deviation from the expected rank of the gene set. I, Corrected cellular proportions of Trm cells in 15 post-ICB treated samples show CD8 Trm infiltration is higher (3.01-fold) in partial responders (PR) versus progressive disease (PD) after treatment. Samples are split by response prior to biopsy. Five of the samples represent a new cohort of post-ICB samples used to verify results from the original scRNA-seq dataset. Wilcoxon test. J, Violin plots of CD8 cells in MO_MCC_079 and MO_MCC_085, two scRNA-seq pre/post sample pairs. Exhaustion score expression was split by treatment timeline and grouped by clonotype expansion showed that expanded cells were more likely to be exhausted pre-treatment, with exhaustion levels significantly decreased after treatment. The exhaustion score(37) was derived from Seurat AddModuleScore.