Unravelling the landscape of skin cancer through single-cell transcriptomics

Highlights • Human skin is a multifunctional and heterogeneous organ with a vast cellular pool.• Single-cell RNA-sequencing provides an unprecedented opportunity to investigate and understand the heterogeneous landscape of skin cancers.• Analysis of single-cell skin cancer transcriptomes opens new avenues of investigation, including the potential for precision-targeted therapy.


Introduction
Cellular diversity and heterogeneity underlie the skin's many functions, which include sensation, stabilization of body temperature, maintaining fluid and electrolyte balance, vitamin and hormone synthesis, protection from external insults, and immunologic defense [1][2][3][4]. The skin is divided into three broad compartments-the epidermis, the dermis, and hypodermis/subcutaneous layer. Each compartment contains specific anatomical structures and diverse cell populations that contribute to skin physiology and homeostasis [1][2][3][4][5][6]. The epidermis is arranged in stratified epithelial layers that establish the skin permeability barrier. While keratinocytes comprise >90% of the epidermal cell population, the epidermis also contains melanocytes, Langerhans cells, and Merkel cells [1][2][3][4][5][6]. By comparison, the dermis is normally paucicellular and characterized by an extracellular matrix of collagen and elastin fibers that provides structural support as well as a niche for sparsely distributed fibroblasts and a variety of immune cells [1][2][3][4]7]. Skin adnexa, blood and lymphatic vessels, and nerve bundles are also housed in the dermis [1][2][3][4]. The hypodermis is primarily composed of adipocytes organized into fat lobules and loose connective tissue; however, it also contains fibroblasts and macrophages (Fig. 1). Crosstalk between different cell populations, such as structural epithelial cells and immune cells, is increasingly appreciated to play an important role in both epidermal homeostasis as well as disease states [1,3,4].
Our aging population, rising global temperatures, and challenges to changing sun-related attitudes and behaviors all contribute to the rising incidence of skin cancer. Keratinocyte carcinomas, namely, basal cell carcinomas (BCC) and cutaneous squamous cell carcinomas (cSCC) [8], account for the majority of these cases and are the most common cancers in the United States [9,10]. While less common, melanoma and Merkel cell carcinoma (MCC) also demonstrate growing incidence rates. Histopathologically, these tumors are characterized by nests or sheets of neoplastic cells within the epidermis and/or involving the dermis as well immune cell infiltrates that reflect the host immune response. Local treatment of small and/or early skin cancers by surgical excision or radiation therapy is often curative; however, effective treatments for advanced or metastatic tumors have met with more measured success. The heterogeneity of tumor cells and cells comprising the tumor microenvironment (TME) often underlies therapeutic resistance, suggesting that a deeper understanding of these cell populations and their interactions might facilitate the development of more effective targeted treatments. Single-cell RNA-sequencing (scRNA-seq) now provides an unprecedented opportunity to study malignant cell heterogeneity and TME complexity at single-cell resolution. A number of recent review articles details how scRNA-seq has recently been applied to study skin homeostasis as well as cutaneous inflammatory diseases and cancers [8,11,12]. Here we discuss how our understanding of skin cancer biology has advanced using bulk transcriptomics as well as single-cell approaches and highlight different methods of scRNA-seq analysis.

