Single-Cell RNA Sequencing Reveals Cellular Heterogeneity in an Acral Amelanotic Melanoma After Immunotherapy Treatment

Introduction

Immune checkpoint molecules, such as programmed death ligand 1 (PD-L1) and cytotoxic T lymphocyte antigen 4 (CTLA-4), negatively regulate T cell function by binding to receptors on activated T lymphocytes and other immune cells.¹ Immune checkpoint inhibitors (ICIs) are commonly used in clinical practice; they enhances anti-tumor immunity by blocking the inhibition of T cell activation pathways.² Furthermore, PD-L1 inhibitors, such as atezolizumab, have been approved as alternative ICIs for advanced cancers, such as urothelial carcinoma and melanoma.^3,4

Patients with metastatic melanoma treated with ICIs still risk developing new primary cutaneous melanoma (CMs).⁵ For instance, in a cohort of more than 2000 patients with metastatic melanoma treated with ICIs, the incidence of new primary CM was 1.9%, the 5-year cumulative risk was 4.9%, and the incidence rate was 1103 cases per 100,000 person-years.⁶ However, new primary CMs occurring during ICI treatment for other malignancies have not been adequately discussed.

Amelanotic melanoma (AM) is a relatively rare subtype accounting for 2–8% of all melanomas,⁷ and it is more rare in patients undergoing ICI therapy. AMs are prone to misdiagnosis owing to the absence of black color, resulting in progression to a more advanced tumor stage with a worse prognosis.⁸ Furthermore, increasing age and ultraviolet exposure predispose individuals to AM.⁹ Therefore, most AMs have been detected on chronically photodamaged skin of older patients.

Herein, we present a case of an acral AM that occurred during anti-PD-L1 therapy for urothelial carcinoma. Then, given the rarity of AM occurring during ICI treatment, we investigated the cellular heterogeneity of the acral AM using single-cell RNA sequencing (scRNA-seq).

Methods

Ethical Approval

This study was conducted following the principles of the Declaration of Helsinki with approval from Peking University First Hospital (Ethics Approval number: 2021[076]). The patient provided informed consent for publication.

Clinical Information

An 82-year-old man presented at our clinic with a 2.0×1.5 cm red lesion on his left palm in January 2021. In December 2018, the patient was diagnosed with advanced urothelial carcinoma and started continuous atezolizumab therapy. In December 2019, the patient noticed a red macule on his left palm. However, it was misdiagnosed as dermatitis and a skin infection; topical steroids and antibiotics were administered without improvement. In November 2020, the patient withdrew from atezolizumab therapy owing to fungal pneumonia. After that, the lesion regionally developed a black color. At this point, the patient came to our clinic, where we performed dermoscopic and biopsy examinations (Figure 1). Finally, the patient was diagnosed with acral AM. The tumor was surgically removed by excising 2 cm of tissue based on the Breslow thickness. A sentinel lymph node biopsy indicated no metastases in the left axillary. The patient was clinically staged as IIB (T3bN0M0). scRNA-seq of the spare resected lesion was also performed. The patient died six months after the surgery from multiple organ metastasis of urothelial carcinoma.

Figure 1 Gross appearance of a red lesion on the left palm of an 82-year-old man who previously underwent immunotherapy; the lesion has a heterogeneous color (black, white, and red).

Dermatopathology Analyses

Standard laboratory protocols were followed when processing formalin-fixed and paraffin-embedded surgical samples. Additionally, routine hematoxylin and eosin stains for the sections were performed. Protein expression was detected with primary antibodies targeting MelanA (IR63361-2, dilution: 1:500, DaKo, Santa Clara, CA, USA), HMB45 (N154587-2, dilution: 1:500, DaKo), and Ki-67 (IR62661-2, dilution: 1:500, DaKo). The secondary antibodies were goat anti-rabbit IgG (ZB-2308, dilution: 1:1000, ZSBIO, Beijing, China) and goat anti-mouse IgG (ZB-2313, dilution: 1:1000, ZSBIO). The Ki-67 positivity rate was assessed by two pathologists independently. The number of positive tumor cells was manually counted and then compared to the total number of tumor cells to calculate the Ki-67 positivity rate.

Melanoma Tissue Processing for Single-Cell Suspension

