Single-cell dissection of Merkel cell carcinoma heterogeneity unveils transcriptomic plasticity and therapeutic vulnerabilities

Abstract

Merkel cell carcinoma (MCC), a rare but aggressive skin cancer, remains a challenge in the era of precision medicine. Immune checkpoint inhibitors (ICIs), the only approved therapy for advanced MCC, are impeded by high primary and acquired resistance. Hence, we dissect transcriptomic heterogeneity at single-cell resolution in a panel of patient tumors, revealing phenotypic plasticity in a subset of treatment-naive MCC. The tumor cells in a "mesenchymal-like" state are endowed with an inflamed phenotype that portends a better ICI response. This observation is also validated in the largest whole transcriptomic dataset available from MCC patient tumors. In contrast, ICI-resistant tumors predominantly express neuroepithelial markers in a well-differentiated state with "immune-cold" landscape. Importantly, a subtle shift to "mesenchymal-like" state reverts copanlisib resistance in primary MCC cells, highlighting potential strategies in patient stratification for therapeutics to harness tumor cell plasticity, augment treatment efficacy, and avert resistance.

Keywords: Merkel cell carcinoma; cellular plasticity; copanlisib; epigenetics; immunotherapy resistance; neuroendocrine carcinoma; single cell RNA-seq; skin cancer; tumor cell state; tumor heterogeneity.

Graphical abstract

Figure 1

scRNA-seq identifies phenotypic plasticity in treatment-naive MCC (A) UMAP visualization of 15 distinct clusters (n = 46,027 [all cells]). (B and C) Distribution of cell types and heatmap of top 20 variable genes. MC, macrophage/monocytes; CAF, cancer-associated fibroblasts. (D and E) (D) UMAP visualization of five tumor cell clusters (n = 22,978 cells) and (E) heatmap of top 20 variable genes in each cluster. (F) Pearson correlation heatmap of pseudo-bulk data from tumor cells with unsupervised hierarchical clustering. (G) Heatmap of differentially expressed genes between these two groups (log₂ fold change >0.25) with unsupervised clustering. (H) Scaled dot plot of epithelial and neuroendocrine marker genes in MCC tumor cells. (I) Hallmark Pathway by GSEA of MCCT_G1 and MCCT_G2 tumor cells. (J) Scaled dot plot of EMT transcription factors (TFs) and mesenchymal marker genes in MCC tumor cells. The vertical dotted line demarcates MCCT_G1 from MCCT_G2 in (H) and (J).

Figure 2

scRNA-seq reveals a well-differentiated neuroepithelial state in ICI-R MCC associated with distinct pathways (A) UMAP clusters of tumor cells from nine treatment-naive and two ICI-R tumors (n = 35,796 cells), split based on treatment status. (B) Heatmap of top 100 variable genes in treatment-naive and ICI-R tumor cells (n = 35,796 cells). (C) Hallmark Pathway by GSEA of ICI-R and MCCT_G1 tumor cells. (D and E) Scaled dot plots of (D) epithelial and neuroendocrine markers, and (E) EMT transcription factors (TFs) and mesenchymal markers in MCCT_G1, MCCT_G2, and ICI-R tumor cells, and the ICI-R dataset as reported by Paulson et al. (UW). (F) Pearson correlation heatmap of pseudo-bulk data from MCCT_G1, MCCT_G2, and ICI-R tumor cells. (G) Volcano plot of differentially expressed genes between ICI-R and MCCT_G2 tumor cells.

Figure 3

scRNA-seq charts heterogeneity of tumor-associated immune cells in MCC (A) Left: UMAP visualization of CD45⁺ cells (n = 12,796 cells) with seven distinct cell populations, namely B cells, CD4 T cells, CD8 T cells, CD4 and CD8 double-negative T cells (DNTC), natural killer cells (NK), dendritic cells (DC), and macrophage/monocytes (MC), and their percentages of distribution (color-matched histogram). Right: split UMAPs depicting distribution of cell types in MCCT_G1 and MCCT_G2 tumors and color-matched histograms of relative abundance. (B) Heatmap of top 20 variable genes across cell types identified in (A). (C) Left: UMAP visualization of nine distinct T cell subtypes in CD4/CD8 T cell population (n = 3,169 cells) from the dataset in (A). Middle: split UMAPs depicting groupwise distribution and color-matched histograms of relative abundance. Right: scaled dot plot of marker genes for the T cell subtypes. (D) Left (bar graph): relative abundance of DNTCs in MCCT_G1 and MCCT_G2 tumors. Right (violin plots): expression of canonical naive, effector, exhaustion, and memory markers in MCCT_G1 and MCCT_G2 tumors. (E) Split UMAP visualization of seven CD45⁺ immune cell subtypes (n = 14,895 cells) and color-matched histograms depicting their relative abundance (bottom) in MCCT_G1, MCCT_G2, and ICI-R tumors. (F) Violin plots of exhaustion-associated immune checkpoint molecules and cytokines in CD3 T cells.