Recent applications of RNA-seq in skin cancer
When compared to hybridization-based microarrays, nextgeneration RNA sequencing (RNA-seq) possesses several strengths, including greater dynamic range and the capacity for de novo transcriptome assembly. In the skin cancer setting, whole transcriptome profiling by RNA-seq has been used to identify novel genes, such as noncoding RNAs and lowly expressed mRNA isoforms/transcripts. More recently, RNA-seq-based approaches have also been used to predict treatment response, nominate new therapeutic targets, and advance our understanding of disease progression as well as pathogenesis. For example, a recent analysis of melanocytic nevi and primary melanomas by Kunz et al. using RNA-seq revealed separate transcriptome signatures, with enrichment of genes driving BRAF/MEK inhibitor resistance in NRAS-mutant nevi and primary melanomas, while genes associated with resistance to PD-1 inhibition were expressed in lesions with wild type NRAS [13]. This work defined two different trajectories for melanoma development, suggesting that treatment resistance is determined early in melanomagenesis [13]. Additional work by Svedman et al. using targeted RNA sequencing (Ion AmpliSeq) to characterize melanomas following immune checkpoint inhibition similarly demonstrated enrichment of key pathways such as DNA replication, chromatin remodeling, and cell cycle, which may predict long-term response to therapy [14]. RNA-seq has also been used to highlight novel targets and signaling pathways in BCC. Wan et al. analyzed the transcriptomes of BCC with high tumor purity compared to patient-matched normal skin and identified enrichment of C2H2-type Zinc Finger genes, including Gli transcription factors as well as many others not functionally characterized in this malignancy [15]. This study, along with a larger pooled analysis of publicly available RNA-seq data by Litvinov et al. confirmed known dysregulated pathways in BCC and also suggested new therapeutic targets in the Wnt/β-catenin as well as IL-17 signaling pathways [15,16]. Other groups have used RNA-seq to elucidate the molecular underpinnings of skin malignancies. Using human specimens that represent increasingly dysplastic states on the cSCC spectrum as well as the equivalent tissues in SKH-1E cSCC-prone mice, Chitsazzadeh et al. distilled four TFs -ETS2, SP1, FOXF2 and AP1 -that appear to regulate the continuum of cSCC development [17]. Similar results were obtained by Das Mahapatra et al., who compared cSCC to unmatched normal skin controls and also reported enrichment of genes regulating immunological pathways, highlighting the importance of immunotherapy in this malignancy [18]. In MCC, which is divided into virus-positive and virus-negative categories based on the presence or absence of Merkel cell polyomavirus (McPyV) transcripts, Starrett et al. used RNA-seq to compare virus-positive to virus-negative MCC transcriptomes and demonstrated that genes regulating the cell cycle were increased in the former while DNA repair genes were downregulated in the latter [19]. These findings demonstrated a role for McPyV in controlling the MCC transcriptome and provide an explanation for the high mutation burden observed in virus-negative tumors [19].
As demonstrated by the aforementioned studies, bulk RNA-seq-based analyses are a valuable mainstay of investigative skin cancer biology. There are many advantages to using this approach, such as detection of mRNA isoforms and RNA species, high read capture capacity which enables detection of low abundance genes, low signal-to-noise ratio, and cost effectiveness. However, when applied to cancer, a few important limitations should be considered. First, details regarding the cellular context of the observed mRNA changes are not readily provided by bulk RNA-seq, although methods now exist to deconvolute this data into known cell populations using reference signatures. Second, bulk sequencing measures the average expression of a pooled population of cells, raising the possibility that transcriptomic changes of rare populations may be masked. This is true of many tissues including skin, where heterogeneity is present not only at the level of diverse cell types and states (Fig. 1), but also evolves from disease progression and treatment. The emergence of therapeutic resistance following cancer chemotherapy, for example, is a setting in which bulk sequencing can fall short due to its inability to distinguish between drug-sensitive and resistant clones. Additionally, the ability to explore RNA trajectory, cellspecific multi-omics or cell-cell communication is limited with bulk RNA-seq. Thus, these and similar questions are increasingly being addressed by single-cell approaches that provide a more granular view of cell types and states, potentially revealing new targetable transcriptomic changes.
The two widely used methods for single cell sequencing are 10X genomics and Smart-seq2. Smart-seq2 has been shown to detect low abundance transcripts and isoforms, while 10X genomics is better at detecting rare cell types [20]. Using these methods and scRNA-seq data analysis approaches, researchers have highlighted cellular/spatial context of skin cancer transcriptome and deciphered novel cell-to-cell communications.