A fresh sample was collected at the Department of Dermatology of Peking University First Hospital. Tumor samples were taken from the operating room, placed in MACS Tissue Storage Solution (cat# 130-100-008, Miltenyi Biotec, Bergisch Gladbach, Germany) on ice, and transported to the laboratory within 45 minutes. Next, whole skin layer sections of melanoma tissues (approximately 2×1 cm) were washed three times with phosphate-buffered saline (PBS). Subsequently, the tissue was digested in 0.125 mg/mL of trypsin at 37 °C for 15 minutes. Then, 2% fetal bovine serum (FBS) in PBS was added to flush the tissue, and all fluid was collected in the centrifuge tube, which was then centrifuged at 1000 rpm for 3 minutes. After discarding the supernatant, the cells were suspended in 2% FBS-PBS. Next, a digestion mixture (2 mg/mL collagenase IV, 2.5 mM CaCl₂, 0.1 mg/mL elastase, 2% FBS, and DNA enzymes) was added to undigested bulk tissue cut into pieces as small as possible. Then, the tissue mixture was placed on a shaker for 30 minutes for digestion. After, the cells were centrifuged and resuspended in 2% FBS-PBS. Finally, the suspension was passed through 100 µm and 40 µm cell strainers in sequence, collecting the entire volume.

The quality was assessed using the Countstar cell analysis systems (Countstar Rigel S2, Alit Biotech, Shanghai) after centrifuging and resuspending the cells. Cell viability was obtained using the trypan blue method with the Countstar machine and intelligent image recognition technology in one step. If the cell viability was above 80%, 10x Genomics sample processing was performed.

Cell Capture and Complementary DNA (cDNA) Synthesis

Using the Single-cell 3′ Library and Gel Bead Kit V3.1 (cat# 1000075 10x Genomics, Pleasanton, CA, USA) and Chromium Single Cell B Chip Kit (cat# 1000074, 10x Genomics), the cell suspension (300–600 living cells per µL determined by Count Star) was loaded onto the Chromium Single Cell controller (10x Genomics) to generate single-cell gel beads in the emulsion according to the manufacturer’s protocol. In short, single cells were suspended in PBS containing 0.04% bovine serum albumin. Approximately 10,000 cells were added to each channel; the target number to recover was around 18,447 cells. Captured cells were lysed, and the released RNA was barcoded through reverse transcription in individual gel beads-in-emulsion (ie, GEMs). Reverse transcription was performed on an S1000TM Touch Thermal Cycler (Bio-Rad, Hercules, CA, USA) at 53 °C for 45 min, followed by 85 °C for 5 min, and then held at 4 °C. The cDNA was generated and then amplified, and quality was assessed using an Agilent 4200 (performed by CapitalBio Technology, Beijing, China).

scRNA-Seq Library Preparation

Following the manufacturer’s instructions, scRNA-seq libraries were constructed using the Single Cell 3′ Library and Gel Bead Kit V3.1 (10x Genomics). The libraries were sequenced using an Illumina NovaSeq 6000 sequencer with a sequencing depth of at least 100,000 reads per cell with a pair-end 150 bp (PE150) reading strategy (performed by CapitalBio Technology).

Cell Ranger

Alignment, filtering, barcode counting, and unique molecular identifier counting were performed with Cell Ranger software (v5.0.1) to generate feature-barcode matrices. The reads were aligned to the human reference genome GRCh38. The Cell Ranger software was obtained from the 10x Genomics website (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger).

R Toolkit: Seurat

Seurat is an R toolkit for single-cell genomics¹⁰ where the SCTransform command replaces the NormalizeData, ScaleData, and FindVariableFeatures commands; this change helps identify rare cell populations by improving the signal-to-noise ratio.¹¹ Cells with ≤400 and ≥6000 genes and a mitochondrial gene percentage of ≥20% were filtered out to include only high-quality true cells. After the quality control assessment, 22,077 high-quality cells were obtained. A principal component analysis was performed on approximately 3000 variable genes. We used the FindClusters function on 50 principal components with a resolution of 0.1 to perform the cluster analysis, annotating each cluster with known markers. For visualization, the dimensionality of each dataset was further reduced using Uniform Manifold Approximation and Projection (UMAP) with Seurat functions RunUMAP. Analysis pipeline has been published in GitHub. The data can be access by the link: https://github.com/TianJie327/acral-amelanotic-melanoma.

Gene Ontology (GO) Biological Process Enrichment Analysis

A GO biological process enrichment analysis implemented in the clusterProfiler (v4.0.5) was used to elucidate the biological properties of the cell clusters at single-cell levels.¹² Differentially expressed genes among three melanoma cell subtypes were detected with the FindMarkers function of Seurat. The significance threshold was set to an adjusted P-value of <0.05 and the top 25% log2-fold change.

Pseudo-Time Cell Trajectories Analysis

Single-cell pseudo-time trajectories were constructed with Monocle 2 (v2.20.0).¹³ The order was based on genes with a q-value of <0.01 that showed expression differences between clusters of melanomas selected via the unsupervised dpFeature procedure. Then, the data space was reduced to two dimensions with the DDRTree method. The cells were ultimately ordered in pseudo-time, and pigmentation melanomas represented the end of the trajectory.

Branched Expression Analysis Modeling (BEAM)

We used BEAM, a function of Monocle 2 (v2.20.0), to present the different fates of melanoma and pigmentation-related (MITF, MLANA, PMEL, and TYR), endosome-related (RAB4A, RAB5A, RAB7A, and RAB11A), and lysosome-related (LAMP1, LAMP2, and M6PR) genes on the branch of fate.

Results

Dermoscopic Examination

The dermoscopic examination showed black globules without a pigment network against a red background in the black region (Figure 2A), a scar-like depigmentation area in the white region (Figure 2B), and a homogeneous red background in the red region (Figure 2C).

Figure 2 Dermoscopic examination results in our 82-year-old male patient. (A) The black region shows black globules without a pigment network against a red background. (B) The white region shows a scar-like depigmentation area. (C) The red region shows a homogeneous red background.

Pathological Examination

Microscopic findings revealed two types of epithelioid tumor cells: pigmented melanoma cells and AM cells (Figure 3A). The pigmented cells were large, contained melanin granules, and were arranged as large tumor nests. The regional nests corresponding to black globules were identifiable on dermoscopic examination. The AM cells without melanin granules were small and scattered throughout the dermis, corresponding to the red background under the dermoscope.

Figure 3 (A) Two types of melanoma cells: pigmented cells containing melanin granules (yellow triangles) and small amelanotic cells distributed throughout the entire dermis (black triangle) (hematoxylin and eosin; 2× magnification). (B) Melan-A-positive melanoma cells (2× magnification). (C) HMB45-positive pigmented melanoma cells (2× magnification). (D) Pigmented melanoma cells (yellow triangle) have a higher KI-67 positivity rate (>10%) than amelanotic melanoma cells (<1%; black triangle) (2× magnification).

All melanoma cells were Melan-A-positive (Figure 3B); however, only pigmented cells were HMB45-positive (Figure 3C). The KI-67 positivity rate was significantly higher in the pigmented cells (>10%) than in the AM cells (<1%). The pigmented cells also demonstrated higher proliferative and invasive abilities than the AM cells (Figure 3D).

scRNA-Seq Profile Findings

UMAP

We identified seven major cell types: Melanoma (MLANA, MITF, TYR, and PMEL), epidermal appendage keratinocytes (KRT8 and KRT18), epidermal keratinocytes (KRT5 and KRT14), endothelial (PECAM1 and VWF), cancer-associated fibroblasts (FAP and MME), fibroblasts (ACTA2, DCN, and LUM), and immune cells (PTPRC) (Figure 4).

Figure 4 A uniform manifold approximation and projection (UMAP) plot shows the unbiased melanoma tissue classified into seven cell types.

GO Biological Process Enrichment Analysis

The scRNA-seq analysis revealed cellular heterogeneity; thus, the melanoma cells were divided into three clusters. GO term enrichment indicated that cluster 1 was active in RNA splicing processes, cluster 2 was active in drug responses and extracellular matrix modulation, and cluster 3 was active in melanin synthesis (Figure 5A). The expression of melanin synthesis-related genes was downregulated in cluster 1 and cluster 2 but upregulated in cluster 3 (Figure 5B). We defined clusters 1, 2, and 3 as amelanotic cell cluster 1, amelanotic cell cluster 2, and pigmented cell cluster, respectively.

Figure 5 (A) Gene ontology (GO) analysis defines three clusters based on the biological processes: 1) melanoma cell cluster 1 is active in RNA splicing; 2) melanoma cell cluster 2 is active in drug responses and extracellular matrix modulation; and 3) melanoma cell cluster 3 is active in the melanin synthesis process. (B) Single-cell transcript expression levels of MITF, MLANA, TYR, and PMEL in the three clusters.

