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

Merkel cell carcinoma (MCC) is a highly aggressive primary neuroendocrine skin cancer, with high rates of metastasis and recurrence (1). There are two subtypes of MCC, with different etiologies and phenotypes. In 60% to 80% of MCCs, the Merkel cell polyomavirus (MCPyV) is mutated and clonally integrated (2). They express obligate viral oncoproteins, small and truncated large T antigens, which drive carcinogenesis in part by inactivating p53 and pRB (3). While virus positive (VP) MCCs inactivate pRB and p53 through the expression of the small and large T oncoproteins, virus negative (VN) MCCs achieve the same molecular endpoint through near universal somatic mutation of both RB1 and TP53 (2, 4).

VN-MCCs are immunogenic because of the high number of tumor neoantigens from high mutation burdens. VP-MCCs have a low mutation burden but are immunogenic because they constitutively express viral oncoproteins (4). Shared antigens in VP-MCCs ease identification of tumor-specific T cells across patients (5).

For metastatic disease, immunotherapy is the only Food and Drug Administration–approved therapy that significantly prolongs survival. Unfortunately, many patients have either primary (∼50%) or acquired resistance (up to 25%) to immune checkpoint blockade (ICB; refs. 6, 7). To elucidate potential mechanisms, we compiled the largest MCC cohort with comprehensive multimodal molecular analysis (n = 116 patients), specifically focusing on the tumor microenvironment with bulk RNA sequencing (RNA-seq; n = 71 patients), single-cell RNA-seq (scRNA-seq; n = 360,069 cells from 64 patients), multiplexed immunofluorescence (n = 41 patients), and spatial transcriptomics (n = 20 patients). A total of 46 patients had more than one sample. These methods complement each other, addressing limitations and offering insights across tumor-intrinsic and -extrinsic factors. Our multiomic studies identified clinically actionable biomarkers and therapeutic avenues for MCCs and potentially for other cancers.

We collected the largest multi-institutional MCC dataset with multimodal molecular characterization of 186 samples across 116 patients including 92 patients with evaluable ICB responses (Table 1; Fig. 1A; Supplementary Tables S1 and S2; “Methods”). These include bulk RNA-seq (n = 85 samples from 71 patients; Supplementary Tables S3–S5; Supplementary Fig. S1A) and a modified form of ECCITE-Seq (n = 54 tumor samples from 48 patients; n = 55 blood samples from 32 patients; Supplementary Tables S6 and S7; Supplementary Fig. S1B), including three samples previously reported (8). Moreover, 44 distinct samples were spatially profiled using multiplex immunofluorescence (mIF) imaging (Akoya Phenoimager), spatial cancer gene or whole transcriptome profiling (Nanostring GeoMx), and single-cell resolution spatial transcriptomics (Nanostring CosMx; Fig. 1B; Supplementary Tables S8–S10; Supplementary Fig. S1C). Overall, 63 patients in our dataset were profiled using two or more modalities (Supplementary Fig. S2).

First, we analyzed the bulk RNA-seq of tumors, pre-ICB, from a cohort of 63 patients with 62 patients that were systemic treatment naive to minimize confounding variables related to prior therapy. Within this dataset, there were 31 responders (partial or complete response) and 32 nonresponders (stable and/or progressive disease; Supplementary Fig. S1A; Supplementary Table S3; “Methods”).

Differential gene expression analysis revealed 800 upregulated genes in responders and 436 in nonresponders (P adj < 0.1, DESeq2.0; Fig. 1C; Supplementary Table S11). Many of the genes are characteristic of CD8 T cells [αβ T-cell receptor (TCR), GZMA, CD8A, CD8B], γδ T cells (γδ TCR), B cells (IG genes, CD79A), and plasmacytoid dendritic cells (TLR7, LILRB4). Several molecules were associated with tumor-specific cytotoxic T cells, including markers of T-cell exhaustion (HAVCR2, LAG3, and TIGIT) and cytotoxic activity (GZMB and PRF1). Collectively, these data suggest an enrichment of preexisting CD8 T cells, B cells, and specific myeloid subtypes in responders.

Pathway analysis with Enrichr (9) showed upregulation of IFNγ signaling (P_adj. = 2.22e−63), IFNα signaling (P_adj. = 8.53e−52), and antigen presentation machinery (P_adj. = 2.36e−17), associated with IFNγ signaling in other cancers including melanoma (Fig. 1D and E; Supplementary Table S12; ref. 10).

Nonresponders had upregulated MCC-specific genes, like PAX6 and FOXA1 (Supplementary Table S11). Pathway analysis revealed additional enrichment of genes associated: (i) IL1 processing (11), including IL1A and IL1B, and (ii) E2F targets, which are normally expressed during cell division (Fig. 1F and G; Supplementary Table S12).

Next, we conducted scRNA-seq initially on 49 tumor samples derived from 43 patients (Supplementary Table S6; Supplementary Fig. S1B). This included eight samples from four patients across two timepoints (including two pairs of samples collected before and after ICB), four samples from two patients spanning different tissues (skin/lymph node pairs), and 37 single samples. Additionally, we analyzed 55 peripheral blood mononuclear cell (PBMC) samples from 32 patients. This included 21 pre/post-ICB pairs (42 samples), two samples from a third timepoint, and 11 samples that were unmatched. Twenty-six of these samples were matched to patients that also had scRNA-seq on tumor tissue (Fig. 1A and B; Supplementary Figs. S3–S5; Supplementary Table S7). A total of 313,757 cells were recovered from these samples using scRNA-seq. An additional cohort of five post-treatment tumor samples with 46,312 cells were integrated and analyzed separately to investigate the cellular landscape after ICB for a total of 360,069 cells. Following quality control, we identified 27 distinct cell types in the tumor immune microenvironment, including CD8 T cells, CD4 T cells, natural killer (NK) cells, γδ T cells, macrophages, dendritic cells, plasma cells, B cells, and plasmacytoid dendritic cell (pDC) clusters (Fig. 2A and B; Supplementary Figs. S3 and S4A; Supplementary Table S13).

First, we analyzed the 60,902 tumor cells, which were marked by the expression of neuroendocrine transcription factor ATOH1 and chromogranin (3). We examined expression of the viral oncogenes small T antigen (sTAg) or large T antigen (LTAg) by scRNA-seq and bulk RNA-seq. Twenty six of 43 patients had samples with multiple reads for both sTAg and LTAg and were thus considered VP. The remainder had no transcriptional evidence of MCPyV sTAg or LTAg expression, which we labeled VN.

Next, we performed differential gene and pathway expression analyses (Supplementary Table S14; Fig. 2C). In responders, we observed upregulation of multiple pathways. Similar to bulk RNA-seq, there was upregulation of genes associated with IFNγ (STAT1, MX1, IRF1) as well as IFNα (OAS1, IFI27, IRF9). Again, many IFNγ-associated genes are implicated in antigen presentation (i.e., HLA-I, HLA-II, ERAP2; Fig. 2D; Supplementary Table S15). Next, we ran SCENIC regulon analysis to identify significant transcription factor–associated transcriptional networks. Top regulons in responders highlighted type I/II interferon responses (IRF7, IRF9, STAT1, IRF1; Fig. 2E; Supplementary Table S16).

We calculated a module score to quantify IFNγ and IFN responses using previously described methods (10, 12). IFNα and IFNγ response module scores were highly correlated (Pearson correlation, R² = 0.928; P < 1e−10). Both signatures were associated with the infiltration of immune cells, suggestive of preexisting immunity (Fig. 2F). The source of type I interferons was not clear because it was not well captured by scRNA-seq. However, subsequent analyses showed that the IFNα module score in multiple cell types (tumor cells, T/NK cells, and stroma) was correlated with the presence of pDCs, which express high levels of type I interferons (Pearson correlation, R = 0.38; P = 0.011 in tumor cells; Supplementary Fig. S6A and S6B; ref. 13). Other possible candidate cell types showed lower correlations (classical dendritic cell (cDC) R = 0.045; macrophage R = −0.0010). The source of IFNγ appeared to be primarily CD8 and other NK/T-cells (Fig. 2G; Supplementary Fig. S7A). The enrichment of type I/II interferon gene sets in responders did not depend on the site of tissue sampling (LN vs. skin) or metastatic status (primary vs. metastasis; Supplementary Fig. S7B).

In nonresponders, we found upregulation of 21 genes and eight Hallmark pathways via Gene Set Enrichment Analysis (GSEA). These include MYC (P_adj = 3.05e−17), mTORC1 (P_adj = 1.68e−10), and E2F (P_adj = 7.62e−4). We hypothesized that the E2F signatures observed in bulk RNA-seq arose from tumor cells. To test this hypothesis, we first identified cycling tumor cells. These cells expressed markers of cell cycle/division, MKI67, TOP2A, CENPF, and E2F targets, such as CCNB2 and CDKN3, varying in proportion across samples from 0% to 5% of the total tumor cells. Because mTORC1 and MYC can be associated with proliferation, we confirmed that they retained significance even when controlled for cell cycle (Supplementary Fig. S8A and S8B; “Methods”).

Consistent with data from other tumor types (14), we observed that proliferating tumor cells was negatively correlated with both overall immune infiltration (Pearson correlation, R = −0.6; P = 1e−5) and the IFNγ response, which we posit is a response to IFNγ producing lymphocytes in the TME (Pearson correlation, R = −0.6, P = 4e−5; Fig. 2H). Finally, to test the effects on survival, we derived a tumor cycling signature and applied it to our bulk RNA-seq, which is significantly larger than our scRNA-seq dataset. Samples from the bulk RNA-seq dataset were scored based on the “tumor cycling” cell type signature derived from our single cell dataset with gene set variance analysis (GSVA) and split into the top or bottom quartile (see “Methods”). Patients with high levels of the tumor cycling signature had worse overall survival post ICB (Fig. 2I; Supplementary Table S17; P = 0.016, Cox proportional hazard).