Figure 4

CellChat analysis reveals stronger interactions of tumor cells with other cell types in MCCT_G1 tumor microenvironment (A) Circle plots representing the number (left) and weight/strength (right) of interactions between the seven cell types in MCCT_G1 (top) and MCCT_G2 (bottom) tumors. Outgoing signals from a cell type are represented by color-matched chords, with chord thickness proportional to the number/strength of the corresponding signal. (B) Chord diagrams of outgoing signals from MCCT_G1 (left) and MCCT_G2 (right) tumor cells to other cell types in the tumor microenvironment. Outgoing signals are color matched to the origin (tumor cells), with each arc representing one pathway and arc length depicting the strength. (C) Comparative bubble plot of communication probability for the top ligand-receptor pairs of outgoing signals from tumor cells to other cell types in the tumor microenvironment. The top three enriched pathways are highlighted in (B) and in the red boxes in (C).

Figure 5

MCCT_G1 tumors with an inflamed phenotype correlate with better ICI response and survival-associated genes (A) Violin plots of the 18-gene IFN-γ gene signature (IFNG18S), survival-associated gene signature (SAG), and death-associated gene signature (DAG) in MCCT_G1, MCCT_G2, and ICI-R tumors. (B) Signature score of IFNG18S, MESI-19 (generated from MCC tumor cells), and SIG-14 (generated from entire cell populations) in the glioblastoma multiforme dataset with known ICI response (R, responder; NR, non-responder). Significance of enrichment between responder (R) and non-responder (NR) was calculated using a non-parametric Mann-Whitney test. Solid black and dotted lines within each violin represent the median and quartiles, respectively. (C) Classification of 102 MCC patient tumor RNA-seq samples based on enrichment of survival-associated gene signature (SAG) and cross-validated against death-associated gene signature (DAG). (D) Enrichment score for IFNG18S, MESI-19, and SIG-14 in SAG_high and SAG_low patients identified in (C). ns, not significant; ∗∗p < 0.005, ∗∗∗p < 0.0005, ∗∗∗∗p < 0.00005 as analyzed by one-way ANOVA for (A) and unpaired non-parametric Mann-Whitney test for (D).

Figure 6

Patient-derived MCC cell lines retain “mesenchymal-like” state and transcriptomic plasticity (A) Pearson correlation heatmap with unsupervised clustering based on sample distance among eight MCC cell lines. (B) Hallmark Pathway enrichment by GSEA of MCC_G1 and MCC_G2 cell lines. (C and D) Heatmaps with unsupervised hierarchical clustering depicting normalized expression across eight MCC cell lines of (C) selected epithelial and neuroendocrine genes and (D) selected EMT transcription factors (TFs) and mesenchymal marker genes.

Figure 7

HDAC inhibition in MKL-1 and MCC-16 cells shifts intrinsic cell state and reverts PI3K therapeutic resistance (A) Immunoblots showing dose-dependent HDAC inhibition in MKL-1 (left) and MCC-16 (right) cells upon domatinostat treatment. Ratios of H3K27-Ac to total H3 are presented as mean ± SD for each dose (n = 3). (B) Hallmark Pathway by GSEA in MKL-1 cells treated with domatinostat by bulk RNA-seq. (C) Relative mRNA expression of EMT transcription factors in MKL-1 and MCC-16 cells receiving either DMSO (vehicle control) or 10 μM domatinostat for 24 h. Data are presented as mean ± SD, and all samples were run in triplicate. ns, not significant; ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001 as analyzed by unpaired Student’s t test. (D) Comparison of dose-dependent drug response in MKL-1 and MCC-16 cells treated with domatinostat or copanlisib alone (dose range 30 nM to 30 μM), or a combination of copanlisib (dose range 30 nM to 30 μM) with fixed dose of 500 nM domatinostat. Data are presented as mean ± SD for each dose, n = 6 per dose, with half-maximal growth inhibitory concentration as analyzed by non-linear regression model using GraphPad Prism.