Cell line treatment model, genomic landscape of melanoma patients and tumor heterogeneity revealed by scRNA-seqWe procured 18 tumor samples from 15 individuals for this study (Supplementary Table 1). Four of these patients had undergone neo-adjuvant treatment with immunotherapy alone (either PD1 or combined PD1 + CTLA4 blockade), while one of the patients underwent neoadjuvant combined therapy with anti BRAF and MEK inhibition, as well as anti PD1 therapy. For two of these patients (1451 and 1511) we used a matched primary (cutaneous) and metastatic tumor (regional lymph node metastasis), while for another patient (1199) we analyzed 2 tumors from different metastatic regional lymph nodes. The remaining samples were either regional lymph node metastases (8 patients) or distant sites (2 adrenal, 1 small bowel, and 1 distant lymph node) of metastatic disease (Fig. 1b and Supplementary Table S1). All tumors were subjected to single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing. Whole exome sequencing (WES) was performed on paired tumor and normal genomic DNA. In addition, cell lines were created and propagated in vitro for three samples (Fig. 1a).Fig. 1Study design, datasets, genomic landscape and cell populations revealed by scRNA-seq. (a) Establishment of melanoma patient-derived cell lines and BRAF-MEK inhibitor resistant cell lines. (b) Landscape of driver mutations, copy number variations (CNV), mutation signatures, melanoma cell percentage, clinical information and data availability across 15 patients, with vertical black lines separating patients and dashed lines separating 2 metastatic samples taken from the same patient. Copy number amplification/gain and copy number deletion/loss are shown in red/light red, blue/light blue respectively. NA = not applicable, P = primary, M = metastasis, CL = cell line, CL_BM = BRAF/MEK inhibitors treated cell line, scRNA-seq = single-cell RNA sequencing, WES = Whole Exome Sequencing. (c) UMAP plot shows integrated tumor and its microenvironment cells from 8 samples. (d) UMAP plot shows integrated tumor and its microenvironment cells from 8 samples as in panel c, colored by main cell types. CAF = cancer associated fibroblasts, DC = dendritic cells, pDC = plasmacytoid dendritic cells, and NK = natural killer cells. (e) UMAP plot shows integrated immune cells and CAFs from 8 samples (as in panel c and d, colored by cell subtypes. Melanoma cells, erythrocytes and gastrointestinal cells have been removed. (f) Bar chart shows the fraction of immune cell subsets for each sample with color indicating cell types.Melanoma exhibits higher mutation rates than most other cancer types11. We detected a median of 164 coding mutations from tumor WES data (ranging from 12 to 523), with the variant allele fraction (VAF) largely consistent between WES and bulk RNA-seq (Fig. 1b and Supplementary Fig. 1). The frequency of hotspot mutations is also concordant with other studies9,12. In our dataset, 65% (11/17) of samples have BRAF V600E mutations and 24% (4/17) have NRAS Q61 mutations9 (Fig. 1b and Supplementary Fig. 1). Somatic copy-number variation profiles identified gains of oncogenes, including MET, MITF, and MYC and loss of tumor suppressor genes, including TP53, PTEN, CDKN2A, and CDKN2B. Mutations and copy number variations were very similar between patient tumors and their corresponding cell lines, suggesting the latter are good representatives of patient tumors. Interestingly, we observed VAF changes during disease progression for several mutations in cancer genes. For example, BRAF-V600E in patient 1451 expanded from 5 to 22%, while CTNNB1-S37Y receded from 15 to 2% from primary to metastasis (Supplementary Fig. 1).Next, we addressed the effect of misrepair of UV-induced DNA damage as a booster of melanoma driver mutations, namely C > T (by UVB) or G > T (by UVA). Out of the 5020 coding mutations, 80% were C > T (75%) or G > T (5%) mutations, highly characteristic of UVB/UVA inducement13. For each sample, we detected a median of 72% C > T (ranging from 17 to 89%) and 6% G > T (ranging from 2 to 20%). Consistent with prior research14, the abundance of C > T mutations was essentially recapitulated in cell lines, being only slightly lower than that in patient tumors (68% to 66% in 762_M, 89% to 84% in 1199_M1, and 88% to 86% in 1199_M2). These estimates were not revealed by previous studies comparing the Cancer Genome Atlas (TCGA) patient melanoma tumors to melanoma cell lines14,15.To examine tumor heterogeneity at single-cell resolution, we isolated 52,384 melanoma cells from scRNA-seq data of 17 patient tumors (almost no tumor cells detected in 1511_M and 1265_M, Supplementary Table S2). UMAP analysis reveals that melanoma cells cluster based on tumor of origin (Supplementary Fig. 2a). Melanoma cells from patient 1199, taken from two different metastatic regional lymph node deposits, not unexpectedly clustered together. We next sought to identify transcriptional differences among tumor subpopulations by calculating Pearson correlation coefficients of average expression of top variable genes across tumor subclusters followed by hierarchical clustering. Compared to subpopulations from different samples, most within-sample tumor clusters were highly correlated (Supplementary Fig. 2b). Notably, patients with the same driver mutations (especially NRAS) tend to be clustered together, highlighting the impact of driver mutations on transcriptome profiles of melanoma cells.To comprehensively investigate cell population composition in melanoma, we integrated 52,739 cells from 9 samples processed by the same kits (Chromium Single Cell 3′ Reagent Kits v2 chemistry). In addition to 24,235 melanoma cells, we identified 17,320 T cells, 6002 B cells, 3761 macrophages or monocytes, 105 dendritic cells (DC), 501 plasmacytoid dendritic cells (pDC), 321 cancer associated fibroblasts (CAFs), and other cell types. In line with previous single-cell melanoma study (4645 cells)10,16,17, we recapitulated that immune cells overlapped well between samples, while tumor cells showed sample specific clusters (Fig. 1c and Fig. 1d). We also observed that low tumor content might slightly affect clustering of melanoma cells. For example, tumor cells from samples with fewer melanoma cells (1232_M and 1239_M) formed clusters adjacent to the clusters composed of samples with high abundance of melanoma cells (762_M and 1143_M). Yet, samples with low tumor content still formed their own independent clusters while melanoma cells of different tumors from the same patients (1199_M1 and 1199_M2) were clustered together, indicating tumor heterogeneity is well captured by our scRNA-seq data regardless of abundance of melanoma cells.Next, to further investigate tumor microenvironment, immune cells and stromal cells were isolated and reclustered with melanoma cells, erythrocytes, and intestinal cells removed (Fig. 1e). Gene signatures of each immune subtype are shown in Supplementary Fig. 2c. Then, we identified 26,490 melanoma cells and 48,503 immune cells from the remainder of samples (Chromium Single Cell 5’ Reagents Cell Kit) by following the same strategy (Supplementary Fig. 2d and Supplementary Fig. 2e). In summary, there were 28,386 CD4 + naive T cells, representing the vast majority of CD4 + T cells (85%) in 18 samples (Fig. 1f). In CD8 + T cells, we identified CD8 + exhausted T cells and CD8 + cytotoxic T cells, occupying 31% and 69% of CD8 + T cells, respectively. In addition, we detected 1,588 monocytes and 4,178 macrophages, with 83% of macrophages being M2 macrophages (Fig. 1f). Different patients have a range of cell type compositions, such as CAFs ranging from 7% in patient 1451 primary tumor to 83% in patient 1511 primary tumor. Also, cell subset frequency varies between tumors with different clinical features. For instance, the proportion of dendritic cells dropped from primary stage to metastatic stage in both patients 1451 (11.96% to 0.73%) and 1511 (8.15% to 0.35%), suggesting a potential critical role of DCs in melanoma tumorigenesis as revealed by previous work18.Single-cell CNV analysis and clonal evolutionTo explore heterogeneity and clonal structure of tumors, we characterized chromosomal copy number variations (CNV) of individual melanoma cells using inferCNV. Previous comparative genomic hybridization (CGH) analysis19 revealed that sites commonly gained among melanomas, including chromosomes 1q, 7p, 7q, 6p, 8q, 20q, and sites commonly lost, including 9p, 9q, 10p, 10q, 6q, 11q, were captured by our single-cell CNV analysis (Fig. 2a). More importantly, tumor subclones carrying different CNVs were observed in multiple samples. For instance, tumor cells from sample 1372_M were composed of multiple subclones: a large tumor subpopulation with gain of 17q and loss of 11q, a small subclone with loss of 10q, and another small subclone with loss of 13p. To further investigate clonal structures, we used Uphyloplot2 to visualize the evolution of subclonal CNV events identified by infercnv hidden Markov models (HMM) subcluster CNV predictions algorithms (https://github.com/harbourlab/UPhyloplot2). Interestingly, we found gains of 1q and 7q often occurred at early stages, being observed in 67% (10/15) and 60% (9/15) metastatic tumors, respectively, and highlighting the importance of these CNV events in melanoma.Fig. 2Single-cell copy number variation analysis of cutaneous melanomas. (a) Single-cell CNV landscape from 17 samples with hierarchical clustering of CNV profiles of melanoma cells within every tumor. Bar plot on the right indicates the number of melanoma cells in each tumor. (b) Clonality trees of driver CNVs for each sample. The length of branches corresponds to the number of cells in the calculated subclone. Copy number gain is labelled in red and copy number loss is labelled in blue.When investigating genetic aberrations associated with metastatic tumors from different sites of the same patient (1199_M1 and 1199_M2), we observed gain of 1q and 7q in early stages of both tumors. Interestingly, gain of 8q and loss of 6q, 9q, and 10q are in all tumor cells of 1199_M1, while these events are only present in subclones in 1199_M2 (branch point B), indicating the heterogeneity and development of CNV subclones between metastatic tumors within a same patient. Overall, this analysis reveals previously unappreciated complexity in CNV clonal development in cutaneous and regional lymph node metastatic melanoma.To investigate genetic changes associated with metastasis and tumor progression, we compared phylogenetic CNV trees between primary and metastatic samples (Fig. 2b). We observed gains of 1q and 15q in the resected metastatic tumor (1451_M) from patient 1451 that were not present in the primary tumor (1451_P), consistent with prior work19. Interestingly, we found evidence that initial gain of 15q was followed by later gain of 1q (1451_M).We next examined how clonal structure transcriptionally evolves from the primary to metastatic tumor using case 1451. To investigate whether metastatic clusters are descended from primary clusters, we first isolated melanoma cells and identified tumor subpopulations within the primary and metastatic tumors (Supplementary Fig. 3a-b). Next, we integrated data from the two time points, investigating how the respective cells cluster. In mapping tumor subpopulations from two timepoints to the integrated data, we found that some cells from the primary tumor formed an isolated cluster (Cluster 3), while others were mixed with tumor cells from the metastasis (Cluster 1, Supplementary Fig. 3c). Specifically, the vast majority of cells from subcluster 1 (1451_P_SC1) and a small percentage of cells from subcluster 0 (1451_P_SC0) formed a distinct cluster (Cluster 3), indicating they are transcriptionally different from tumor cells at the metastatic stage. Interestingly, most cells from subcluster 0 (1451_P_SC0) and subcluster 2 (1451_P_SC2), as well as some from cluster 1 (1451_P_SC1) in the primary stage were mixed with subcluster 2 in the metastatic tumor (1451_M_SC2), meaning these tumor subpopulations share similar transcriptome profiles.Next, we sought to delineate this transition by trajectory analysis (Supplementary Fig. 3d). Consistent with our observation in Supplementary Fig. 3c, initial cell populations (pseudotime point 1) are composed of two subsets, the one dominant clone originating from subcluster 1 in the primary tumor (1451_P_SC1) and the other dominantly derived from subcluster 0 and 2 (1451_P_SC0 and 1451_P_SC2). At pseudotime points 2 and 3, cells developed along different branches. Most cells from clone 1451_P_SC1 became the dominant clones observed, with their evolution ceasing. Proceeding further past pseudotime point 3, the remaining primary tumor cells are slowly overtaken by metastatic tumor cells (1451_M_SC2), which becomes the new dominant clone. These observations collectively suggest that 1451_M_SC2 likely descended from 1451_P_SC0 and 1451_P_SC2.Finally, to investigate genes associated with clonal changes from primary to metastasis, we compared expression profiles of different clusters in integrated data by conducting an unbiased differential expression gene analysis (Supplementary Fig. 3e). Interestingly, we found that differentially expressed genes (DEGs) in cluster 2, mainly composed of cells from 1451_M_SC1, are involved in cell cycle and mitosis (q value = 2.99e-10), while cluster 4, derived from 1451_M_SC3, has high expression of keratin related genes involving the keratinization pathway (q value = 0.0076). Of note, pathway analysis revealed cluster 1 is significantly involved in transcriptional regulation by AP-2 (q value = 0.0099). We next sought to examine the expression of AP-2 related genes, including TFAP2A, TFAP2B, and KIT in clusters 1 and 3, finding that cluster 1 has higher expression of these genes (Supplementary Fig. 3f., q values are < 2.23e-308, < 2.23e-308, and = 1.16e-20 respectively). More importantly, we found AP-2 related genes are downregulated in the metastatic tumor as compared to the primary tumor (Supplementary Fig. 3 g, q values are < 2.23e-308 for all 3 genes). Previous work has revealed that loss of AP-2 expression occurred with malignant transformation and tumor progression in cutaneous malignant melanoma20, which aligns with our observations in this study.Patient-derived cell lines and anti-BRAF/MEK treatment mechanismTo investigate mechanisms of BRAF/MEK inhibitor resistance in vitro, we created 3 patient derived cell lines from 3 tumors (2 patients, 762 and 1199) with BRAF mutant metastatic melanoma (Fig. 1a). As noted in Table S1, patient 1199 received no therapy during the year prior to diagnosis, while patient 762 had been on anti CTLA4 + PD1 therapy in the year prior to resection.Cell lines were first compared to their parental tumors. Interestingly, the immune responsive fraction seems to disappear in the cell lines as compared to the patient tumors, as shown by loss of tumor cells involved in TCR activation and the PD1 signaling pathway (Fig. 3a; 1199_M2 not included due to low number of tumor cells on single cell analysis). In other words, this reflects a mechanism whereby the loss of the immune microenvironment in cell lines induces the disappearance of the cell populations regulating PD1 signaling. Conversely, there is a relative enrichment of actively cycling cells in the cell line as compared to the matched patient tumor.Fig. 3Transcriptional changes of melanoma cells from patient tumors, to patient tumor-derived cell lines, to BRAF/MEK inhibitors resistant cell lines. (a) UMAPs show all cell populations in 762_M and 1199_M1 tumors (left), melanoma cells in 762_M and 1199_M1 with cells colored by related pathways (middle), and melanoma cells in tumor-derived cell lines (762_M_CL and 1199_M1_CL) with cells colored by related pathways (right). (b) UMAPs show integrated melanoma cells from tumor-derived cell lines (762_M_CL, 1199_M1_CL, and 1199_M2_CL) and BRAF/MEK inhibitor resistant cell lines (762_M_CL_BM, 1199_M1_CL_BM and 1199_M2_CL_BM), with cells colored by samples (left) or clusters (right). c. Violin plots show normalized expression of ALDOA and PGK1 in each cluster for sample 762_M (top), 1199_M1 (middle) and 1199_M2 (bottom). Clusters mainly composed of cells from BRAF/MEK inhibitors treated cell lines are highlighted in black boxes.Using these cell lines, we then demonstrated the ability of BRAF/MEK inhibition to cause melanoma tumor cell death in BRAF mutant tumors. We subsequently subjected the lines to 1 month of BRAF/MEK treatment in vitro, creating a treatment “resistant” line. This line was no longer susceptible to death with BRAF/MEK treatment (Supplementary Fig. 4a-d). We next looked at the characteristics of these treated cell lines to see if we could target any resistance mechanisms. By integrating cells from initial cell lines and BRAF/MEK treatment resistant lines, we found that some cells from resistant lines were mixed with cells from the initial tumor derived cell lines, while other cells were clustered separately. This suggested unique transcriptome features of the resistant cells (Fig. 3b).Taking case 762_M as an example, some BRAF/MEK treatment resistant cells (pink) were clustered together with cells from the initial cell lines (purple), suggesting transcriptional similarity. By contrast, most resistant cells have different transcriptome features from the initial lines, forming distinct clusters, including clusters 0, 3, 6, 7, and 9. Interestingly, we found two glycolysis related genes, aldolase A (ALDOA) and phosphoglycerate kinase 1 (PGK1), were highly expressed in the BRAF/MEK treatment resistant clusters, as compared to other clusters (Fig. 3c, q values are < 2.23e-308 and 1.28e-113, respectively, for 762_M). ALDOA might play an important role in developing melanoma chemoresistance by promoting the effect of angiopoietin-like 4 (ANGPTL4) on melanoma cell invasion21. PGK1 has also been revealed as a target sensitizing BRAF V600E mutant melanoma cells to BRAF inhibition22. Notably, the overexpression of ALDOA and PGK1 in clusters composed of resistant lines, is also observed in cell lines from 1199_M1 (q values are 2.14e-99 and < 2.23e-308, respectively) and 1199_M2 (q value are < 2.23e-308 for both ALDOA and PGK1). Our observation indicates that BRAF/MEK resistant cells are characterized by high expression of ALDOA and PGK1, and therefore the glycolysis pathway might be associated with the development of resistance.Targeting BRAF/MEK inhibitors resistant cell lines with metforminMetformin can suppress tumor cell proliferation and promote apoptosis of cancer cells by inhibiting the glycolysis energy mechanism23,24. We thus treated the BRAF/MEK resistant lines with metformin, which demonstrated partial resolution of treatment resistance (Fig. 4a). Briefly, for cell lines from patient 762, at both lower metformin concentrations, significantly increased cell death was seen in the BRAF/MEK resistant cell line as compared to the wild type cell line. For cell line 1199_M1 all concentrations of metformin led to significantly increased cell death in BRAF/MEK resistant cell lines as compared to the wild type line.Fig. 4Transcriptional changes of BRAF/MEK resistant cells before and after metformin treatment. (a) Bar graphs show metformin related cell death on initial lines (762_M_CL and 1199_M1_CL) and BRAF/MEK inhibitors treated cell lines (762_M_CL_BM and 1199_M1_CL_BM). 762_M is shown at the top and 1199_M at the bottom. CL_BM_MF = metformin treated cell lines. mM = millimolar, ** = p < 0.01, *** = p < 0.001. (b) UMAPs show integrated melanoma cells from tumor-derived cell lines (762_M_CL, 1199_M1_CL), BRAF/MEK inhibitors resistant cell lines (762_M_CL_BM, 1199_M1_CL_BM) and metformin treated cell lines (762_M_CL_BM_MF, 1199_M1_CL_BM_MF), with cells colored by samples (left) or clusters (right). (c) Violin plots show normalized expression of ALDOA and PGK1 in each cluster for sample 762_M (left), 1199_M1 (right). Clusters mainly composed of cells from BRAF/MEK inhibitors treated cell lines and/or metformin treated cell lines are highlighted in black boxes. (d) Violin plots show normalized expression of ALDOA and PGK1 of cells from 762 and 762_M_CL_BM_MF in boxed clusters of panel c. * = p < 0.05, ** = p < 0.01, *** = p < 0.001 (e) Violin plots show normalized expression of ALDOA and PGK1 of cells from 1199_M1_CL_BM and 1199_M1_CL_BM_MF in boxed clusters of panel c. * = p < 0.05, ** = p < 0.01, *** = p < 0.001 (f) Violin plots show normalized expression of representative DEGs for cell populations disappearing post metformin treatment, such as cluster 9 in 762_M and cluster 1 in 1199_M1, highlighted in black boxes. (g) Violin plots show normalized expression of representative DEGs for cell population emerging post metformin treatment, such as cluster 10 in 762_M, highlighted in black boxes.Next, we sought to investigate the characteristics of metformin treated lines by integrating cells from post metformin treatment, with cells from initial lines and BRAF/MEK resistant lines (Fig. 4b). Consistent with observations in Fig. 3c, we found ALDOA and PGK1 were highly expressed in clusters mainly composed of cells from resistant lines and metformin treated lines (Fig. 4c and Supplementary Fig. 5 a-b, q values are < 2.23e-308 and = 1.28e-113 respectively for 762_M; 2.72e-182 and 6.17e-279 respectively for 1199_M1). More importantly, we observed ALDOA and PGK1 were significantly downregulated in metformin treated lines as compared to the resistant lines in most tumor subpopulations for case 762_M (Fig. 4d) and in all subpopulations for case 1199_M1 (Fig. 4e). These observations suggest the exciting possibility that metformin could overcome BRAF/MEK inhibitor resistance, potentially by downregulating glycolysis related genes such as ALDOA and PGK1.In addition, we found that BRAF/MEK treatment resistant cells in 1199_M1 cluster 1 vanished after metformin treatment (Fig. 4b, Supplementary Fig. 5b). Likewise, cluster 9 in case 762_M represents a BRAF/MEK resistant population that vanished after metformin treatment (Fig. 4b, Supplementary Fig. 5a). To investigate the characteristics of these populations, we performed DEG analysis and found that Annexin A1 (ANXA1), Brain-expressed X-linked protein 1 (BEX1), and Leupaxin (LPXN) are specifically highly expressed in BRAF/MEK resistant populations, but disappeared after metformin treatment (Fig. 4f, q values are 9.28e-38, 1.03e-22, and 2.54e-70, respectively, in 762_M and < 2.23e-308, 5.56e-197 and 2.07e-70, respectively, in 1199_M1). Previous work has suggested ANXA1 promotes melanoma dissemination in primary tumors25, and has been shown to help sustain tumor proliferation and metabolism in gastric cancer26. Although there are no studies discussing the roles of BEX1 and LPXN on melanoma cell proliferation or drug resistance, downregulated LPXN expression was shown to suppress malignant proliferation in a human acute monocytic leukemia cell line27. Our findings identify these 3 genes as potential prognostic markers of BRAF/MEK resistant cells that could be targeted by metformin.Lastly, we investigated the characteristics of cell populations that emerged after metformin treatment, such as cluster 10 in case 762_M (Fig. 4b, Supplementary Fig. 5a). Interestingly, we found that MYC, a proto-oncogene, and several heat shock proteins, including DNAJB1, HSPA1A, and HSPA1B, were upregulated in cluster 10 in 762_M (Fig. 4g, q values are 7.45e-39, 4.79e-62, 6.17e-10, and 1.69e-12, respectively). Studies have found c-myc overexpression drives melanoma metastasis28, heat-shock protein (HSP) 70 mRNA encoded by HSPA1A and HSPA1B promotes tumor cell growth29, and HSP70 inhibition enhances the response to melanoma treatment with BRAF inhibitors30. DNAJB1, a member of the HSP40 protein family, could interact with HSP70 and induce its ATPase activity31. However, the overexpression of these genes was not observed in newly emerged clusters after metformin treatment in another case, cluster 3 of 1199_M1. This might be due to different genomic alterations and different treatment regimes between patients 762 and 1199. Together, our investigations characterized cell populations induced by metformin treatment and highlighted the important roles of MYC, DNAJB1, HSPA1A, and HSPA1B on melanoma, which should be further explored.