Cellular context
A key goal of scRNA-seq is to characterize the cell types present within a sample, as this allows researchers to understand cell-type specific gene expression. An important first step when analyzing scRNA-seq data is to consider batch effects, or gene expression patterns caused by non-biological factors that can cause erroneous interpretation of cell types. Removing this technical noise for proper biological interpretation is referred to as integration. Methods like Seurat CCA, Harmony, scVI, scanorama, and scMC have been developed to integrate cells from multiple technical conditions [21][22][23][24][25]. Although integration methods help remove technical noise, when applied incorrectly they can overcorrect and remove biological variations.
Many approaches have been developed to identify cell types. Manual annotation is a popular technique combining dimensionality reduction and unsupervised clustering followed by manual inspection of marker genes. This approach works best for smaller studies or when most of the cells being analyzed are expected to have few known marker genes. For larger studies of samples with previously profiled cell types, automated classification methods can greatly speed up annotation. Some of these approaches, such as SingleR, CellAssign, and Garnett rely on transcriptomic reference profiles or gene signatures to label cell types [26][27][28][29][30]. By contrast other methods, including Azimuth, ProjecTILs, and scArches projects new cells to a reference atlas of previously annotated cell types [27,31,32]. Large consortium efforts such as Human Cell Atlas are currently cataloging major cell types and generating reference atlases that can be used to map and classify data from new studies [33].
After identifying broad cell types, further clustering is often performed on cells within a specific lineage (e.g., tumor cells, fibroblasts, Tcells) to identify subpopulations of cells with transcriptionally distinct states. This analysis is usually exploratory in nature and thus performed manually unless prior studies have already determined marker genes for distinct cell states. In the cancer setting, identification of cell types and subclusters is often used to characterize intra-and intertumoral heterogeneity and can also be applied to immune and stromal cells to study the transcriptional heterogeneity of the TME.
For example, in the first scRNA-seq study of skin cancer, Tirosh et al. analyzed 4645 cells isolated from 19 patients with melanoma and showed that across all tumors, malignant cells from the same melanoma contained two distinct transcriptional states not distinguishable by bulk RNA-seq that were characterized by either high expression of MITF or AXL. Cells with high AXL expression exhibited drug resistance to RAF/ MEK inhibition, suggesting they drive tumor recurrence [34]. This work also highlighted differential T-cell activation, expansion, and clonal exhaustion programs across different patients [34]. By sequencing single-cells from short-term culture of three melanomas with differing BRAF/NRAS genotypes, Gerber et al. derived gene expression modules from different cell subpopulations that included those with high MITF or AXL and also revealed targetable upregulation of the CDK4 and CDK2 cell cycle-dependent kinases responsive to palbociclib [35]. Yao and colleagues, identified similar patterns of intratumoral heterogeneity in BCC by utilizing scRNA-seq to characterize tumor cells based on expression of nuclear myocardin-related transcription factor (nMRTF), which correlates with Smoothened (SMO) inhibitor resistance and is targetable using AP-1 inhibitors [36]. In cSCC, Ji et al. also identified four different keratinocyte subpopulations, including tumor-specific keratinocytes (TSKs). TSKs are observed at the leading edge of the tumor and are thus positioned to act as a hub in cSCC tumor-stromal interactions [37]. Tumor subpopulation clustering was also used by Paulson et al. to elucidate mechanisms of treatment relapse in McPyV+ MCC by demonstrating tumor cells adapt to the immunological pressure created by autologous McPyV-specific CD8 + T cells and immune checkpoint inhibitors by suppressing the HLA specific to the targeted McPyV epitope [38].
When applied to cancer-associated immune and stromal cells, the identification of cell types and subclusters has broadened our understanding of the transcriptional heterogeneity of the TME. Several studies have examined the landscape of T-cells in skin cancers, finding a diverse range of cells including cytotoxic, effector, exhausted, regulatory, and helper T-cells. Yost et al. coupled scRNA and T-cell receptor sequencing (scTCR-seq) to study advanced BCC following PD-1 blockade and showed that novel T-cell clones infiltrate the TME of advanced BCC upon treatment, replacing their pre-existing exhausted counterparts rather than reinvigorating them [39]. Deng et al. performed a meta-analysis of 59 melanomas and catalogued the transcriptional states of CD8 + T-cells. They defined seven cytotoxic, exhausted, and naive/memory subpopulations, including a subset of exhausted T-cells associated with poor prognosis and characterized by high expression of PMEL, TYRP1, and EDNRB [40]. Frazzette et al. recently used scTCR-seq and gene expression data to compare the T-cell landscape in cSCC from skin cancer-prone immunocompromised organ transplant recipients (OTR) to that of immunocompetent individuals [41]. This effort highlighted several differences between the tumor-infiltrating lymphocytes (TILs) of these two groups, with a reduction in cytotoxic T-cell number, TCR clonotypes, and clonal expansion observed in OTR [41]. Other studies, such as that by Davidson and colleagues, have focused on characterizing cancer-associated fibroblasts (CAFs) [42]. Using the murine B16 melanoma model, these authors defined three functionally distinct CAF subpopulations that are positioned to play roles in immune crosstalk, extracellular matrix remodeling, and cytoskeletal reorganization [42].

