BCC is the most common NMSC and numerous studies have demonstrated the existence of intra-tumoral heterogeneity (ITH) within BCC, encompassing tumor cells, immune cells, and stromal cells of BCC (49–52). For example, Huang et al. extensively depicted the dynamic cell heterogeneity among tumor epithelial cells, immune cells, myeloid cells, endothelial cells, and stromal cells of infiltrative BCC and revealed the inflammatory characteristics and CD4⁺ Treg-derived impaired immunity of BCC by scRNA-seq (49). Ganier and colleagues also conducted a comprehensive comparison of the cell composition differences between BCC healthy skin from different anatomical sites and found that there were significant expansions of RGS5 ⁺ pericytes and POSTN ⁺ cancer-associated fibroblasts (CAFs) in BCC (52), the latter was in accord with the result of Chen et al. (49). Another prominent finding was that their data suggest that malignant cells of BCC could originate from cells of the inner and outer hair follicle bulb, given that they analyzed keratinocytes in the interfollicular epidermis and in follicular structures using scRNA-seq (52). However, due to differences in classification criteria or the selection of canonical marker genes, the heterogeneity of some cells cannot be compared across studies, such as fibroblasts (49, 51). In such instances, the sharing and integration of scRNA-seq data become of even greater importance.

The complicated cell-to-cell communication within the TME plays a pivotal role in regulating tumor progression. Through integrating single-cell and spatial transcriptomics, Yerly et al. uncovered the epithelial collective migration characteristics of tumor cells and extracellular matrix (ECM) remodeling features of CAFs from the invasive niche of invasive BCC, and Activin A likely promotes tumor cell collective migration by activating fibroblasts to remodel the ECM (50). Furthermore, this large degree of ECM remodeling in BCC is likely driven by the expression of collagen- and metalloproteinase-coding genes (51). Additionally, they also found tumor cells responded to the sudden burst of WNT5A-drived CAFs-specific inflammatory signaling pathways by producing heat shock proteins (51). In addition to the crosstalk between tumor cells and CAFs, a strong crosstalk between ANGPT2 ⁺ lymphatic endothelial cells and leukocytes was observed, demonstrating a strong ability to promote leukocyte migration and activation of this lymphatic endothelial cell (49).

Advanced BCCs exhibit a high mutational burden and ITH, and typically occur drug resistance or relapse after Smoothened (SMO) inhibitors or ICI treatments (31, 53). Yao and colleagues employed scRNA-seq to determine the molecular characteristics of SMO inhibitor-resistant tumor cells, which showed high expression of LY6D, LYPD3, and TACSTD2 (54). Mechanically, AP-1 and TGF-ß cooperativity drive the nuclear myocardin-related transcription factor (nMRTF) resistance pathway, which amplifies non-canonical GLI1 activity in the LY6D ⁺ resistant tumor cells (54). Indeed, the plasticity of tumor cells represents a significant contributing factor to drug resistance. BCC transdifferentiation into squamous cell carcinoma is a distinct and important form of resistance for BCC, discarding the dependency on the HH pathway (55). Interestingly, the LY6D ⁺ SMO inhibitor-resistant tumor cells, extensively cover the morphological and transcriptional cell states of basal to squamous cell carcinoma transition (56). With the transition processing, the expressions of LYPD3, TACSTD2, LY6D, CD44, and MUC1 were increased, while the expression of GLI1 was decreased, indicating the occurrence of drug resistance (56, 57). As a supplement, the same research team also demonstrated the mechanism regarding LY6D^– tumor epithelial proliferation. A TREM2⁺ VCAM1⁺ macrophage population was identified and spatially and functionally defined by integrating the scRNA-seq, the co-detection by indexing (CODEX) multiplex system, and the cytometry by time of flight (CyTOF) approach, which was found to be adjacent to the highly proliferative LY6D^– lower tumor epithelium area and promoted LY6D^– tumor epithelial proliferation via secretion of the ligand oncostatin-M (OSM) (58). In addition to the resistance to the SMO inhibitor, Pich-Bavastro et al. also found that Activin A secreted by tumor cells would mediate the polarization of CAFs and macrophages, thereby facilitating the resistance of tumor cells to ICI by scRNA-seq (59).