To find additional drivers of primary resistance, we analyzed transcription factors upregulated in tumors from nonresponders. The most significantly upregulated gene in nonresponders was OTX2, a transcription factor involved in neural cell development. SCENIC analysis (15) confirmed that OTX2 target genes are also upregulated (Fig. 2E). The analysis identified additional transcription factors associated with neural stem cells (HOXA10, GATA2, and SIX1; Supplementary Table S16; refs. 16–18).

Collectively, tumor cell analysis identified genes and signatures predictive of the response to ICB that correlate with the numbers of immune cells and effector molecules produced by them in the tumor microenvironment. To assess this further, we analyzed the immune infiltrates.

We turned next to CD8 T cells because they express the highest levels of IFNγ in the tumor microenvironment (Fig. 2G; Supplementary Fig. S7B). Single-cell analysis identified 15 distinct clusters that were organized into two distinct metaclusters (Fig. 3A). The first, which we termed circulating T cells (Tcirc), included CD8 clusters C00_CD8, C01_CD8, C02_CD8, C09_CD8, C11_CD8, and C14_CD8. This cluster is enriched in T cells normally found in circulation, e.g., naïve, central memory, and terminally differentiated effector memory (Temra) T cells. Cells within this metacluster have upregulated genes or proteins associated with naïve and central memory cells. This includes markers of stemness, e.g., IL7R, TCF7, and LTB (19) markers of blood circulation, e.g., CCR7, SELL, S1PR1, and KLF2 (19). Changes in a subset of genes was confirmed at the protein level including CD127, CD62L, and CD45RA, which are upregulated in naïve, central memory, and Temra cells (Supplementary Table S18; Fig. 3B–D).

The second metacluster (comprising CD8 clusters C03_CD8, C04_CD8, C05_CD8, C06_CD8, C07_CD8, C08_CD8, C10_CD8, C12_CD8, and C13_CD8) showed the upregulation of genes associated with tissue-resident T (Trm) differentiation. These include adhesion genes (ITGAE, ITGA1, CD69), chemokine receptors normally found in skin-resident T cells (CXCR6), transcription factors (PRDM1 and RUNX3; ref.20), and skin Trm-specific metabolic genes (FABP5; Fig. 3D; Supplementary Table S19; ref. 21).

Trm cells displayed significantly increased expression of genes reflecting recent activation (41BB, CD38), as well as markers of chronic activation/exhaustion (LAG3, HAVCR2, PDCD1, ENTPD1, TOX; Fig. 3C; χ²; P < 2.2e−16). Based on the antibody-derived tags (ADT), the most common surface phenotype of Trm cells in MCCs was CD103⁺, CD137⁺, CD39⁺, CD279⁺, TIGIT⁺, CD8 T cells (Fig. 3B). With the addition of RNA markers, the majority of Trm cells also expressed LAG3 and HAVCR2 (Fig. 3C and D). Trm cells had over 2-fold increased expression of IFNγ compared to Tcirc cells, with clusters C03_CD8, C05_CD8, C07_CD8, and C13_CD8 showing the highest levels (Supplementary Table S20; Supplementary Fig. S9). They also showed the highest expression of cytotoxic enzymes including GZMB and PRF1 (Fig. 3B–D).

To elucidate the relationship between Trm and Tcirc metaclusters, we performed two analyses. First, diffusion mapping revealed two distinct developmental lineages (Fig. 3E), which largely mirrored Trm and Tcirc metaclusters. Second, we performed clonotype linkage analysis, which identifies CD8 clusters that have shared TCR clones and are thus developmentally linked. Hierarchical clustering highlighted the two Trm and Tcirc metaclusters (Supplementary Table S21; Fig. 3F). Analysis of individual TCRs underscored the lack of inter-mixing, as 88% of expanded T-cell clones were exclusively restricted to either the Trm or Tcirc metacluster (Fig. 3G), which is unlikely to occur by chance alone (Monte-Carlo simulation, P < 1e6; Supplementary Fig. S10A and S10B).

To confirm that Tcirc cells were derived from the blood, we compared the clonotypic overlap in patient-matched scRNA-seq from the blood. Consistent with our expectation, we found that Tcirc cells shared 4.95-fold more clonotypes with CD8 T cells found in the blood than Trm cells (P < 2.2e−16, χ²; Fig. 3H).

Pseudotime analysis of tissue CD8 T cells provided further insights into the molecular characteristics of these trajectories. We found that Tcirc cells differentiated from a naïve state to a circulating effector memory state, exhibiting decreased expression of naïve markers IL7R and CCR7, increased expression of GZMH, and a subsequent spike in GZMK expression, coinciding with JUND expression (Supplementary Fig. S11A and S11B). Trm cells acquired exhaustion markers during differentiation (TOX, PDCD1) and lost the expression of genes such as BTG1, which may function as a negative regulator of effector activity (Supplementary Fig. S11C and S11D; ref. 22).

We hypothesized that Trm CD8 T cells would be enriched for tumor antigen-specific cells because they showed the highest expression of both activation/exhaustion markers (proteins 41BB, CD38, PD-1, and the TOX gene) and cytotoxic enzymes (GZMA and GZMB; Fig. 3B–D). Moreover, the Trm lineage showed significantly higher levels of clonal expansion (Supplementary Table S22; Fig. 3I and J; Supplementary Fig. S12A and S12B).

To test this, we measured the NeoTCR8 score, a gene signature previously found to be enriched in tumor neoantigen-reactive T cells in multiple solid tumor types (23). The highest NeoTCR8 scores were found in the Trm clusters, with all but one of the Trm clusters ranked higher than any of the Tcirc clusters (Fig. 4A; Supplementary Fig. S9). Importantly, we found that this signature was similar across VP and VN-MCCs and was independent of the site of tumor sampling or stage (Supplementary Fig. S12C and S12D; Supplementary Fig. S13A–S13C).

VP MCCs showed larger TCR clone sizes among CD8 T cells compared to VN MCC. This may be because of a restricted tumor antigen space due to uniform expression of immunogenic viral proteins and few, if any, private somatic mutations/neoantigens (4). This feature facilitated functional validation. We transduced Jurkat-76 cells (αβ TCR knockout) with an NFAT-GFP reporter with four highly expanded TCRs (three Trm-restricted and one Tcirc-restricted) and co-cultured them with autologous lymphoblastoid cell lines (LCL) pulsed with peptide pools. These pools contain 116-amino acid peptides that tile along both the small T antigen and the large T antigen of MCPyV (Supplementary Table S23; ref. 24). All Trm-restricted TCRs demonstrated antigen-specific activation in response to a unique peptide antigen, confirming viral specificity (Fig. 4B and C; Supplementary Fig. S14A–S14D). The expanded Tcirc TCR did not show any specific activation (Supplementary Fig. S14E). Structural modeling supported the antigen specificities of the TCRs tested. When the TCR sequences were modeled with the experimentally validated peptides and the associated predicted MHC, 2/3 Trm TCRs had a TCRmodel2 score >0.8 out of 1, indicating high confidence of binding of the TCR to the functionally validated peptide-MHC complex (Fig. 4D; Supplementary Fig. S14F; ref. 25).

To link TCRs to public Merkel cell polyomavirus antigens within and across patients with VP-MCCs, we deployed GLIPH2, a computational pipeline that identifies biochemical similarities across complementarity-determining regions 3 (CDR3; ref. 26). We found 10 GLIPH2 motifs among 23 unique TCRs that were highly restricted to the VP cohort (>85%) and were highly expanded (>20 clones), suggesting that they respond to tumor-specific Merkel cell polyomavirus antigens (Fig. 4E). These TCRs were also highly enriched in the Trm cells (Supplementary Table S24). In contrast, we found fewer GLIPH motifs (n = 5) enriched in the VN cohort. Unlike for VP MCC, none of these motifs was expanded (i.e., clones of more than one cell) in multiple patients (Supplementary Table S24). This is consistent with the relative absence of public antigens in VN-MCCs where private mutations are created by UV-mediated mutagenesis.

To further explore the antigen specificity of MCC tumor-infiltrating lymphocytes (TIL), we compared our TCR sequences with known TCRs epitopes from three TCR databases that have compiled TCRs with known specificity to common viral pathogens that are not associated with MCC (27–29). We found 128 clonotypes that were predicted to bind to public viral antigens not expressed by Merkel cell carcinomas (including cytomegalovirus, influenza, and SARS-CoV-2 but not MCPyV). Consistent with the absence of these antigens in MCC tumor cells, these TCRs showed minimal expansion in the tumor microenvironment (Supplementary Table S25). When present, these clonotypes were significantly enriched in the Tcirc cell trajectory (χ² = 19.4; P = 0.00001; Supplementary Fig. S15A).

Finally, we assessed the clinical significance of Trm cells. First, the prevalence of Trm cells was correlated with response to ICB (P = 0.018; Supplementary Fig. S15B). This was corroborated by multiplex immunofluorescence (mIF) data performed on 44 samples from 38 patients (24 with assessable response), of which 17 had associated tumor scRNA-seq. No one marker set defines Trm cells; however, our single-cell data suggest that CD103 is expressed in 38.6%, a plurality of Trm cells. CD3⁺CD8⁺CD103⁺ skin-resident T cells were observed infiltrating the tumor parenchyma (Fig. 4F), and tumors from responders had significantly higher densities of CD3⁺CD8⁺CD103⁺ cells as compared to nonresponders in our mIF panel (Fig. 4G).