Pseudo-Time Trajectories Analysis

scRNA-seq and the pseudo-time cell trajectory analysis revealed the cluster’s transition. Amelanotic cluster 1 was in the initial part of the trajectory, and the pigmented cluster was at the end, transitioning from amelanotic cluster 2 (Figure 6).

Figure 6 Pseudo-time trajectory showing the distribution of cells from melanomas based on the cell type (left). The pseudo-time trajectory is divided into three different states by Monocle 2 (middle). The pseudo-time color gradient transitions from dark (right; start) to light blue (end).

BEAM

To elucidate the molecular basis of the transformation, we explored the genes that determine the cell fate of the pigmented melanoma and AM clusters (Figure 7A). The melanin synthesis-related genes (MITF, MLANA, PMEL, and TYR) were downregulated, but lysosome and endosome-related genes (LAMP1, LAMP2, RAB4A, and M6PR) were upregulated in the amelanotic clusters compared to the pigmented cluster. Furthermore, cell cycle genes (RB1, PCNA, CDK2, CDK4, and CDK6) were upregulated in the pigmented cluster compared to the amelanotic clusters, consistent with the histological finding that the pigmented melanoma cells acquired higher division activity and were more aggressive (Figure 7B).

Figure 7 (A) Branched expression analysis modeling of the pigmented and amelanotic melanoma clusters. (B) Single-cell transcript expression levels of cell cycle genes. RB1, PCNA, CDK2, CDK4, CDK8, and CENPF were significantly upregulated in the pigmented cell cluster.