Collectively, most studies focused on the infiltrative BCC and identified high heterogeneity within BCC, including tumor cells and immune cells, indicating the complexity of BCC ( Table 1 ). For this reason, with the increase in the number of cells, a large number of cell subpopulations have been identified, but their exact functions in promoting or inhibiting cancer require further investigation. The investigation of drug resistance mechanisms in BCC is a highly active area of research, with notable progress having been made in understanding the role of LY6D^+/– tumor cells in BCC resistance ( Table 1 ). However, the comparability of cell subpopulations across different studies needs to be further clarified, either by using identical marker genes, or the same classification criteria or by integrating and reanalyzing different datasets. Moreover, the multi-omics approach is gradually demonstrating its unique potential in deciphering the heterogeneity and drug resistance of BCC (52).

cSCC is the second most common NMSC after BCC. As the first scRNA-seq study in cSCC, Ji et al. extensively dissected the ITH of cSCC using 10 tumors and matched adjacent tissues, and identified a cluster of tumor-specific keratinocytes (TSKs) (60). The TSKs exhibited high expression of MMP10, PTHLH, LAMC2, and SLITRK, localizing to a fibrovascular niche at the leading edge of the tumor as their adjacent cells were enriched for CAFs and endothelial cells. Furthermore, they also exhibited an abundant epithelial-mesenchymal transition (EMT) signature which serves as a hub for intercellular communication with stromal and immune cells, contributing to cSCC progression, immunosuppression, and ITH. Compared with the 10× Genomics Chromium platform, the Smart-seq2 strategy shows the capacity to detect a greater number of genes within each cell. Based on this, Yan et al. obtained 350 single cells from six primary UV-induced cSCC and determined the gene expression and chromosomal copy number variation pattern in cSCC (61).

In addition to focusing on cSCC itself, the recurrence of cSCC also warrants attention. The treatment of recurrent cSCC is usually more challenging, therefore exploring the recurrence mechanisms of cSCC can help prevent its recurrence and facilitate its treatment. Recurrent cSCC trended to exhibit worse immunosuppressed microenvironment featured by T cell exclusion, CD8⁺ T cell exhaustion, and enrichment of pro-tumor SPP1 ⁺ tumor-associated macrophages (TAMs) (62). Interestingly, a cluster of TSKs identified by Ji et al. in the primary cSCC (60), was also studied in the recurrent cSCC. Although the total number of cell communication events was reduced in the recurrent cSCC (possibly due to the smaller cell number in the recurrent cSCC), TSKs still demonstrated a robust communication with inflammatory CAFs (iCAFs) and IL7R ⁺ CAFs which may be associated with cSCC recurrence (62). Actually, IL7R ⁺ CAFs may be a subtype of iCAFs, as they expressed high levels of IL6, CXCL1, CXCL8, and other cytokines (62). Those studies highlight the ITH of cSCC and also emphasize the pivotal role of TSKs in cSCC progression and recurrence, through extensive communication with different kinds of CAFs.

cSCC can progress from normal skin under persistent UVR, or from AK, a precancerous lesion, or the in situ carcinoma, Bowen’s disease (BD). In recent years, there has been a growing body of research investigating this progression spectrum of cSCC, attempting to prevent and intercept at the early stages of cSCC. Generally, there was no significant difference in the proportion of different kinds of keratinocytes (63), or other immune and stromal cells between AK and normal skin (64), indicating a small change during the progression from normal skin to AK. Even so, some AK-specific keratinocytes or genes were still identified including ALDH3A1, IGFBP2 from Zou et al. (63), and IL1R2, WFDC2 from Li et al. (64). Similarly, frequent cell communications between AK-specific keratinocytes and papillary fibroblasts were observed in AK, especially through the ANGPTL4–ITGA5 interaction, indicating the central role of the ANGPTL pathway in AK development (64). Additionally, the number of proliferative basal cells was found to be significantly increased in BD, compared with AK (63). During the progression of cSCC, malignant cells activate fibroblasts into tumor-promoting CAFs. Thus, Schütz et al. demonstrated the heterogeneity of CAFs between BD and cSCC (65) and found that iCAFs and myofibroblastic CAFs (myCAFs) were dominantly present in BD and cSCC, respectively, in contrast to their absence in AK. Further analyses showed that iCAFs may derive from pro-inflammatory FBs playing immunoregulatory function, while myCAFs may originate from mesenchymal FBs which are associated with ECM organization (65). However, although most studies emphasize the importance of different types of CAFs during cSCC progression and recurrence, It remains an unsolved mystery regarding fibroblasts transform into CAFs at which stage during the progression of cSCC, and through what mechanisms.