Strikingly, the ratio of Trm to Tcirc cells was also predictive of the clinical outcomes. The ratio of Trm:Tcirc was higher in responders, with all nonresponders having equal or more Tcirc cells compared to Trm cells. (Fig. 4H). This was confirmed in the bulk RNA-seq using GSVA on a Trm signature we derived from our single-cell dataset (see “Methods”). Patients with tumors with high expression of the Trm signature (the top quartile) had significantly improved median survival from the start of ICB compared to those with low Trm signatures (the bottom quartile; median survival not reached vs. 578 days, respectively, P = 0.038, Cox proportional hazard; Fig. 4I). Collectively, these data show that the presence of Trm cells is associated with both response to ICB and longer survival.

Next, we focused on samples from patients who responded to immunotherapy despite having a low IFNγ signature (Fig. 5A). This analysis identified non-conventional cytotoxic T cells as potential prognostic biomarkers. Single-cell analysis identified 17 clusters of γδ and NK cells (Fig. 5B). To unequivocally identify γδ T cells, we utilized TCR data obtained through our modification of ECCITE-seq (30), which enriches γδ TCR sequences. With one exception, most patients had low levels of NK cells (<2% of cells). In contrast, we found that γδ T cells were enriched (5.7%–25.0% of all cells) in 7/49 of tumor samples, including four samples from responders with a low IFNγ signature (Fig. 5A). We observed that tumors with low MHC-I expression were enriched for expression of TRDV1, which encodes Vδ1 TCR (Fig. 5C).

γδ T cells have three developmentally distinct cell types distinguished by the choice of Vδ segment during VDJ recombination (31). Using our γδ-specific TCR enrichment protocol, we found that Vδ1 cells in MCCs preferentially paired with Vγ3 and Vγ4 chains, and Vδ2 cells predominantly paired with the Vγ9 chain. We identified nine clusters of γδ cells: seven Vδ1 and two Vδ2, with few Vδ3 cells primarily found in the Vδ1 clusters (Fig. 5B–D).

The predominant γδ cells among the MCC TILs were Vδ1 cells, consistent with previous reports (Fig. 5D; ref. 32). Similar to CD8 Trm cells in MCC, Vδ1 and Vδ3 cells are tissue resident. In contrast, Vδ2 cells are the primary γδ cell in the circulation. Vδ1 and Vδ3 cells (but not Vδ2 cells) appeared to undergo antigen-dependent activation. Vδ1 cells expressed multiple markers of co-inhibition (LAG3, PDCD1, TIM3), activation (41BB, KLRC2), and effector functions (GZMA, GZMB, FASLG), as well as multiple inhibitory killer cell immunoglobulin-like receptors (KIR; KIR2DL1, KIR2DL3), and transcription factors associated with CD8 exhaustion (EOMES, BATF) or tissue residence (RUNX3). Many of these were significantly enriched compared to those in Vδ2 cells (Supplementary Tables S26–S28; Fig. 5E–G). Conversely, Vδ2 cells showed high expression of circulatory markers (CCR7, CD127, CD62L) along with higher expression of AP-1 family members, thus broadly matching the expression patterns observed in CD8 Tcirc cells (Fig. 3D).

Using CITE-seq antibody-based analysis, tumor-resident Vδ1 cells expressed less surface CD3 compared to circulating Vδ1 cells and/or Vδ2 cells in the blood or tumor, reflective of the receptor downregulation following activation (Supplementary Table S28; Fig. 5F). For Vδ1 cells, there were 397 unique clonotypes, of which 128 (32.2%) were expanded, suggesting higher clonal proliferation and clonal focusing than in Vδ2 cells (372 total clones, 85 expanded). The median clonal size was significantly higher among Vδ1 (17 vs. 3; P < 2.2e−16, Wilcoxon test; Supplementary Fig. S16A). Almost all the expanded TCRs from 128 clones had biochemically similar δ-chain CDR3s (see “Methods”), suggesting that they responded in a clonotype-specific manner to a public antigen, whereas the rest of the Vδ1 CDR3s were not expanded and were not similar, suggesting a diverse background of TCRs that may not recognize the same antigen (Fig. 5H).

Pseudotime analyses suggest that Vδ1 cells start from a naïve-like population (cluster C08_GD), which consisted of T-cell clones that did not proliferate (as determined by clonotype analysis) and enriched for the expression of markers associated with naïve γδ T cells, including CD27, CD28, CD45RA, CD127, CCR7, and TCF7 (Supplementary Fig. S16B; ref.33). The Vδ1 trajectory shifted primarily toward clusters C02_GD, C15_GD, and C10_GD, matching the clonal expansion observed previously (Fig. 5I; Supplementary Fig. S16C). These differentiated cells had higher levels of cell cycle markers (MKI67), exhaustion markers (ENTPD1, LAG3, HAVCR2, TIGIT, BATF), cytotoxic effectors (GZMH and GZMB), and cell migration proteins (CXCR6 and VCAM1; Supplementary Table S29; Fig. 5J). These data suggest that Vδ1 cells are CD8 T cell-like, in that they differentiate to a chronically activated phenotype in response to a CDR3-dependent antigen (Fig. 3D). Unfortunately, there were insufficient cells to perform pseudotime analysis of Vδ2 or Vδ3 cells.

To orthogonally validate the functional relevance of γδ T cells in MCCs, we performed mIF on seven patients (eight samples), including four with RNA-seq signatures suggestive of γδ T–cell infiltration. For a subset of MCCs, γδ cells were also detected in the tumor parenchyma expressing cytotoxic markers (Granzyme B; Fig. 5K). Interestingly, we found no statistically significant differences in Vδ1 infiltration between skin and LN samples in either our single-cell data or bulk RNA-seq data. All highly enriched Vδ1 samples were found in the skin, but the numbers of samples preclude statistical significance (Supplementary Fig. S13).

To resolve the TME spatially, we first used the GeoMx Digital Spatial Profiler, which identified 180 regions of interest (ROI) across 15 samples from 12 patients. These 12 patients included four complete responders to ICB after biopsy and eight nonresponders. The varied ROIs allowed for a comprehensive overview of the tissue landscape (Supplementary Figs. S17A, S17B, and S18A).

Differential expression analysis (DE-seq) was performed on these ROIs. Within the tumor region of complete responders (CR), we observed both T-cell genes (LAG3, GZMB, TRDC) and interferon response genes (GBP4), while in nonresponders (all with progressive disease), there was prominent expression of tumor markers (keratins, CD44, HOX genes) and proliferation markers (TYMS, MKI67; Supplementary Table S30; Supplementary Fig. S18B).

We detected a strong IFNγ response signature (GBP4, HLA-I, TAP, and proteasome genes) in the stroma of CR samples. Additionally, T-cell markers (CD27 and CD7), B-cell markers (CD79A and POU2AF1), and macrophage markers (C1QB, CD68, CD163) were identified, confirming the localization of these cells in the stroma. Conversely, in progressive disease samples, we identified tumor markers (keratins, PRAME) and proliferation markers (TYMS, MCM, E2F1; Supplementary Table S31; Supplementary Fig. S18B). We utilized the Trm and Tcirc gene signatures that we had already developed to spatially examine these specific T-cell lineages. We observed that while both Trm and Tcirc were upregulated in the stroma of responders, only Trm was upregulated in the tumors of responders. Nonresponders failed to show enrichment of Trm or Tcirc genes, regardless of the location (Fig. 6A).

Next, for single-cell resolution, we performed CosMx Spatial Molecular Imager analysis. We selected 122 fields of view (FOV) sized 0.9 × 0.7 mm across the tumor–stromal interface from 13 samples across 12 patients (six responders, four nonresponders). We segmented 614,209 cells (Supplementary Table S32; Supplementary Fig. S19A–S19D; Fig. 6B; “Methods”). Next, we deployed SpatialTIME, a permutation-based method that we developed to generate a degree of clustering (DoC) score to rank whether cell populations are more colocalized versus more dispersed than a distribution based on complete spatial randomness across all possible binary comparisons at specified neighborhood radii (34). We then used the presence of elevated DoC to identify key cellular immune circuits associated with ICB response. To do this, we performed unsupervised hierarchical clustering to identify four cell types (T cells, B cells, pDC, and DC) exhibiting high DoC across all FOVs and the entire range of radii r = 250 to 500 (90–180 µm diameter; Supplementary Table S33; Supplementary Fig. S20A–S20D; Fig. 6C). The highest DoC corresponding to the highest levels of colocalization (>4.5 × 10⁵) was seen across two circuits: T cell–B cell–DC and T cell–B cell–pDC. An example of this is seen in an FOV highlighting T cells, B cells, and pDC all colocalized at the tumor–stroma interface (Fig. 6D). We further sought to identify whether these distributions differed between responders and nonresponders. Among responders, B cells, DCs, and pDCs showed 61.9-, 1.8-, and 14.3-fold higher counts, respectively, with a r = 300 (108 µm diameter) in the neighborhood around T cells (P < 2.2e−16, χ²; Fig. 6E; Supplementary Fig. S20E). Conversely, T cells were preferentially colocalized with 1.7- and 13.5-fold more tumor and tumor cycling cells, respectively, in the same sized neighborhoods (r = 300; 108 µm diameter), suggesting a paucity of intact immune circuits in nonresponders (Fig. 6E; Supplementary Fig. S20E). We interpret this to suggest that the presence of a key immunological triad of B cells, T cells, and DCs or pDCs is correlated with ICB response.