Discussion

Several reports have described BRAF inhibitor-associated cutaneous squamous cell carcinoma and T-cell lymphoma during checkpoint inhibitor therapy.^14–16 Furthermore, literature reviews were performed to determine if AM has been reported during immunotherapy.^17,18 However, none have investigated cellular heterogeneity in acral AM during or after immunotherapy treatments.

Cancer cells exhibit different transcriptome characteristics, which is one type of intra-tumoral heterogeneity in response to systemic therapy.¹⁹ Our scRNA-seq profile provided a model to interpret how cancer heterogeneity contributes to immunotherapy resistance.²⁰ Our cell trajectory results identified amelanotic cell cluster 1 at the beginning of the trajectory. Therefore, we propose that this was a group of melanoma cells active in the RNA splicing process. Thus, in response to the PD-L1 drug, an immune escape effect was exerted by modulating cellular transcriptomic changes and decreasing MITF expression through active RNA splicing, thereby reducing melanoma antigen expression (eg, TYR, MLANA, and PMEL).^21–23 A by-product of this response is the development of a non-pigmented clinical presentation. As a result, some researchers are exploring how RNA splicing processes affect cancer immunotherapy resistance and are using spliceosome-targeted therapies combined with immunotherapies to try and reduce tumor immunotherapy resistance.^24,25 Furthermore, after RNA splicing, cluster 1 transformed into cluster 2, and a large population of amelanotic cells exhibiting drug resistance emerged (ie, cluster 2 with GO biological process enriched in drug responses and extracellular matrix modulation).²⁶ Finally, the pigmented cell cluster appeared at the end of the trajectory. After discontinuing the medication, the melanoma cells transitioned from a therapy-resistant state to a proliferative and differentiated state (ie, the pigmented cell cluster).²⁷ The clinical course of our patient’s lesion supports this theory; after discontinuing the drug, an unrestricted tumor cell state and color appeared. Thus, the pseudo-time cell trajectory analysis was consistent with the clinical changes.

In our patient, the AM cell cluster underwent phenotypic transition to pigmented melanoma cell cluster after immunotherapy discontinuation. Therefore, we suspect the tumor had phenotypic plasticity. Phenotypic plasticity refers to a wide range of phenomena, including dedifferentiation and transdifferentiation (ie, a direct fate change to another differentiated cell type).²⁸ Phenotypic plasticity in cancer cells is one form of intra-tumoral heterogeneity in response to systemic therapy, which has been broadly reported.¹⁹ Although we did not confirm the exact association between phenotypic transition and immunotherapy, we found downregulated expression of melanin synthesis-related genes and upregulated expression of lysosome-related genes in AM cells. Melanosome is a lysosome-related organelle whose synthesis process is closely related to lysosome function, and PD-L1 antibody (H1A, STM108) induces PD-L1 degradation in lysosomes.^29,30 Thus, we hypothesize that increased PD-L1 degradation after medication use contributes to the upregulated expression of lysosomal function genes in AM cells. However, these hypotheses should be confirmed by cell and functional experiments in the future.

In our case, the acquired proliferative ability accompanied by discontinuation of the medication uncovered unanticipated biological principles dictating the aggressive disease state of the pigmented melanoma cells. The pigmented melanoma cells in this acral AM presented higher proliferative capacity, such as high Ki-67-positive cells and upregulated cell cycle gene expression. This finding contradicted the previous viewpoint of higher proliferation and invasiveness of AM compared to pigmented melanomas.^8,31 However, we do not think this result is inconsistent because more sophisticated tools, such as single-cell sequencing, were used to identify intra-tumor heterogeneity. Different environments in which the tumor exists (eg, continued or discontinued immunotherapy) might cause intra-tumor proliferative heterogeneity. However, this is not equivalent to inter-tumor heterogeneity.³²

The present study has some limitations. Melanomas are highly heterogeneous; thus, more samples should be accumulated to dissect the heterogeneity of ecology systems of melanoma occurring after immunotherapy exposure. Furthermore, we could not obtain lesion tissue before our patient discontinued atezolizumab; therefore, further studies are needed to confirm the findings and uncover the cancer immune editing mechanisms.

Conclusions

For patients exposed to immunotherapy, coexisting amelanotic and pigmented melanoma cells indicate prominent cellular heterogeneity in acral AM. Furthermore, the pigmented melanoma cells acquired higher proliferative ability than the AM cells in the evolution sequence. The atypical appearance of this tumor should remind clinicians, both oncologists and dermatologists, to monitor persistent skin lesions occurring during tumor immunotherapy, including non-skin cancers.