Cell-cell communication
The availability of expression profiles from multiple cell types within a single sample allows detection of intercellular signaling via ligandreceptor expression analysis [43,44]. Several computational methods have been developed, including CellPhoneDB and CellChat [45,46]. These tools utilize databases with prior knowledge of ligand-receptor pairs or complexes to infer cell-cell communication based on expression of ligands and receptors from annotated cell groups (Fig. 2). This approach has proved to be especially useful for investigating interactions between tumor cells and their microenvironment as well as immune-stromal crosstalk in skin cancers (Fig. 2).
Within the TME, CAFs are well known to modulate tumor growth and progression. In their work identifying the TSK keratinocyte subcluster in cSCC, Ji and colleagues used NicheNet [37] to highlight TSK signaling to CAFs through ligand-receptor pairing of MMP9-LRP1 and TNC-SDC1 37 . Similarly, work by Davidson et al. used CellPhoneDB to identify a subset of CAFs in the murine melanoma TME that express the immunomodulatory C3 protein and are in close contact with C3aR-positive myeloid cells; disrupting this interaction reduced tumor growth and resulted in fewer myeloid cells that suppress CD8 + T-cells [42]. Deng et al. also applied CellPhoneDB to their re-analysis of publicly available single-cell melanoma transcriptomes to show that ligand-receptor interactions involving exhausted T-cells, such as CCL5-CCR5 and CD74-MIF, are central to communication between different CD8 + subpopulations [40]. Guerrero-Juarez et al. used CellChat to investigate tumor-stromal interactions in BCC and identified that the nearby fibroblasts act as an inflammatory hub regulating BCC growth through induction of heat shock proteins [47]. Knowledge of ligand-receptor pairing was also essential to work performed by Miao and colleagues, who performed adoptive cytotoxic T-cell transfer (ACT)-based immunotherapy in a mouse model of cSCC to investigate mechanisms of immune evasion [48]. Using scRNA-seq and lineage tracing, these authors identified a population of CD80 + tumor-initiating stem cells that are refractory to ACT and engage with CTLA4-expressing cytotoxic T-cells, reducing their activity and promoting tumor relapse [48].

Copy number alterations
Genetic variation is traditionally assayed using exome or genome sequencing, which provides data on mutations, structural variants, and genomic copy number. When applied to cancer, these methodologies allow researchers to query the clonal structure of tumors, estimate tumor purity and ploidy, and distinguish neoplastic from healthy samples. Single cell DNA-sequencing (scDNA-seq) enables elucidation of genetic heterogeneity at higher resolution than bulk exome/genome approaches; however, it is not used as widely as scRNA-seq due to its high cost, lack of commercial kits, and data quality issues [49]. Method developers have attempted to circumnavigate these issues by estimating single-cell copy number profiles from scRNA-seq data with tools such as inferCNV, honeyBadger, copyKat, and CaSpER [50][51][52]. These approaches have been used to distinguish neoplastic from normal cells and understand clonal lineages within tumors. For example, in their scRNA-seq analysis of melanoma, Tirosh et al. classified melanocyte-like cells as malignant, non-malignant, or intermediate based on a copy number alteration score [34]. Only malignant cells were included in subsequent analyses of tumor gene expression programs, thus avoiding contamination from normal melanocytes that can occur when analyzing bulk RNA-seq data [34]. Yost and colleagues demonstrated concordance between copy number variants inferred from single-cell transcriptomes and those called using whole exome sequencing data in BCC [39].