Immunosuppression is a key risk for cSCC. Frazzette et al. found that Immunocompromised organ transplant recipients exhibited reduced proportions of cytotoxic and naïve tumor-infiltrating lymphocytes through scRNA-seq, which may be the cause of increased cSCC risk (66). Immunotherapy has been widely used in the treatment of high-risk cSCC. Adoptive cytotoxic T cell transfer-based immunotherapy (ACT) has a good initial therapeutic effect on cSCC, but there is a high probability of recurrence. Using cSCC mouse models, Miao et al. found that tumor-initiating stem cells selectively express CD80, which inhibits T cell activity upon engagement with cytotoxic T lymphocyte antigen 4 (CTLA4), suggesting that blockade of these cells may repress cSCC recurrence (67). Additionally, using scRNA-seq and scTCR-seq, Yost et al. found that cSCC and BCC that respond to anti-PD-1 therapy are derived from their own characteristics, continuously recruiting new T cells, rather than reinvigorating pre-existing tumor-infiltrating T cells through scRNA-seq (68).

Taken together, similar to BCC, the heterogeneity was uncovered first in cSCC and highlighted the important role of TSKs in cSCC progression and recurrence ( Table 2 ). Additionally, developing novel immunotherapeutic strategies that target the immunosuppressive microenvironment and prevent recurrence will be essential to improve patient outcomes. In recent years, an increasing number of studies have focused on the progression from normal skin to cSCC. The transition from AK to cSCC is marked by the activation of fibroblasts into tumor-promoting CAFs. However, the transformation of fibroblasts into CAFs during cSCC progression remains unclear.

MCC is an aggressive skin cancer typically caused by the MCPyV and UVR. Similar to other NMSC, the prominent ITH was observed in MCC. Das et al. identified two distinct MCC tumor cell states: a “mesenchymal-like” state characterized by an inflammatory phenotype and a “well-differentiated neuroepithelial” state, which is associated with acquired resistance to ICIs (69). Additionally, the study demonstrated that manipulation of the intrinsic tumor cell into a “mesenchymal-like” state could lead to increased sensitivity to PI3K inhibitors, suggesting a potential non-genomic resistance mechanism (69). Furthermore, it seems that the heterogeneity of MCC cells was not due to their viral status, as one study reported the co-clustering of MCC cells from virus-positive and virus-negative tumor samples (70). Specifically, MCC cells can be classified into two subtypes according to the expression of neuroendocrine (NE) related genes, and YAP1 suppresses cell-cycle progression in the MCC cells with high expression of NE genes in part through indirect downregulation of MCPyV large T antigen expression (70). The similarities and differences in gene signature between the “mesenchymal-like” and “high NE” states of MCC cells should be further investigated. Another study also reported the plasticity of MCC cells and they found that the inhibition of MCPyV T antigen expression could revert MCC cells to neuron-like cells revealing the potential conversion of an aggressive cancer phenotype to a differentiated state (71).