We then leveraged the droplet-based scRNA-seq dataset to identify receptor-ligand interactions because of higher numbers of transcripts per gene, focusing on B cell and myeloid populations (including dendritic cell subsets; Supplementary Fig. S6A). Using CellChat (35), we found that B cells, DCs, and pDCs interacted with Trm T cells via multiple chemokine–chemokine receptor interactions (i.e., CXCL16 to CXCR6), costimulatory ligand–receptor interactions (i.e., 41BBL (TNFSF9) to 41BB (TNFRSF9), CD80/CD86 to CD28), and activating cytokines [i.e., IL15 to IL15 receptor, IFNA2/IFNW1 to IFNAR1/2 (IFNα receptor); Supplementary Table S34; Fig. 6F].

To examine the effects of ICB, we examined a cohort of seven responders and seven nonresponders from which we sampled matched tumor tissue pre- and post-therapy. Among responders, 1,560 genes were upregulated and 995 downregulated after therapy. (DE-Seq2, P_adj. < 0.1; Supplementary Table S35). In contrast, only two differentially expressed genes were identified among nonresponders (Fig. 6G).

Overall, responders showed upregulation (1.34–1.44 fold; P < 0.05) of multiple inflammatory signatures after treatment including our Vδ1 signature and an IFNG signature (10). There was a strong overlap in the pathways and the genes upregulated in pre-ICB responders and those upregulated in responders post-ICB (Fig. 6H; Supplementary Table S36; Supplementary Fig. S21A and S21B). These include markers of B cells (e.g., immunoglobulins), pDC/DC markers (e.g., PTGDS), T-cell markers (e.g., CD3, TCR), chemokines/cytokines (e.g., IFNG, IL15, CXCL16), CD8 Trm markers (e.g., CD69, PRF1, CD8B, CXCR6), and multiple interferon response genes (e.g., HLA, IRF1, IFI16) and gene signatures of Trm cells (but not Tcirc cells; 1.44-fold; P = 0.031; Wilcoxon; Supplementary Fig. S21B). These data collectively suggest that ICB increases the recruitment or expansion of immune cells that were present prior to ICB. Cell cycle and cell proliferation pathways were downregulated in responders post treatment. These included markers of tumor proliferation (FOXM1, TOP2A, MKI67, E2F7) and MCC-specific markers (e.g., ATOH1, NEFM), suggesting tumor clearance. In contrast, nonresponders show little to no change after ICB treatment (2%–9% change, P > 0.5) and were thus lower than temporally matched responders pre- or post-ICB (Supplementary Fig. S21).

Next, we examined a cohort of scRNA-seq samples only acquired from patients after immunotherapy. These 15 samples included five additional samples that were projected on the original single-cell dataset (see “Methods”). Consistent with our bulk RNA-seq data, responders (n = 5) had 11.1-fold lower levels of tumor cycling cells and 3.01-fold higher CD8 Trm cells (3.01-fold increase) following ICB (P < 0.05, Wilcoxon; Fig. 6I; Supplementary Fig. S21C).

Next, we examined antigen receptor sequences to study clonal dynamics in the TME after ICB in our bulk RNA-seq data using MIXCR (36). As suggested in our differential expression analysis, we observed 60- to 400-fold more total TCR or BCR clones in responders versus nonresponders (Wilcoxon; P < 0.05). Nonresponders had no expanded TCRs and BCRs prior to treatment, while 3/7 responders had expanded TCR clones and 5/7 responders had expanded BCR clones prior to treatment. After treatment, 1/7 (TCR) and 2/7 (BCR) nonresponder samples showed evidence of new, expanded clones. In contrast, all of the responders had expansion of newly recruited and preexisting TCR and BCR clones (Supplementary Fig. S22A–S22D, see “Methods”). Thus, both preexisting immunity and the recruitment of new, reactive T cells and B cells are implicated in response.

As MIXCR data can be biased by differences in TCR mRNA between cells, we used our matched single-cell data to more precisely assess CD8 T-cell dynamics. We examined two paired sets of ICB responder samples with pre- and post-therapy scRNA-seq from both the blood and tissue. One sample pair had matched bulk RNA-seq and one pair was new. After treatment, we saw the same results as our MIXCR data, where both samples had expansion of preexisting clones, but also recruitment of expanded clones (Supplementary Fig. S23A and S23B). Expanded cells from both patients were enriched for putative Trm-like signatures. They expressed CD103 and exhaustion markers [exhaustion score (37), PD-1, and LAG3]. Following anti-PD1 therapy, expanded cells had significantly lower levels of exhaustion markers, including a reduced exhaustion score and a 2.11 to 5.44-fold decrease in PD-1 and LAG3 expression, respectively (P < 1e6; Wilcoxon; Fig. 6J; Supplementary Fig. S23C), which may reflect increased functionality after ICB. Full cell type analysis of pre-post matched scRNA-seq from the blood (21 pairs) did not reveal any significant cell clusters correlated with response or resistance. Tumor-specific T cells comprised a maximum of 0.5% of the circulating T cells, but the number of cells per sample was not high enough to reliably capture a cluster of those cells (38).

In this study, we compiled the largest database to date of MCC tumor samples and matched blood with multimodal profiling of a total of 116 patients. The breadth and depth of our studies have allowed us to confirm findings across multiple orthogonal datasets, highlight the importance of tumor sampling (over blood), and identify multiple novel candidate biomarkers that may be clinically actionable. For instance, the negative impact of tumor cell proliferation on ICB response (Figs. 1F, G, 2H, and I) may be validated in CLIA settings using assays assessing Ki-67 expression. Moreover, next-generation sequencing may similarly enable prospective validation of IFNγ and/or IFNα signatures to predict responses to immunotherapy (Figs. 1D, E, 2D, and F; ref. 10).

Our data emphasize the role of tissue-resident lymphocytes in MCCs. Despite representing a minority (∼45%) of the total CD8 T cells, tissue-resident CD8 T cells appear to be the primary tumor-specific population responsible for tumor recognition and killing. Moreover, the prevalence of this tissue-resident population (but not the circulating T-cell population) is clinically important. In contrast to some other cancer types (39), their presence in the tumor prior to immunotherapy predicted responses (Fig. 4H and I).

Our study provides novel insights into γδ T cells in the TME. While others have reported significant infiltrates of γδ T cells (32), our modified ECCITE-seq allowed unique γδ TCR sequencing with protein markers, which enables linkage of VDJ, transcriptome, and proteomics. This analysis indicated that skin-resident Vγ3Vδ1 and Vγ4Vδ1 cells were the tumor-reactive population (and not the circulating Vγ9Vδ2 cells; Fig. 5E–J).

Vδ1 cells in MCCs expressed costimulatory and inhibitory receptors commonly expressed by innate NK cells including killer cell lectin-like receptor (KLR)/KIR family members (Fig. 5E and J). However, our analysis highlighted deep commonalities with CD8 Trm cells. Like CD8 Trm cells, they express acute activation markers (such as 41BB), undergo clonal expansion, and acquire markers associated with T-cell exhaustion (e.g., PD-1, LAG3, TIGIT). The pseudotime analyses paralleled those of CD8 Trm cells, wherein they both started in a TCF7+ cell state and ended in an exhaustion state and loss of BTG1, a transcriptional repressor of lymphocyte effector function (Fig. 5J; ref. 22). Taken together, these skin-resident Vδ1 cells appear to display a form of functional convergence with skin-resident CD8 Trm cells in driving antitumor immunity (Supplementary Fig. S24A and S24B).

Consistent with the recent reports, these Vδ1 cells are enriched in tumors with altered MHC-I expression (40). In contrast to previous reports, our TCR sequencing data suggest that the CDR3 region is critically important for Vδ1 cells (Figs. 4D and 5I; refs. 40, 41). Through single-cell γδ CDR3 sequencing, we demonstrated clonal expansion of Vδ1 cells that express biochemically similar CDR3s, highlighting the possibility of shared reactivity to a common (as yet unidentified) antigen (Fig. 5I).

Importantly, multiple spatial approaches suggest that T cells in the MCC TME do not function independently. Our data highlights the importance of cellular neighborhoods in response to immunotherapy (Fig. 6A–E).

Finally, we identified novel biomarkers, such as IL1, type I interferon and neural stem cell markers. Both IL1 and type I interferon have been described as having opposing protumor and antitumor roles in cancer (42–44). However, our data suggest that exogenous type I interferon and/or IL1 blockade may enhance immunotherapy efficacy. OTX2 and other neural stem cell markers may represent novel targets to improve response to ICB in MCC (Supplementary Fig. S25). MCCs may have unique, yet undiscovered, properties that highlight the biological importance of these molecules, specifically for this cancer type. Alternatively, this unique genomically homogeneous dataset (particularly VP-MCCs) may provide statistical power to make discoveries that are not possible for genetically heterogeneous cancers of other types. In other cancer types, these signals may be subsumed, for example, by the tumor mutation burden, which is required to make the tumor immunogenic. Thus, we believe that many of these results are applicable to other cancers, especially skin cancers.

A total of 116 patients with MCC treated at three institutions (Moffitt Cancer Center, Northwestern University Lurie Cancer Center, University of Washington) were included for analysis. The cohort was 79% male, median age 73 years (range, 18–93 years). Disease stage and immunocompromised status was reported in Table 1 and Supplementary Table S1. Inclusion criteria were histologically confirmed diagnoses of MCC made within the Department of Pathology at the corresponding institution and treatment with systemic immunotherapy. The clinical response to ICB was defined as complete response, partial response, stable disease, or progressive disease. However, for the purposes of all molecular analyses, characterization of clinical response to ICB was dichotomized to responder (complete and partial response) versus nonresponder (stable and progressive disease) status. These assessments were made on manual review of all electronic medical records on the basis of a combination of histopathologic evaluation (for neoadjuvant therapy), clinical assessment (clinical assessment of skin and/or lymphadenopathy), and serial radiographic assessment (computed tomography or positron emission tomography—computed tomography).