Gene trajectory
Many biological processes such as differentiation create cells that exist along a continuum rather than distinct transcriptional states. Trajectory inference methods have been developed to infer the relationship between cells across a continuous process [53,54]. RNA velocity analysis is a related technique that uses the proportion of unspliced and spliced mRNA molecules within a cell to provide temporal ordering along the trajectory [55] (Fig. 3). When applied to skin cancer, trajectory analysis has also been used to study tumor evolution. Wouters et al. employed Monocle-2 to characterize phenotype switching in melanoma [56]. Their findings demonstrate that melanoma cells transition from a melanocytic transcriptional state to a mesenchymal-like one through a stable intermediate state characterized by unique chromatin features and transcriptionally regulated by SOX6, NFATC2, EGR3, ELF1 and ETV4 56 . Another trajectory study by Su and colleagues found that BRAF V600E melanomas exhibit a bifurcated developmental trajectory in response to BRAF inhibition that is dependent on MITF expression [57]. Trajectory analyses have also provided new insight into the TME response to treatment. Using scRNA and T-cell receptor sequencing, Yost et al. showed that novel T-cell clones infiltrate the TME of advanced BCC following PD-1 blockade, replacing their pre-existing exhausted counterparts. The authors used trajectory analysis to analyze the evolution of T-cells in response to anti-PD-1 therapy and discovered they follow a bifurcated trajectory, becoming either terminally activated or exhausted, with the latter state marked by increased expression of PDCD1 and HAVCR2 39 . Deng and colleagues observed similar findings in melanoma, with their trajectory analysis suggesting that CD8 + T-cells convert to either exhausted or cytotoxic terminal states [40].

Bulk deconvolution
While scRNA-seq offers new opportunities to study skin cancer, bulk RNA-seq remains more common due to its lower cost and ease of use. Transferring findings from single-cell experiments to bulk datasets to take full advantage of available data resources is thus an important part of single-cell analysis. The existence of several publicly available skin cancer datasets with bulk microarray or RNA-seq data allows researchers to validate findings from single-cell experiments in larger cohorts of bulk data such as The Cancer Genome Atlas (TCGA) [58]. One common approach is to derive cell-type or tumor subpopulation-specific gene signatures from single-cell data and then score bulk data using these signatures, allowing signature activity to be correlated with clinical outcomes such as overall survival or metastasis [34]. More advanced deconvolution algorithms such as CIBERSORTx, MuSiC, and others use single-cell or bulk reference profiles to estimate cell-type abundance and gene expression from bulk tumor samples [59][60][61]. In the skin cancer setting, deconvolution has successfully been applied to study tumor-stromal interactions. Tirosh and colleagues used gene signatures derived from single-cell melanoma transcriptomes to deconvolute 471 TCGA melanoma samples and showed that CAF abundance correlates with AXL expression [34]. They also identified a set of CAF-expressed genes that correlates with T-cell infiltration [34]. Similarly, after deriving signatures for T-cell exhaustion and cytotoxicity from single-cell data and scoring TCGA melanoma samples with these signatures, Deng et al. identified cytotoxic and exhausted subpopulations that predict prognosis [40]. Other studies have used deconvolution algorithms to link levels of infiltrating immune cells with survival and nominate genes that may be involved in this process [62][63][64]. Deconvolution was also been applied to cSCC by Ji et al., who linked TSK-like expression with CAF activity in several TCGA cancer types [37].