The immune cells and fibroblast states have also been dissected by scRNA-seq. Gherardin et al. identified a cluster of γδ T cells in MCC, exhibiting proinflammatory phenotype, which may be important effector cells and may serve as prognostic biomarkers and potential therapeutic targets for MCC (72). Previous studies have characterized tumor-specific T cell responses to MCPyV-positive MCC, which have typically been CD8⁺ T cells (73), but little is known about the T cell response to MCPyV-negative MCC. Therefore, Church et al. identified and characterized tumor-specific Th1-skewed CD4⁺ T cells targeting multiple neoantigens in a MCPyV-negative MCC patient who achieved a profound and durable partial response to anti-PD-L1 therapy (74). Fan et al. demonstrates that MCC-derived exosome-shuttle miR-375 induces fibroblast polarization by inhibiting RBPJ and TP53 expression, leading to a pro-tumorigenic microenvironment (75).

Therapeutic resistance has been studied by many studies. MHC class I molecules (MHC-I), encoded by the human leukocyte antigen (HLA-A, -B, and -C) genes, are the primary molecules for presenting endogenous tumor antigens to CD8⁺ T cells. Genetic loss of a single HLA allele has been proven to be associated with ICI resistance in MCC (76), as it results in an inability to effectively present tumor antigens, preventing CD8⁺ T cells from recognizing and killing tumor cells. Intriguingly, Paulson et al. found that tumor cells evaded recognition by CD8⁺ T cells by selectively downregulating the expression of specific HLA genes, resulting in acquired adoptive cellular immunotherapy and ICI therapy resistance, through scRNA-seq in two metastatic MCC patients (77), while the histone deacetylase inhibitor domatinostat can reverse this effect by restoring HLA gene expression on MCC cells (78), and cause decreased cell viability and apoptosis through suppressing the transcription factor HES1 (79). Meanwhile, the tumor-associated macrophages (TAMs) were the predominant myeloid subset in MCC, exhibiting an immunosuppressive gene signature, and were also associated with ICI resistance (80). Furthermore, Pulliam et al. found that the frequency of tumor-specific CD8⁺ T cells in the blood was positively correlated with ICI response and exhibited a senescent phenotype, while tumor-specific CD8⁺ T cells in the tumor showed an exhausted phenotype (73); and it was pointed out that the resistance to anti-PD-L1 originates from the decrease of HLA-I on tumor cells (73), which was consistent with the result of Paulson et al. (77). This suggests that ICI resistance at treatment re-initiation was not due to a shortage of tumor-specific T cells in the blood, but rather to the downregulation of HLA-I on tumor cells.

Overall, the ITH of MCC including MCPyV-positive and MCPyV-negative has been investigated through scRNA-seq technology ( Table 3 ). At least two states of MCC tumor cells have been identified, and one of them is sensitive to therapeutic treatments. The reversion or conversion of different subtypes indicates the plasticity of MCC cells. Furthermore, the mechanism of ICI resistance should be paid more attention in the future.

EMPD is one rare type of NMSC, mainly localizing to the vulvar, perianal, scrotal, penile, and axillary regions. To date, only one scRNA-seq study has been performed in EMPD. Song et al. revealed significant epithelial and immune cell heterogeneity using tumor and peri-tumoral skin samples from one scrotal EMPD patient through scRNA-seq (81). Furthermore, this study found that MSI1, an RNA-binding protein, is overexpressed in the basal epithelial cells of human EMPD skin, and ectopic overexpression of Msi1 in the epidermal basal layer of mice results in a phenotype resembling human EMPD at histopathological, single-cell, and molecular levels (81), which supports the in situ transformation theory, i.e. keratinocyte-to-Paget-like cell conversion of disease pathogenesis. Mechanically, Msi1 activates the Her2-mTOR signaling pathway by directly suppressing Cmtm5 and Msi1 can also directly repress Pten to activate mTOR signaling. In conclusion, this study provides compelling evidence that the Msi1-mTOR signaling pathway plays a critical role in the pathogenesis of EMPD and that targeting mTOR signaling may offer a novel therapeutic approach for treating this disease. Although only one scRNA-seq study was carried out because of the rarity of EMPD, it is valuable to explore the recurrence mechanism of EMPD. Similarly, multi-omics at the single-cell level should be considered to uncover the pathogenies of EMPD.