Complete responses were defined as complete radiographic and clinical response, or, where applicable, concordant pathological complete response. In the case of adjuvant immunotherapy, complete responses were noted in situations where surgery with curative intent had evaluable clinical or radiographic residual disease burden that was subsequently treated with ICB with clearance of disease. Any intermediate responses with residual disease burden were considered partial responses.

Progressive disease was defined as enlargement of overall burden of disease as noted by a combination of clinical and radiological assessment and/or the emergence of new metastases. In the case of adjuvant ICB administered following surgery with or without radiation, any recurrence or progression was considered progressive disease (n = 3). One patient on adjuvant ICB had stable disease with neither disease resolution nor progression. One mixed response was noted and considered an overall partial response in terms of overall tumor burden but the lesion used for molecular analysis showed no evidence of response to therapy and was considered as a progressive disease (MO_MCC_057). Death or date of last follow up was the clinical endpoint used for survival analysis.

Fresh tissue and blood specimens were collected from patients with MCC treated and managed at the H. Lee Moffitt Cancer Center Department of Cutaneous Oncology following written informed consent under the IRB-approved Total Cancer Care Protocol (MCC 14690, 94962) and performed in accordance with the Declaration of Helsinki. Collected specimens were processed and used under IRB-approved MCC Protocol 19191. Freshly collected tumor specimens were dissected for single-cell dissociation, formalin fixation with subsequent paraffin embedding, or flash frozen. Additional archival tissues at Moffitt were processed under waiver of consent (MCC 19191).