Spatial context
Single-cell approaches are increasingly being combined with spatial transcriptomics (ST) to gain further insight into the complex nuances of different cell types and states. By providing tissue context, ST enables researchers to understand gene expression as it relates to cell position within a tissue. The two widely used methods for ST are imaging-based spatially resolved transcriptomics and in situ barcoding-based spatially resolved transcriptomics. Imaging-based spatially resolved transcriptomics methods such as MERFISH, seqFISH and ABER-FISH relies on multiplexed fluorescence in situ hybridization (FISH) and expansion microscopy to detect several RNA molecules [65][66][67][68]. In situ barcoding-based spatially resolved transcriptomics methods such as 10X Visium, Slide-seq and high-definition spatial transcriptomics (HDST) utilizes DNA barcoding to map out gene expression [68][69][70][71]. Although these methods can interrogate the entire transcriptome, they are limited by their tissue resolution, which varies from 55 μm to 2 μm. ST methods have been applied to skin cancer to localize tumor subpopulations and may shed light on their role in cancer progression (Fig. 4). Using a highly multiplexed and sensitive platform called CEL-Seq2 Baron and colleagues used scRNA-seq to characterize zebrafish melanomas and identified a tumor subpopulation that expresses a stresslike transcriptional program [72]. ST was used to demonstrate that these cells are enriched in tumor regions, but not in the normal surrounding tissue [72]. After confirming the presence of stress-like cells in human melanomas and other cancer types, the authors showed that the stressed state efficiently seeds new tumors in zebrafish and is associated with resistance to MEK and BRAF inhibitors. In another study, Thrane et al. used the 10X Visium ST platform to characterize melanoma lymph node metastases, highlighting differential gene expression profiles as well as multiple melanoma signatures across the same region in a single tumor biopsy [73]. Poor patient survival was associated with greater expression heterogeneity in the area of transition between melanoma and lymphoid tissue [73]. Ji et al. similarly used the 10X Visium ST platform to demonstrate that cSCC TSKs and basal cells localize to the tumor leading edge and confirm the presence of vascular and CAF-enriched transcripts in the surrounding fibrovascular niche [37]. By incorporating a novel multi-parameter tissue imaging workflow named Pick-Seq, Nirmal et al. showed that invasive cutaneous melanoma possesses a unique cellular microenvironment at the tumor-stromal boundary that harbors different cell types supported by cytokine gradients [74].
Some of the limitations of ST are that it lacks single-cell resolution and the read quality as well as the number of transcriptomes captured is often tissue-dependent. HDST notably enables users to reach a resolution of 2µM [68,70], providing the possibility of subcellular resolution. This method currently requires specialized analytics as well as bioinformatic expertise, and is limited by low sensitivity of mRNA capture [75]. Moreover, both scRNA-seq and ST do not consistently detect genes with lower expression. To overcome these barriers, approaches such as RNA in situ hybridization, which offers a targeted but more sensitive alternative, can be used to generate spatial information at single-cell resolution. A recent study by Tran et al. developed a multimodal strategy, Spatial Transcriptomic and RNA in situ Hybridization (STRISH), to study cancer-immune cell crosstalk at the genome-wide level in keratinocyte cancers [76]. This approach starts by inferring ligand-receptor interactions in the cancer-immune cell context using scRNA-seq and ST. Cell-cell interactions are then visualized by RNA in situ hybridization and quantitated using digital droplet PCR [76]. Using this analysis pipeline, the authors demonstrated co-expression of the IL34-CSF1R and THY1-ITGAM ligand-receptor pairs in cancer cell nests and areas of immune cell infiltrate [76]. We predict that similar strategies combining different technologies to overcome the limitations of each approach in isolation will be increasingly adopted to discover and confirm critical cell-cell interactions in the cancer setting.

Conclusions
Bulk RNA-seq remains a valuable method for defining the coding and non-coding transcriptome; however, when applied to the skin and its attendant malignancies, scRNA-seq has revealed robust expression differences between cell types and states, providing new insight into the cellular diversity and heterogeneity that occurs in these settings (Table 1). While this approach is accompanied by its own limitations, most notably dropout events and the scalability of existing data analysis methodologies, the aforementioned studies demonstrate how scRNA-seq has enabled an improved understanding of molecular events regulating skin cancer progression. We predict that future studies will utilize singlecell transcriptomics to interrogate larger cohorts of treatment-resistant and relapsing skin cancers, with an eye towards predicting therapy response. As sequencing costs fall, scRNA-seq may also increasingly be used in the clinical setting to devise targeted personalized treatment strategies for advanced skin cancers that cannot be managed surgically.

Declaration of Competing Interest
None.