At Northwestern, fresh samples were provided by the Skin Biology & Diseases Resource-Based Center Translational & Experimental Skin Testing & Immune Tracing Core and were obtained through the IRB-approved Dermatology Tissue Acquisition and Repository (IRB # STU9443) or through the IRB-approved Dermatopathology Tissue Repository for Research (IRB # STU24696) and performed in accordance with the Declaration of Helsinki. Freshly collected tumor specimens were dissected for single-cell dissociation, formalin fixation with subsequent paraffin embedding, or flash frozen.

Blood was collected via phlebotomy following informed consent as above and PBMCs isolated under Ficoll gradient and cryopreserved in ∼1 to 2e6 cell aliquots.

Human skin tumor biopsies or portions of surgical resections were dissociated into single-cell suspensions using either the Miltenyi Biotec human tumor dissociation kit (enzymes H, R, A) or a mixture of RPMI, DNAse (0.1 mg/mL), collagenase IV (2 mg/mL), and hyaluronidase (0.1 mg/mL). Biopsies were gently minced into 1 to 2 mm pieces, transferred to an enzyme mix as above, and incubated at 37°C for 30 to 60 minutes at 800 rpm in a thermomixer or by gentle rocking. Following incubation, cells were filtered, washed, and cryopreserved.

Cryopreserved, digested MCC tumor samples were thawed and processed following a modified version of the 10× Demonstrated Protocol “Thawing Dissociated Tumor Cells for Single Cell RNA Sequencing” and BioLegend TotalSeq C protocol. The samples were thawed, labeled with TotalSeq C antibodies, and stained with 7-AAD viability staining solution and Apotracker Green (BioLegend). Cells were then subjected to fluorescence-activated cell sorting (FACS) to sort the CD45⁺ and CD45⁻ cell populations for target enrichment of 2/3 CD45⁺ cells. After sorting, the final cell preparation was resuspended in RPMI + 10% FBS + SUPERase•In RNase Inhibitor (Sigma), filtered, and counted to ensure an optimal concentration of 1e6 cells/mL. Finally, the cell suspension was processed for encapsulation and reverse transcription using the 10× 5′ v2 Chromium platform, according to the manufacturer’s instructions.

Cryopreserved PBMCs were thawed, labeled, and prepared for 10× droplet encapsulation following a modified protocol adapted from the 10× and BioLegend TotalSeq C protocols. Cells were blocked using Human TruStain FcX Fc blocking reagent (BioLegend) and stained with TotalSeq C and hashing antibodies (BioLegend). Live cells were sorted using FACS with a gate set for CD45⁺ and CD3⁺ cells. The sorted cells were combined with CD3⁻ cells to maintain a 1:3 ratio and cryopreserved in Bambacker preservation medium with SUPERase•In RNase Inhibitor (Sigma). Equal volumes of four to five hashed samples were combined to ensure that samples stained with the same hashing antibody were not mixed. The cells were slow-frozen at −80°C before being transferred to liquid nitrogen for long-term storage. Thawing and final sample preparation for encapsulation were performed immediately before loading onto the 10× 5′ v2 Chromium platform. Sequencing parameters were set for the 5′ cDNA, αβ VDJ, γδ VDJ, and ADT libraries according to the targeted cell number for scRNA-seq.

αβ/γδ TCR libraries were prepared following the 10× Genomics Single Cell V(D)J protocol. TCR γδ libraries were prepared following the 10× Genomics protocol to generate TCR αβ libraries with modifications. Specific primers were designed, including R1_hTRDC (5′AGCTTGACAGCATTGTACTTCC), R1_hTRGC (5′TGTGTCGTTAGTCTTCATGGTGTTCC), R2_hTRDC (5′TCCTTCACCAGACAAGCGAC), R2_hTRGC (5′GATCCCAGAATCGTGTTGCTC), and the forward primer (5′GATCTACACTCTTTCCCTACACGACGC), which were resuspended in 100 µmol/L EB buffer. The library preparation involved two PCR reactions: PCR1, using Human γ/δ mix 1 (R1_hTRDC + R1_hTRGC), and PCR2, using Human γ/δ mix 2 (R2_hTRDC + R2_hTRGC), substituting the primers as detailed. PCR1 was performed using 2 to 5 µL full-length cDNA, 50 µL 2× Kapa HIFI Master Mix, 5 µL cDNA additive, 2 µL 20 µmol/L forward primer, and 2 µL 20 µmol/L human γ/δ mix 1, with a reaction volume of 100 µL using water. The PCR was performed for 14 cycles. PCR2 was performed with 35 µL bead elution, 50 µL 2× Kapa HIFI Master Mix, 5 µL cDNA additive, 2 µL 20 µmol/L forward primer, and 2 µL 20 µmol/L human γ/δ mix 2, with a reaction volume of 100 µL using water. The PCR was performed for 12 cycles. Cleanups and PCR conditions were identical to the 10× protocol but γ/δ libraries required additional amplification cycles because of the rarity of γδ T cells compared to αβ T cells in most cell populations. Finally, TCR γ/δ-enriched libraries were processed according to the 10× Genomics Single Cell V(D)J protocol.

Formalin-fixed, paraffin-embedded (FFPE) tissue blocks were sectioned at a thickness of 5 μm. RNA was extracted from FFPE tissue sections using the RNeasy FFPE Kit (Qiagen) according to the manufacturer’s instructions. The extracted RNA was treated with an RNase-Free DNase Set (Qiagen) to remove any genomic DNA contamination. RNA quality and quantity were assessed using a Bioanalyzer (Agilent) and Qubit fluorometer (Thermo Fisher Scientific), respectively.

RNA sequencing libraries were prepared using the TruSeq Stranded mRNA LT Sample Prep Kit (Illumina), according to the manufacturer’s instructions. Briefly, 500 ng total RNA was used for poly(A) mRNA selection, fragmentation, and cDNA synthesis. The resulting cDNA was end-repaired, adenylated, and ligated using Illumina-specific adapters. PCR amplification was performed to enrich adapter-ligated fragments. The final library was purified and validated for size distribution using a Bioanalyzer (Agilent) and quantified by qPCR using a Kapa Library Quantification Kit (Kapa Biosystems). Libraries were pooled at equimolar concentrations and sequenced using an Illumina NovaSeq 6000.

FFPE tissue samples were immunostained using the AKOYA Biosciences OPAL TM 7-Color Automation IHC kit (Waltham) on a BOND RX autostainer (Leica Biosystems). The OPAL 7 color kit uses tyramide signal amplification (TSA) conjugated to individual fluorophores to detect various targets within the multiplex assay.

The sections were baked at 65°C for 1 hour then transferred to BOND RX (Leica Biosystems). All subsequent steps (deparaffinization and antigen retrieval) were performed using an automated OPAL IHC procedure (AKOYA). OPAL staining of each antigen was performed as follows: heat-induced epitope retrieval (HIER) was achieved with EDTA buffer (pH 9.0) for 20 minutes at 95°C before the slides were blocked with AKOYA blocking buffer for 10 minutes. The slides were then incubated with the primary antibody, GZMB (CST, D6E9W, RRID:AB_2799313, 1:75, dye 570) at RT for 30 minutes, followed by incubation with OPAL HRP polymer and one of the OPAL fluorophores during the final TSA step. Individual antibody complexes were stripped after each round of antigen detection in an optimized order and were determined independently for each panel. After the final stripping step, DAPI counterstain was applied to the multiplexed slide and removed from the BOND RX for coverslipping with ProLong Diamond Antifade Mountant (Thermo Fisher Scientific). All slides were imaged using the Vectra3 Automated Quantitative Pathology Imaging System.

Multilayer TIFF images were exported from InForm (AKOYA) and loaded into HALO (Indica Labs) for quantitative image analysis. A classifier was trained to identify areas of the tumor, stroma, or non-tissue regions. Pan-cytokeratin was used to train tumor regions because it is a masking marker for tumor cells. A classifier was created and tested on various images of the image set. The tissue was segmented into individual cells using the DAPI marker, which stains the cell nuclei. For each marker, a positivity threshold within the nucleus or cytoplasm was determined per marker based on published staining patterns and intensities for that specific antibody. After setting a positive fluorescent threshold for each staining marker, the entire image set was analyzed using the created algorithm. The generated data included positive cell counts for each fluorescent marker in the cytoplasm or nucleus and the percentage of cells positive for the marker. The data was sorted by classification to include the entire tissue, stroma, or tumor regions. Along with the summary output, per-cell analysis was exported to provide the marker status, classification, and fluorescence intensities of every individual cell within an image.

Local density analysis was performed by outlining the bulk tumor based on para-nuclear dot-like accentuation of panCK expression. Epidermal/follicular expression of panCK was excluded from analysis. Following the definition of the tumor-stromal interface, bands of 40 µm widths were defined as +500 microns (out of tumor into stroma) and −500 microns (into tumor interior), and the local density of cell counts according to marker expression was calculated per mm².

Tissue slices of 15 MCC tumor samples from 12 patients were placed on charged glass slides and stained with four immunofluorescent markers to assist in the selection of 12 ROIs per sample. The tissue was selected from a cohort of patients with MCC who received immunotherapy (four responders and eight nonresponders). FFPE human skin samples were baked for 1 hour at 65°C and processed on the Leica automation platform in three major steps:(i) slide baking and deparaffinization, (ii) antigen retrieval for 20 minutes at 100°C, (iii) 1.0 μg/mL Proteinase K treatment for 15 minutes. Slides were incubated overnight with GeoMx CTA or WTA assay probe cocktail. The following day, the slides were washed and subjected to morphological marker incubation before loading onto a GeoMx machine for fluorescence scanning and ROI selection. Barcodes were collected and the subsequent barcoding reading was completed on the Illumina NGS platform.

The markers used in the ROI selection were SYTO13 (DNA/nuclei), KRT18 (keratin/cancer cells), PTPRC (CD45/lymphocytes), and CD3E (T cells). Reagent: Syto13; Manufacturer: Thermo Fisher; ID: 121301310. Protein: PanCK; Manufacturer: Novus; Clone: AE1/AE3; ID: NBP2-29429, RRID:AB_3068002. Protein: CD45; Manufacturer: CST; Clone: D9M8I; ID: 13917BF, RRID:AB_2750898. Protein: CD3; Manufacturer: Origene; Clone: UMAB54; ID: UM500048, RRID:AB_2629062.

The classifications assigned to the ROIs were based on manual examination by the pathologist in our team, aided by DAPI and PanCK immunofluorescent stain, as presented in Supplementary Fig. S18A. Selection of the ROIs was done by the pathologist before RNA quantitation by GeoMx. Those ROIs in which the pathologist observed a predominance of epithelial cells were classified as tumor ROIs. Contrarily, those ROIs showing the absence of epithelial cells were classified as stroma. Variations of these ROIs were also included, in which tumor or stroma ROIs could also present infiltration of immune cells (CD45 or CD3 stain). ROIs of arbitrary shapes were selected to represent the tumor and stroma, including areas close to the tumor–stromal interface and central tumor areas corresponding to at least 150 cells, assigned to five categories with respect to their spatial location and level of immune infiltration:(i) “immune cold” stroma (paucity of immune cells), (ii) immune-excluded stroma (immune cells distant from tumor nodule), (iii) immune-infiltrated stroma (stroma adjacent to tumor infiltrated by immune cells), (iv) tumor close to immune cells (adjacent tumor areas to immune-infiltrated stroma), and (v) inner tumor (central areas of tumor nodules typically devoid of significant immune infiltrates; Supplementary Fig. S17A). The slides were profiled for gene expression using the Nanostring Cancer Transcriptome Atlas panel (CTA, nine samples/108 ROIs) and Whole Transcriptome Atlas panel (WTA, six samples/72 ROIs). The number of reads per ROI varied from 1.1 × 10⁶ to 142.9 × 10⁶ (Supplementary Fig. S17B).

The resulting .dcc files were processed using the GeoMx Analysis Suite v2.5. Probe counts were filtered using the default values and aggregated at the gene level. Given the technical effects associated with the panel used (WTA or CTA), we applied ComBat correction as implemented in the R package sva. We applied third quartile normalization of the counts, followed by log transformation. Normalized counts were used to detect differentially expressed genes between immune infiltration categories and patient therapy responses within each spatial niche (stroma/tumor). Differential gene expression was performed via pairwise t tests in the scran package, and P values were adjusted for multiple testing using false discovery rate (FDR).

To study the spatial patterns in gene expression of immune and tumor cells, we profiled six slides containing MCC tissue using the 960-gene Immune Oncology panel Nanostring’s Spatial Molecular Imaging (SMI) platform. Tissue samples were obtained from 12 patients, six of which showed a positive response to therapy, four were nonresponders and two were not accessible. FFPE tissue sections were prepared for CosMx SMI profiling as described by He and colleagues (45). Morphological markers used are as follows.

Protein: CD298; Manufacturer: Abcam; ID: ab167390, RRID:AB_2890241. Protein: B2M; Manufacturer: Abcam; ID: ab214769, RRID:AB_1266996; Protein: PanCK; Manufacturer: Novus; ID: NB2-33200AF532, RRID:AB_2924722. Protein: CD45; Manufacturer: Novus; ID: NBP2-34528AF594, RRID:AB_960384; Protein: CD3; Manufacturer: Origene; ID: UM500048, RRID:AB_2629062.

A total of 122 FOVs were profiled (2–16 per patient), which were chosen by inspecting adjacent hematoxylin and eosin–stained tissue slices (Supplementary Fig. S19A) and reconfirmed by imaging with morphological markers. The raw dataset consisted of 712,617 image-segmented cells, and transcripts falling within the cells were summarized at the gene level. Cells without any recorded transcripts were removed from the downstream analyses (815 cells). The SMI platform includes 19 negative probes for assessment to assist in the removal of background levels for transcript detection and quality control. We removed from the data set cells with more than 5% negative probe transcripts as well as cells with less than 50 expressed genes. To reduce the effect of segmentation artifacts (i.e., two or more cells merged into a single segment—“doublets”), we removed cells with more than 500 expressed genes or 1,250 transcript counts. The filtering procedure resulted in 614,209 cells being used in downstream analyses. We applied SCTransform normalization to the resulting counts. Corrected counts were input to InSituTypeML, a Bayes classifier to transfer the cell types identified in the tissue single-cell dataset to the SMI cells. To account for the background signal, InSituTypeML used the average negative probe count per cell. We also inferred Louvain clusters and tested for differentially expressed genes among clusters using Seurat to manually curate annotations generated by InSituTypeML. We confirmed that DE genes from the assigned CosMx clusters significantly overlapped with the DE genes from the higher resolution single-cell dataset (Supplementary Fig. S19D).

To quantify the degree of dispersion or colocalization of specific cell types within the spatial immune landscape in the MCC samples, we used spatialTIME (34) to measure the DoC across all possible binary combinations of cell subsets based on scRNA-seq data. The method implemented in spatialTIME computes the bivariate Ripley’s K statistic and assumes that cells occur at the same frequency across an FOV (i.e., complete spatial randomness) and is further normalized for area and cell numbers. Because FOVs often show non-uniform cellularity, spatialTIME compares the observed K to the exact K estimate based on a spatial random point process. This difference is the derived measure called the DoC.

Autologous LCLs were generated using a modified version of the Current Immunology protocol 7.22. PBMCs were thawed and washed with prewarmed complete RPMI medium. After counting and assessing viability, cells were centrifuged and resuspended in complete RPMI with 1.1 µg/mL cyclosporin A concentration of 2.2 × 10⁶ cells/mL. The cell suspension was incubated for 1 hour at 37°C before adding the thawed Epstein-Barr virus (EBV) viral supernatant (ATCC) to achieve a final cell concentration of 2 × 10⁶ cells/mL, with the virus constituting 10% of the final volume. The culture was monitored daily, and when the medium turned yellow, it was gently aspirated and replaced with fresh RPMI medium containing 1 µg/mL cyclosporin A. After 2 to 3 weeks, large cell clumps appeared and cyclosporin was no longer added to the culture. LCLs were maintained at high densities (>1 × 10⁶cells/mL) and thrived in slightly acidic (yellowish) media, with pipetting used to break up large cell clumps, as needed. After expansion, the LCLs were cryopreserved for further use.

The patient’s LCLs were thawed and treated with mitomycin C to inhibit proliferation, followed by washing and plating. Jurkat-76 cells previously transduced with a CD8-NFAT-GFP reporter construct (RRID:Addgene_153417; ref.46) were lentivirally transduced with transgenic TCRs (TWIST Biosciences) and selected in 1 μg/mL puromycin. A total of 116, 13-mer peptides were designed to tile along the small T and large T antigens of MCPyV with 8-mer overlap (Genscript; Supplementary Table S23; ref. 24). Peptides were pooled into 22 pools, such that each peptide was present in two pools. On day 2, the peptide pools were added to the LCLs and incubated for an hour at 37°C. Peptide-loaded cells were washed with cRPMI. Jurkat reporters were added to each LCL well, mixed, and centrifuged before incubation for 20 to 24 hours at 37°C. The next day, the remaining cells were stained with CD69 PE (BioLegend Cat# 310905 (also 310906), RRID:AB_314840), HLA-DR APC-Cy7 (BioLegend Cat# 307618 (also 307617), RRID:AB_493586), CD3 BV421 (BioLegend Cat# 344834 (also 344833), RRID:AB_2565675) antibodies for 30 minutes. The cells were washed, resuspended with 7AAD, and analyzed using a Fortessa HTS flow cytometer to measure Live, HLA-DR-, CD3⁺ CD69⁺ GFP + cells. Jurkat-76 cells were received as a gift from Mirjam H. M. Heemskerk on September 25, 2018 (47). They were authenticated as TCR-negative but no additional authentication was performed. Upon initial culture, the cells were periodically tested for mycoplasma with the MycoAlert kit (Lonza). Cultures were passaged up to five times before experiments.

The GRCh38.p13 primary assembly from Ensembl Release 98 and the GENCODE v32 primary assembly were modified as follows. CellRanger 6.1.0 was used to create a reference for scRNA-seq alignment. The MCPyV sequence was then added to the primary and final GTF files.

FASTQ files containing gene expression and feature barcodes (ADT) were aligned to the human hg38 and MCPyV genomes using CellRanger 6.1.0. Ambient RNA was removed using the CellBender pipeline (48). Processed read counts for each sample were imported as Seruat objects, after which TCR-seq data were added. We then filtered the data based on mitochondrial reads (25% cut-off) and minimum RNA genes (200 genes). Doublets were removed using DoubletFinder with the default pN of 0.25 and the optimal pK value calculated with the “find.pK” function. Cells were annotated with ProjectTILS and SingleR using Monaco and HPCA reference datasets. Data were merged into a single Seurat object. We performed standard Seurat (4.0) clustering, Harmony batch correction, and manual cell cluster annotation based on marker genes, including transcription factors, effector molecules, and previously annotated cell types. The data were subclustered by major cell types. Contaminating cells were filtered based on marker gene expression. Differential gene expression, gene enrichment, and GSVA were performed using the default Seurat pipeline, followed by comparison of cell type proportions, analysis of TCR clonotype expression, and TCR linkage analysis.

Protein-coding genes were subset upregulated or downregulated genes as indicated and input into the Enrichr pipeline using the enrichR package (3.0) and querying all available databases (9). The resulting enriched gene sets were grouped by database and then sorted by “combined score” which was labeled the Enrichr score.

The cells of interest were subset from the full scRNA-seq object. Cells were then grouped by sample into a pseudobulk matrix using Seurat’s “AggregateExpression” function. Next, this matrix was used to input into differential expression analysis using edgeR (3.36.0) using the “glmLRT” function.

Genes from differential expression analysis were ranked by signed P value. These genes were used as input for the fgsea function from the fgsea package (1.24). Reference genesets were derived from all human geneset in the msigdbr package (7.5.1). To control for cell cycle and tumor cycling cells, GSEA was run on two different pseudobulk DE results on tumor cells. The first result was without any tumor cycling cells and the second was controlling for cell cycling scoring as described in the pseudotime analysis and covariate DE-seq.

In the absence of specific hypothesis testing, sample size rationales were not calculated a priori. The major components of data were bulk RNA-seq and scRNA-seq data. These samples numbered 85 (63 for pre-ICB analysis) and 109, respectively. For bulk RNA-seq, primary conclusions were based upon signature based analyses and the power calculation for the RNA-seq experiment was conducted using package “ssizeRNA”. With an assumption of 10,000 tested genes (including the top 100 prognostic genes), a minimum fold change criterion of 2, and accounting for dispersion parameters ranging from 0.1 to 0.4, 63 samples (comprising 32 responders and 31 nonresponders) achieves an estimated probability of more than 80% (desired power) to reject the null hypothesis assuming equality in expression means between the two groups, with an associated FDR of 0.05. Additionally, we performed power calculations based on the gene set–based signature, specifically GSVA scores. Assuming a normal distribution of GSVA scores, our sample size of 63 provides a power of 90% to detect effect sizes (Cohen’s d) at 0.86 and a power >80% to detect effect sizes at 0.75, while assuming a type I error rate of 0.05. In our scRNA-seq experiments, where each sample comprises an average of 5,129 cells. Assuming that each subpopulation contains a minimum of 10 cells, our approach achieves a 95% probability of success to detect cell subpopulations with frequency as rare as 0.003.

Cells in the scRNA-seq dataset were assigned a module score using the “AddModuleScore” function in the Seurat package. This function was run on the combined Seurat object with all samples. In some cases, the cells were grouped by sample and/or cell type, and the average module score was reported. The following gene sets were used: Hallmark interferon gamma signature (MSIGDB), Hallmark interferon alpha signature (MSIGDB), exhaustion score from the cited paper (37), and NeoTCR8 score from the cited paper (23).

Following loom file conversion from the subset Seurat object, we executed the standard pySCENIC pipeline (v 0.12.1; ref. 15) in three major steps: gene regulatory network inference, context-dependent enrichments using cis-regulatory motifs, and activity scoring of regulons in individual cells using the AUCell algorithm. Subsequent to the SCENIC analysis, we integrated the results into our Seurat object. We used the FindAllMarkers function, ensuring even sample representation with equal number of cells, to find differentially expressed regulons between responders and nonresponders.

Cells were reclustered with cell cycle genes regressed using Seurat’s “CellCycleScoring” and “ScaleData” functions. Batch-corrected principal components were input into the DiffusionMap function of the Destiny package (3.14).

Cells were subset and reclustered with cell cycle genes regressed using Seurat’s “CellCycleScoring” and “ScaleData” functions. After cell cycle regression, pseudotime was determined using the default pipeline in monocle3 (1.0.0). The “root” node for pseudotime was determine by expression of naïve like markers including CD45RA, TCF7, and non-expression of exhausted or memory markers including PD-1, CD45RO, TOX. Genes were then linearly correlated with pseudotime using Pearson correlation. The next genes were nonlinearly correlated to pseudotime using the “fitGAM” function in the tradeSeq (1.8.0) package. The expression levels of significant genes were scaled across pseudotime and smoothed using the “smooth.spline” function with five degrees of freedom from the stats package (4.1.1). Genes were ordered by maximum value across pseudotime and then smoothed/scaled expression was plotted as a heatmap using the ComplexHeatmap package (2.13.2, RRID:SCR_017270).

All cells from the tissue scRNA-seq dataset were made into a pseudobulk average expression Seurat object, grouped by sample and cell type. This was input into the CellChat (1.4.0) pipeline using the human database accessed on April 29, 2023.

Five new post-ICB samples were harmonized using the Harmony package. All cells were mapped to the original Uniform Manifold Approximation and Projection (UMAP) and clusters from the original 49 sample tissue MCC scRNA-seq data using the Seurat “FindTransferAnchors” and “MapQuery” functions.

FASTQ files for both αβ/γδ VDJ sequencing were input into the Cellranger (3.0.2) vdj pipeline. Both contig annotations and clonotype annotations were added as metadata to the scRNA-seq Seurat objects for downstream analysis.

Identical clonotypes were matched between cell clusters and the likelihood ratio of clonotype matching was calculated for each cluster pair. P-values were determined using Fisher’s test and adjusted using the Benjamini–Hochberg procedure. Finally, log-likelihood ratios were subjected to unsupervised clustering and plotted using the Pheatmap (v1.0.12) package.

All CD8 T cells were grouped according to a unique TCR clonotype. Each cell was randomly assigned either Tcirc or Trm. The number of clonotypes that were 100% restricted to either Trm or Tcirc was counted. This simulation was repeated 1e6 times. The P value was determined by the proportion of simulations that had equal or greater numbers of restricted clonotypes than what was observed in the dataset.

TCR clonotype information was extracted from the CD8 T cells. Additionally, samples were HLA types using Polysolver, Optitype, ArcasHLA, and HLA-HD from bulk and single RNA-seq data. These HLA types and the extracted clonotype information were input into the GLIPH2 online portal (http://50.255.35.37:8080) with the following parameters: GLIPH2 algorithm, reference version 2.0, CD8 reference, all_aa_interchangeable as “YES.” The results were filtered using the Fisher_score <0.05, number_unique_CDR3 >number_subject, and Freq > number_unique_CDR3.

TCR clonotype information was extracted from CD8 positive T cells: the TRB CDR3 and TRBV genes. TCRs were considered a match if both the TRB CDR3 and TRBV genes matched a TCR found in the database. The VDJdb, McPAS-TCR, and TRAdb databases were accessed on January 2, 2023 (27–29).

The delta CDR3 sequences were extracted from the scTCR dataset. All CDR3s were aligned pairwise to each other using the “pairwiseAlignment” function from the Biostrings package (2.62.0) with the following settings: BLOSUM62 substitution matrix, gap opening = −10, gap extension = −1. The pairwise alignments were then clustered using unsupervised hierarchical clustering.

We implemented the MIXCR software (ref. 36; 4.1.0) to extract T/B cell receptor sequences from the RNA-seq data. The command used was mixcr to analyze RNA-seq-tcr-full-length, with the parameters: species hsa to specify the human species and rna to indicate RNA input data. We included the -keep-non-CDR3-alignments flag to retain alignments that did not include the CDR3 region, which provides a more comprehensive perspective on the TCR repertoire. The input data consisted of paired-end RNA-seq reads. Expanded clones were defined as clones with greater than 1.5 times the reads of the least frequent clone per sample and per receptor type.

Sequence reads were trimmed and filtered using Trim Galore! (v0.6.7, RRID:SCR_011847) to remove low-quality reads, adapter sequences, and bases with low Phred scores. The reads were then aligned to the reference human genome (GRCh38.p13 primary assembly) augmented with the Merkel cell polyomavirus transcriptome using GENCODE v32 transcript annotation with STAR (v2.7.10b, RRID:SCR_004463) in two-pass mode. The dataset was quality controlled to remove any significant outliers. Differential gene expression analysis was performed using DESeq2 (1.28.1; RRID:SCR_000154) or other appropriate software packages. The limma package (RRID:SCR_010943) was used to fit the differential expressions to a generalized linear model. Pathway analysis was conducted with GSVA and GSEA.

The patient metadata associated with the corresponding RNA-seq were associated with a DESeq object (RRID:SCR_000154). Raw read counts were filtered to remove low-expression genes, defined as having more than 20% of the total samples containing zero reads. A variance-stabilizing transformation was applied to the DESeq object and a PCA plot was visualized, showing clustering by virus status, tissue of origin, and patient sex. Differential expression was determined by comparing the immunotherapy response with the covariates virus status, tissue of origin, and patient sex. For pre-post ICB DE-seq, only patients were used as covariates with samples compared by time point.

The Wilcoxon markers of the cell type of interest versus all other cells were compared to the Wilcoxon markers of that cell type compared to the next most similar cell type. Trm cells were compared to Tcirc cells and tumor cycling cells were compared to all other tumor cells. These markers were intersected and filtered for a log₂ fold change >0.5 and a P value < 0.05. The genes are listed in Supplementary Table S17.

The cell-type signature genes were used as a gene set for single-sample GSVA. The corresponding GSVA scores were then classified into quartiles. The Kaplan–Meier curve for immunotherapy-specific survival compared to the top 25% of the cell-type score against the bottom 25% cell-type score among patients who had an evaluable immunotherapy response.

TCR-pMHC complex models were generated using TCRmodel2 (25) and AlphaFold (ref. 49; v.2.3.0) deep learning structural prediction methods. Structural prediction using TCRmodel2 was performed via its online server (https://tcrmodel.ibbr.umd.edu/), with structural refinement. AlphaFold 2.3.0 was downloaded from the AlphaFold Github repository (https://github.com/deepmind/alphafold) in February 2023 and run on a local computing cluster, generating 25 predictions per complex, with refinement performed on the top-ranked predictions. No template date cutoff was applied during TCRmodel2 or AlphaFold modeling. For the MO_MCC_095c1 TCR-peptide-MHC complex (Fig. 4D), models from TCRmodel2 and AlphaFold were pooled and ranked using the AlphaFold model confidence score, which is associated with TCR-pMHC model accuracy (25). The AlphaFold-generated model was ranked the highest and was selected for further inspection. For all other TCR-peptide-MHC complexes, models were generated by TCRmodel2 and ranked by model confidence score.

Bulk RNA-seq and scRNA-seq data are available at GEO accession GSE235093. Raw sequencing data are available at dbGAP accession phs003629 (https://www.ncbi.nlm.nih.gov/gap/sstr/report/phs003629.v1.p1).

Z.Z. Reinstein reports grants from NCI during the conduct of the study; in addition, Z.Z. Reinstein has a patent for Novel cytokines, chemokine and receptors to augment Merkel cell carcinoma immunotherapy response pending. Y. Zhang reports grants from Alpha Omega Alpha Carolyn L. Kuckein Fellowship and ISID LEO Foundation Travel Grant during the conduct of the study. S. Chandra reports other support from Bristol Myers Squibb and Regeneron during the conduct of the study; other support from Bristol Myers Squibb, Regeneron, Novartis, and Pfizer outside the submitted work. J.B. Cheng reports grants from LEO Foundation and Veterans Health Administration Office of Research and Development during the conduct of the study; grants from Sun Pharma, Regeneron, Pfizer, Bristol Myers Squibb, and Janssen outside the submitted work. C.H. Chung reports personal fees from Fulgent, GenMab, AVEO, Seagen, Regeneron, Exelixis, and Bicara during the conduct of the study. N.I. Khushalani reports grants and non-financial support from Bristol Myers-Squibb, grants, personal fees, and non-financial support from Regeneron and Replimune, grants and personal fees from Merck, personal fees from Immunocore, other support from Incyte, and AstraZeneca, personal fees and non-financial support from Iovance Biotherapeutics, grants and personal fees from Novartis, personal fees from Nektar, personal fees and non-financial support from Castle Biosciences, personal fees from Instil Bio, grants from GlaxoSmithKline, HUYA Bioscience International, Celgene, and Modulation Therapeutics, other support from Bellicum, Asensus Surgical, Amarin, and T-Knife Therapeutics outside the submitted work. A.A. Sarnaik reports personal fees from Guidepoint, Gerson Lehrman Group, second City LLC, MJH Holdings, Clinical Education Alliance, Bluprint Oncology Concepts, Boxer Capital, Rising Tide Foundation, Istari Onc, and Keyquest Health outside the submitted work; in addition, A.A. Sarnaik has a patent for 61/973,002 issued, licensed, and with royalties paid from Iovance, a patent for 61/955,970 issued, licensed, and with royalties paid from Iovance, a patent for 62/612,915 issued, a patent for 14/974,357 issued, a patent for 16/157,174 issued, and a patent for 17/279,327 issued. V.K. Sondak reports personal fees from Bristol Myers Squibb, Genesis, Iovance Biotherapeutics, Merck, Novartis, Regeneron, and Helix Biopharma, grants from Neogene Therapeutics, Skyline Dx, and Turnstone Biologics outside the submitted work. C. Vaske reports other support from NantOmics LLC during the conduct of the study. S. Kim reports grants from Bristol Myers Squibb during the conduct of the study; grants from AstraZeneca outside the submitted work. A.S. Brohl reports personal fees from Deciphera and Bayer, and grants from Merck outside the submitted work. B.L. Fridley reports grants from National Institutes of Health during the conduct of the study. K.Y. Tsai reports non-financial support from NantOmics, grants and other support from Moffitt Cancer Center Pinellas Partners and Moffitt Cancer Center Donald A. Adam Melanoma & Skin Cancer Center of Excellence, and grants from NIH during the conduct of the study; personal fees from NFlection Therapeutics, WorldCare Clinical, and Verrica Pharmaceuticals outside the submitted work; in addition, K.Y. Tsai and J. Choi has a patent for methods to augment immunotherapy response pending to Northwestern University, Moffitt Cancer Center. J. Choi reports grants from V Foundation and Leukemia and Lymphoma Society during the conduct of the study; other support from Moonlight Bio outside the submitted work; None. No disclosures were reported by the other authors.

Z.Z. Reinstein: Data curation, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. Y. Zhang: Data curation, formal analysis, visualization, writing–original draft, writing–review and editing. O.E. Ospina: Data curation, formal analysis, visualization, writing–original draft, writing–review and editing. M.D. Nichols: Data curation, visualization. V.A. Chu: Investigation, visualization. A. de Mingo Pulido: Investigation. K. Prieto: Investigation. J.V. Nguyen: Formal analysis, investigation, visualization. R. Yin: Formal analysis, writing–original draft. C. Moran Segura: Formal analysis, investigation. A. Usman: Resources. B. Sell: Investigation. S. Ng: Investigation. J.V. de la Iglesia: Investigation. S. Chandra: Resources. J.A. Sosman: Resources, writing–review and editing. R.J. Cho: Methodology. J.B. Cheng: Methodology. E. Ivanova: Methodology. S.B. Koralov: Methodology. R.J. C. Slebos: Investigation. C.H. Chung: Investigation. N.I. Khushalani: Resources. J.L. Messina: Resources. A.A. Sarnaik: Resources. J.S. Zager: Resources. V.K. Sondak: Resources. C. Vaske: Investigation. S. Kim: Resources. A.S. Brohl: Resources. X. Mi: Formal analysis. B.G. Pierce: Formal analysis, supervision. X. Wang: Formal analysis. B.L. Fridley: Formal analysis, supervision. K.Y. Tsai: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing. J. Choi: Conceptualization, resources, formal analysis, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

This research was made possible through the Division of Research Pathology, Total Cancer Care protocol and supported by the Advanced Analytical Digital Laboratory, Tissue Core, Flow Cytometry Core, and Molecular Genomics Core at the H. Lee Moffitt Cancer Center & Research Institute, an NCI-designated Comprehensive Cancer Center (P30-CA076292). Funding support was provided by the Donald A. Adam Melanoma & Skin Cancer Center of Excellence (K.Y. Tsai), the Barry S. Greene Fund (K.Y. Tsai, A.S. Brohl), and the Moffitt Distinguished Scholar Award (N.I. Khushalani). The authors acknowledge the Nanostring Technology Access Program for access to both the GeoMx and CosMx platforms to generate spatial transcriptomic data. We would like to thank the Northwestern Dermatology Clinical Trials Unit for their assistance with tissue acquisition. This work was supported by the Northwestern University Flow Cytometry Core Facility and the Cancer Center Support Grant (NCI CA060553). The research reported in this publication was supported by the Northwestern University Skin Biology & Diseases Resource-Based Center of the National Institutes of Health under the award number P30AR075049. Funding for this study was provided by the National Cancer Institute grant T32 CA009560 and F30CA278298 (Z.Z. Reinstein), National Cancer Institute grant T32 CA233399 (O.E. Ospina, B.L. Fridley), the 2019 AOA Carolyn E. Kuckein Fellowship (Y. Zhang), the V Foundation for Cancer Research grant T2021-019 (J. Choi), the Bakewell Foundation (J. Choi), National Institutes of Health grants DP2 OD024475-01 (J. Choi), the National Comprehensive Cancer Network YIA (J. Choi), R35 GM144083 (B. Pierce), and the Leukemia and Lymphoma Society grant 1377-21 (J. Choi).