
==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322
Nature Publishing Group UK London

39300183
72255
10.1038/s41598-024-72255-9
Article
Treatment resistance to melanoma therapeutics on a single cell level
Yao Lijun 12
Krasnick Bradley A. 3
Bi Ye 3
Sethuraman Sunantha 12
Goedegebuure Simon 3
Weerasinghe Amila 12
Wetzel Chris 3
Gao Qingsong 12
Oyedeji Abimbola 3
Mudd Jacqueline 3
Wyczalkowski Matthew A. 12
Wendl Michael 12
Ding Li lding@wustl.edu

1245
Fields Ryan C. rcfields@wustl.edu

34
1 https://ror.org/01yc7t268 grid.4367.6 0000 0004 1936 9350 Department of Medicine, Washington University in St. Louis, St. Louis, MO 63110 USA
2 https://ror.org/01yc7t268 grid.4367.6 0000 0004 1936 9350 McDonnell Genome Institute, Washington University in St. Louis, St. Louis, MO 63108 USA
3 grid.4367.6 0000 0001 2355 7002 Department of Surgery, Washington University School of Medicine, St. Louis, MO USA
4 https://ror.org/01yc7t268 grid.4367.6 0000 0004 1936 9350 Siteman Cancer Center, Washington University in St. Louis, St. Louis, MO 63110 USA
5 https://ror.org/01yc7t268 grid.4367.6 0000 0004 1936 9350 Department of Genetics, Washington University in St. Louis, St. Louis, MO 63110 USA
19 9 2024
19 9 2024
2024
14 2191528 1 2024
5 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Therapy targeting the BRAF-MEK cascade created a treatment revolution for patients with BRAF mutant advanced melanoma. Unfortunately, 80% patients treated will progress by 5 years follow-up. Thus, it is imperative we study mechanisms of melanoma progression and therapeutic resistance. We created a scRNA (single cell RNA) atlas of 128,230 cells from 18 tumors across the treatment spectrum, discovering melanoma cells clustered strongly by transcriptome profiles of patients of origins. Our cell-level investigation revealed gains of 1q and 7q as likely early clonal events in metastatic melanomas. By comparing patient tumors and their derivative cell lines, we observed that PD1 responsive tumor fraction disappears when cells are propagated in vitro. We further established three anti-BRAF-MEK treatment resistant cell lines using three BRAF mutant tumors. ALDOA and PGK1 were found to be highly expressed in treatment resistant cell populations and metformin was effective in targeting the resistant cells. Our study suggests that the investigation of patient tumors and their derivative lines is essential for understanding disease progression, treatment response and resistance.

Subject terms

Cancer
Computational biology and bioinformatics
National Cancer Institute, USAT32CA009621 U2CCA233303 Krasnick Bradley A. issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Melanoma affects over 90,000 patients per year in the United States and is implicated in nearly 10,000 deaths (Seer.gov). Targeted therapeutics involving the BRAF-MEK pathway for patients harboring BRAF mutant tumors has led to improved survival for advanced melanoma1,2. Yet, treatment resistance and progression remain expected long term outcomes, with four out of five patients progressing by five years post treatment initiation3. Mechanisms to overcome resistance have relied upon global melanoma signatures4,5.

The relatively recent addition of immunotherapeutics to the melanoma treatment armamentarium has increased treatment options, further prolonging survival for advanced melanoma patients6. This work is founded upon murine systems, but in the drive to understand mechanisms of human melanoma resistance to immunotherapy, human melanoma cells and cell lines are increasingly being utilized to validate findings in murine and other models7,8.

Bulk sequencing by TCGA facilitated classification of melanoma into distinct groupings9. More recently, single cell sequencing has enabled researchers to investigate cell to cell differences in gene expression profiles of both melanoma cells and infiltrating cells, with further therapeutic implications10. Further work is needed to leverage scRNA sequencing data to effectively target select tumor populations in human tumors.

Here, we combine bulk and single cell RNA sequencing to investigate profiles of individual tumors and infiltrating cells for patients with advanced melanoma. We then analyzed how the addition of single cell sequencing of patient derived cell lines may help us to rationally overcome treatment resistance.

Results

Cell line treatment model, genomic landscape of melanoma patients and tumor heterogeneity revealed by scRNA-seq

We 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. 1 Study 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 evolution

To 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. 2 Single-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 mechanism

To 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. 3 Transcriptional 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 metformin

Metformin 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. 4 Transcriptional 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.

Discussion

We utilize multiple melanoma tumor samples from 15 patients to perform a comprehensive bulk and single cell genomic analysis. As seen in prior work, melanoma tumors commonly harbor mutations, with our dataset tallying ~ 164 coding mutations per sample, with the range of 12 to 523. Tumor purity could affect variant calling accuracy and the number of mutations tends to be higher in samples with high purity that that in samples with low purity32. We noted a few samples (1199_M2, 1511, 1232, 1239, 1294) having only a few mutations detected, which could be mainly due to low tumor purity given low abundance of melanoma cells based on scRNA-seq (Table S2). We have tried to improve variant detection accuracy by using both WES and RNA-seq, however, variant calling results from low tumor purity samples should be carefully interpreted. In terms of mutation signatures, as expected, the majority of them corresponded to UVA/B induced damage. Notably, based on UMAP, despite the relatively low cell counts, 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 inter-patient tumor heterogeneity is well captured by our scRNA-seq data (Fig. 1c-d). Using single cell data, we showed that diverse patient tumors tended to cluster most strongly by patient of origin (Supplementary Fig. 2a). In addition, there was independent clustering based on tumor driver mutations (e.g. NRAS) (Supplementary Fig. 2b).

Clonality analysis based on cell-level copy number variants demonstrated that gains of 1q and 7q occurred early in the progression of metastatic melanoma tumors (Fig. 2b). 1q and 7q represented an early stage event in clonal evolution prior to significant divergence in over 67% and 60% of metastatic tumors, respectively. Studies have speculated on the timing of general genomic events in melanoma33, but the single-cell resolution of our CNV findings reveals the potential sequence of the development of CNV subclones that may contribute to tumor progression. We observed significant clonal variation within the same patient, demonstrating the complexity of the clonal evolution process. We also looked at the expression profiles of tumor subpopulations at the primary and metastatic stage. We observed expansion and reduction of tumor subclones, as well as emergence of new tumor subclones in the metastatic tumor (Supplementary Fig. 3). Interestingly, and consistent with prior literature, we used single cell trajectory analysis to show that AP-2 related genes are downregulated in metastatic melanoma as compared to the primary tumor20.

We next showed the ability to create patient derived cell lines, which genomically were similar to their parental tumor (Fig. 1b). Interestingly, the PD1 responsive tumor fraction that is crucial in vivo, especially in the era of immunotherapeutics, disappears when cells are propagated in vitro (Fig. 3a). This likely means that cell systems propagated outside of intact immune systems are no longer relevant to the tumor-immune system relationship. It appears that while the immune responsive fraction of cells seen in the cell lines went away, there was a corresponding increase in cells engaged in cell division and mitosis, a signal of upregulation of the respective pathways. Although intuitive, these findings we believe are novel.

For the ~ 50–60% of melanoma patients harboring a BRAF mutation, BRAF + MEK inhibition has been groundbreaking, with the vast majority of patients showing treatment response1,3. Unfortunately, 50% of patients progress by 6 months post-treatment initiation, 75% progress at 3 years, and 80% at 5 years3. Some of these patients have previously failed or are not candidates for immunotherapeutics targeting PD1 + /− CTLA434. Overall, a majority of patients will eventually develop progressive disease, with limited treatment options outside of clinical trials. Over the past decade different epigenetic and genetic mechanisms of resistance to anti BRAF therapy have been identified35.

Herein, we demonstrated the ability to utilize scRNA-seq to look at individual melanoma tumor populations resistant to therapy, and then target them rationally based on RNA expression profile. We saw upregulation of glycolysis related genes ALDOA and PGK1 in BRAF/MEK resistant tumor subpopulations (Fig. 3c). This supports prior work that suggested the potential roles of ALDOA and PGK1 on the chemoresistance mechanism21,22. Both ALDOA and PGK1 catalyze key steps in the glycolysis pathway. Metformin has previously been shown to inhibit ALDOA and PGK1, and has been found to be potentially beneficial in a variety of malignancies, both invitro and in vivo, including in melanoma36–41. In our samples, metformin treatment led to cell death in these treatment resistant cells beyond what was seen in wild type cells (Fig. 4a).

Excitingly, we found decreased expression of ALDOA and PGK1 in metformin treated samples, suggesting that metformin treatment can potentially rescue BRAF/MEK resistance by downregulating ALDOA and PGK1 (Fig. 4d and e). This provides rationale to utilize scRNA-seq to rationally develop therapeutic strategies in the setting of treatment resistance. This paves the way for real time treatment strategies using scRNA-seq in the setting of melanoma resistance to targeted therapeutics. We also highlighted other tumor subpopulations that disappeared with metformin treatment (populations overexpressing ANXA1, BEX1 and LPXN), as well as newly emerging populations after metformin treatment (overexpressing MYC, DNAJB1, HSPA1A, and HSPA1B). This provides both further prognostic markers for metformin treatment, as well as subpopulations that may emerge and prove to be resistant to metformin therapy. Metformin has been shown to be safe for and is already approved for patient use. We feel that this lays the groundwork to begin testing our hypothesis in a cohort of patients resistant to BRAF/MEK, and may prove to be a novel treatment option for this hard to treat group of patients.

Future studies could expand our findings by including more patients with longitudinal samples for tumor transcriptional trajectory analysis along the metastasis. Also, it would be beneficial to expand the number of patient-derived cell lines to validate and further investigate the BRAF/MEK resistance mechanisms. In regards to pretreatment heterogeneity of our samples, one of our cell lines highlighted for the paper underwent neoadjuvant therapy with combined anti-PD1 and CTLA4 blockade, while the other two had no neoadjuvant treatment. This difference in neoadjuvant therapy may affect future behavior of the tumor cells in-vitro, but nevertheless we found consistent trends when it comes to ALDOA and PGK1 expression, as described above. Given the rather high tumor mutation burden in cutaneous melanoma42, it is worthwhile to investigate cell-level mutation clonality and evolution by single-cell DNA-sequencing and scRNA-seq43,44. We have mapped somatic mutations to individual single cells using an in-house tool (https://github.com/ding-lab/10Xmapping), but no mutation-specific clusters were observed due to the low coverage of detected mutations. Lastly, further endeavors could integrate single-cell proteomics with current single-cell transcriptional level investigation, which will enable researchers to comprehensively study metastasis and drug resistance mechanisms from multi-modal perspectives.

In summary, we utilized bulk sequencing and scRNA-seq to highlight tumor heterogeneity, the mutational signatures and CNV clonality of melanoma. By taking advantage of scRNA-seq in combination with in-vitro cell treatment, we investigated mechanisms of BRAF/MEK treatment resistance and partially overcame this with metformin.

Methods

Patient cohort and dataset

Patients over 18 years of age with advanced/ metastatic melanoma were approached for consent for tumor collection and processing. All work was done under Washington University in St. Louis approved IRB research protocol, entitled Advanced Genetic and Molecular Analsysis of Solid Tumors. For all patients verbal and written informed consent was obtained. Patients undergoing resection of primary and/or metastatic melanoma at Barnes Jewish Hospital/ Washington University in St. Louis by surgeon and senior author RCF were approached for tumor collection and consent. This was done under our institutional review board approved protocol. Overall there were 15 total patients included, with 18 melanoma tumors. See Table S1.

Tumor collection and processing

Resected tumor samples were collected in RPMI 1640 media (Gibco, USA) on ice and processed immediately to dissociation for the purpose of preparing single cell suspension. Briefly, for single cell sequencing, tumor tissue was mechanically processed into small pieces followed by enzymatic digestion. Enzymatic dissociation media (EDM) was prepared as follows: 1 g collagenase (Sigma, MO, USA), 0.1 g DNAse I (Type IV; Sigma), 10 mL HEPES (10 mM), 2500U/1L Hyaluronidase (Sigma) in 1L RPMI. The suspension was filtered over a 100 micron filter (Sigma). The resulting cells were counted and resuspended as our single cell suspension.

For whole exome and RNA sequencing, tumor samples were collected and kept in RPMI collection media supplemented with Penicillin and Streptomycin immediately following resection/biopsy, and placed on ice. Samples were processed into multiple 1 × 1 mm chunks, followed by snap freezing of tissue chunks in liquid nitrogen and storage at – 80 °C. Total time from resection to freezing tumor pieces was limited to < 60 min. Alternatively, tumors were directly processed into a single cell suspension (see above).

Blood collection

Patients consented for tumor collection, were also consented for collection of peripheral blood. 30 ml of peripheral blood was collected and processed into peripheral blood mononuclear cells (PBMCs) via ficoll density gradient centrifugation. PBMCs were stored in liquid nitrogen, and subjected at a later date to DNA isolation followed by whole exome sequencing.

Cell lines and treatment

Approximately 1 million cells from total single cell digest as described above were plated on sterile 25cm2 tissue culture flasks (Sigma). Non-adherent cells were removed 24 h after initial plating. Every 2–3 days non adherent cells were removed, and the cells were then passaged in 1–3 weeks using 0.25% Trypsin–EDTA (Gibco, MA, USA) once ~ 75% confluent, for 5–10 total passages. RPMI (Gibco) + 10% FBS (Sigma) was used for cell culture media.

Cell lines (762_M, 1199_M1, 1199_M2) were subjected to BRAF/MEK inhibitors treatment for 4 weeks. Vemurafenib (anti-BRAF) and Cobimetinib (anti-MEK) were obtained (MedChemExpress, NJ, USA), and reconstituted per manufacturer instructions. To create resistant cell lines, Vemurafenib was used at a concentration of 0.25 mM and Cobimetinib at 0.005 mM. The concentration was determined by the combined dose that resulted in approximately 50% death of the non-resistant lines (762_M_CL as an example, Supplementary Fig. 4a-c). After the first treatment, subsequent weekly doses of Vemurafenib and Cobimetinib were changed to 0.125 mM and 0.0025 mM in order to promote BRAF MEK resistance, but also allow for some cellular proliferation for subsequent passaging. After 4 weeks of constant treatment, changing media and drug weekly, “resistant lines” were created.

For BRAF + MEK cell treatment of resistant lines (to confirm resistance, as in Supplementary Fig. 4d), Vemurafenib was used at a concentration of 0.25 mM and Cobimetinib was used at a concentration of 0.005 mM. Treatments were done for 96 h. BRAF-MEK resistant lines derived from 762 and 1199_M1 lines were further subjected to Metformin treatment (Med Chem Express) at increasing concentration levels ranging from 0.1 to 1 mM (Fig. 4a). Metformin treatment courses were 96 h. For post metformin treatment cell lines subjected to scRNA-seq analysis, a metformin concentration of 0.1 mM was used to treat 762_M_CL_BM, and 0.5 mM was used to treat 199_M1_CL_BM (these were selected due to their similar efficacy at causing melanoma cell death, ~ 20–40%). After 96 h of treatment, viable cells were counted and sent for scRNA-seq analysis.

Cell viability assay

To assess in-vitro cell viability a CellTiter-Glo Luminescent Viability Assay was employed (Promega, WI, USA). 10,000 cells per well were plated in a 96 well plate (Fisher Scientific, MA, USA). Treatments were employed at the prespecified concentrations and treatment time course. Cell viability was assessed after the treatment was complete, utilizing the kit instructions. Briefly, an equal volume of reagent was added to the cell media in each well. Plates were incubated for 10 min, and cell luminescence recorded. Luminescence is generated in this assay by Adenosine Triphosphate (ATP), which is proportional to the number of live cells in the assay.

Single cell library prep and sequencing

Single cell suspensions or cell lines (as described above) were counted using a hemocytometer. Cells were washed and resuspended in PBS (Sigma) + 0.04% bovine serum albumin (BSA; Fisher Scientific) at 1000 cells/ul.

Utilizing the 10 × Genomics (CA, USA) Chromium Single Cell 3’ v2 or 5' Library Kit and Chromium instrument, approximately 17,500 cells were partitioned into nanoliter droplets to achieve single cell resolution for a maximum of 10,000 individual cells per sample. The resulting cDNA was tagged with a common 16nt cell barcode and 10nt Unique Molecular Identifier during the reverse transcription reaction. Full length cDNA from poly-A mRNA transcripts was enzymatically fragmented and size selected to optimize the cDNA amplicon size (approximately 400 bp) for library construction (10 × Genomics). The concentration of the 10 × single cell library was accurately determined through qPCR (Kapa Biosystems) to produce cluster counts appropriate for the HiSeq 4000 or NovaSeq 6000 platform (Illumina). 26 × 98 bp (3' v2 libraries) or 2 × 150 bp (5' libraries) sequence data were generated targeting between 25 K-50 K read pairs/cell, which provided digital gene expression profiles for each individual cell.

Nucleic acid extraction and melanin removal

Total RNA and genomic DNA was co-purified from frozen melanoma specimens using the Qiagen AllPrep DNA/RNA Kit (catalog #80204) according to the manufacturer's instructions (Qiagen, Valencia, CA). Melanin, a known inhibitor of enzymatic reactions, coprecipitates with the RNA, therefore the RNA required further purification. Performed as described by Gariboldi, modified with a RNeasy (Qiagen, Valencia, CA) column-based clean-up to remove the additives used to bind the melanin45.

IDT exome

A 700 ng aliquot of the existing WGS library was used for the exome capture. Five libraries were pooled at an equimolar ratio yielding a ~ 3.5 µg library pool prior to the hybrid capture. The library pools were hybridized with the xGen Exome Research Panel v1.0 reagent (IDT Technologies) that spans a 39 Mb target region (19,396 genes) of the human genome. The concentration of each captured library pool was accurately determined through qPCR (Kapa Biosystems) to produce cluster counts appropriate for the NovaSeq6000 platform (Illumina). 2 × 15 bp sequence data was generated ~ 50 Gb per library targeting a mean depth of coverage of 500x.

RNA-seq

Total RNA was isolated from ~ 700 K cells utilizing the AllPrep DNA extraction kit (Qiagen). ERCC RNA Spike-In Mix 1 was added to 100-250 ng of total RNA as outlined by the manufacturer (Ambion, Life Technologies). The ERCC control mix is a set of external RNA controls that enable performance assessment for gene expression experiments. The cDNA library was prepared with the TruSeq Stranded Total RNA Sample Prep with Ribo-Zero Gold kit (Illumina). The concentration of each cDNA library was determined through qPCR (Kapa Biosystems). 2 × 150 reads were generated on the HiSeq4000/NovaSeq6000 instrument (Illumina) generating ~ 83 million read pairs/sample.

scRNA-seq data quantification preprocessing

For single cell RNA-seq analysis, we used Cell Ranger v2.1.1 from 10 × Genomics for de-multiplexing sequence data into FASTQ files, aligning reads to the human genome (GRCh38), and generating gene-by-cell UMI count matrix.

Seurat v3.0.046,47 was used for all subsequent analysis. First, a series of quality filters were applied to the data to remove those barcodes which fell into any one of these categories recommended by Seurat: too few total transcript counts (< 300); possible debris with too few genes expressed (< 200) and too few UMIs (< 1,000); possible more than one cell with too many genes expressed (> 50,000) and too many UMIs (> 10,000); possible dead cell or a sign of cellular stress and apoptosis with too high proportion of mitochondrial gene expression over the total transcript counts (> 20%). Finally, predicted doublets were also removed using scrublet V0.2.3.

We constructed a Seurat object using the unfiltered feature-barcode matrix for each sample. Each sample was scaled and normalized using Seurat’s ‘SCTransform’ function to correct for batch effects (with parameters: vars.to.regress = c(“nCount_RNA”, “percent.mito”), return.only.var.genes = F).

scRNA-seq cell type annotation

Cell types were assigned to each cluster by manually reviewing the expression of marker genes. The marker genes used were MLANA, SOX10, MITF, AXL, GPR137B, S100B (Melanoma cells) CD79A, CD79B, MS4A1 (B cells); CD8A, CD8B, CD7, CD3E (CD8 + T cells); CD4, IL7R, CD7, CD3E (CD4 + T cells); NKG7, GNLY (NK cells); MZB1, SDC1, IGHG1 (Plasma cells); FCGR3A (Macrophages); CD14, LYZ (Monocytes); FCER1A, CLEC10A (Dendritic cells); IL3RA, CLEC4C (Plasmacytoid Dendritic cells); AHSP1, HBA, HBB (Erythrocytes); FAP, PDGFRA, ACTA2, TND (Cancer-associated fibroblasts, CAFs) and PHGR1, FABP2, FABP1, ALDOB (Intestinal cell).

Somatic mutation detection

Somatic variants were called by our SomaticWrapper pipeline, which includes four established bioinformatic tools, namely Strelka, Mutect, VarScan2 (2.3.83), and Pindel (0.2.54)48–51. We retained SNVs and INDELs using the following strategy: keep SNVs called by any 2 callers among Mutect, VarScan, and Strelka and INDELs called by any 2 callers among VarScan, Strelka, and Pindel. For these merged SNVs and INDELs, we applied coverage cut-offs of 14X and 8X for tumor and normal, respectively. We also filtered SNVs and INDELs with a high-pass variant allele fraction (VAF) of 0.05 in tumor and a low-pass VAF of 0.02 in normal. The SomaticWrapper pipeline is freely available from GitHub at https://github.com/ding-lab/somaticwrapper.

Copy number detection

We used CNVkit (v0.9.4)52 to measure copy number alterations (CNAs) in matched tumor-normal samples. Low quality CNAs were filtered based on the coverage (< 20), the number of probes (< 10), and length (< 5 kb). To define absolute copy numbers from CNVkit, the threshold is as follows: −t −1.3, −0.4, 0.3, 0.9. Deletion, loss, neutral, gain, and amplification of segment or gene-level defined as 0, 1, 2, 3, > 5 in absolute copy number.

Analysis of bulk RNA-seq data

Transcript quantification was performed using Kallisto v0.44.053, against the GENCODE transcript reference (release 29, GRCh38). Subsequent analysis was performed using R v3.6.0. R package ‘tximport’ v1.12.054 was used to import and aggregate transcript level data to the gene level.

scRNA-seq data integration

We used the “merge” function in Seurat to combine the Seurat objects from multiple samples after quality control. Cell types were assigned based on manual review of marker gene expression (as described above). Cells with inconsistent cell type assignments between the integrated and individual analyses were filtered out. In some cases, the inconsistencies arose from evident clustering issues (for example, when reviewing marker gene expression, two sub-clusters were obvious within one cluster). Such instances were manually resolved and the cells were rescued.

Single cell CNV detection and clonality analysis

To detect large-scale chromosomal copy number variations using single-cell RNA-seq data, we used inferCNV v1.2.0 (https://github.com/broadinstitute/inferCNV) to obtain relative expression intensity of melanoma cells in comparison to a set of reference “normal” cells, including B cells, T cells, Erythrocytes, NK cells, etc. To investigate clonality within each patient tumor, cutoff = 0.1, analysis mode “subcluster” and HMM was used for revealing cnv signals. Arm-level CNV was determined by converting from gene-level CNV based on GRCh38 cytoband information (http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/). Finally, we visualized the subclonal CNV evolutions within each patient tumor by UPhyloplot2 (https://github.com/harbourlab/UPhyloplot2). We used −c 10 to remove subclones comprising < 5% of cells for plotting.

Pseudotime-based trajectory analysis

Trajectory analysis was performed by Monocle 255. Taking patient 1451 as an example, melanoma cells from primary tumor and metastatic tumor were extracted and imported into Monocle 2. Parameters for the analysis were consistent with the tutorial (http://cole-trapnell-lab.github.io/monocle-release/docs/#constructing-single-cell-trajectories), except that different samples were set as the variable for differential expression test. To choose genes for ordering, we set 1e-10 as the q value cut-off. Finally, we visualized melanoma subclusters projection in the trajectory with the function “plot_cell_trajectory”.

Differential expression analysis

Differential expression analysis was performed in order to compare tumor subpopulations, such as Fig. 3c, Fig. 4c-g and Supplementary Fig. 3f., using the default test (Wilcoxon Rank Sum test) of function FindMarkers (from the Seurat package) with the specified parameters: min.pct = 0.25, logfc.threshold = 0.25, and only.pos = T.

Pathway analysis

We performed pathway enrichment analysis using ReactomePA (available at: https://github.com/YuLab-SMU/ReactomePA) on the differentially expressed genes of each tumor sub-cluster in 762_M and 1199_M1 (Fig. 3a). pvalueCutoff = 0.05.

All methods were performed in accordance with the guidelines and regulations of under Washington University in St. Louis approved IRB research protocol, entitled Advanced Genetic and Molecular Analsysis of Solid Tumors.

Supplementary Information

Supplementary Figures.

Supplementary Table 1.

Supplementary Table 2.

Supplementary Information

The online version contains supplementary material available at 10.1038/s41598-024-72255-9.

Acknowledgements

melanoma patients, families, and professionals who have contributed to this study. R.C.F. and L.D. are supported by NCI U2CCA233303 and institutional funds. R.C.F. is supported by U54CA224083 and R01CA248277 as well. L.D. is also supported by U24CA211006 and R01HG009711. B.A.K. credits support from the Washington University School of Medicine Surgical Oncology Basic Science and Translational Research Training Program grant T32CA009621, from the National Cancer Institute (NCI).

Author contributions

R.C.F. and L.D. led project design. B.A.K. led experimental design, sample collection and processing. L.Y. led data analysis, interpretation and figure generation. B.A.K, Y.B., J.M., S.G. , C.W. and A.O. collected samples, performed treatment experiments and generated sequencing data. L.Y., S.S., A.W. and Q.G. developed data processing and analysis pipelines. L.Y. and B.A.K wrote the manuscript. M.A.W helped polish figures. M.W helped polish the manuscript. R.C.F and L.D. reviewed the manuscript.

Data availability

The sequence data generated in this study has been submitted to the NCBI BioProject database PRJNA742837 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA742837?reviewer = f949g4roaer33omqln329bhr7i).

Code availability

For mutation calling, code can be accessed at https://github.com/ding-lab/somaticwrapper

Competing interests

The authors declare no competing interests.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Lijun Yao and Bradley A. Krasnick.
==== Refs
References

1. Larkin J Combined vemurafenib and cobimetinib in BRAF-mutated melanoma N. Engl. J. Med. 2014 371 1867 1876 10.1056/NEJMoa1408868 25265494
Larkin, J. et al. Combined vemurafenib and cobimetinib in BRAF-mutated melanoma. N. Engl. J. Med. 371, 1867–1876 (2014).25265494
2. Long GV Adjuvant Dabrafenib plus Trametinib in stage III BRAF-mutated melanoma N. Engl. J. Med. 2017 377 1813 1823 10.1056/NEJMoa1708539 28891408
Long, G. V. et al. Adjuvant Dabrafenib plus Trametinib in stage III BRAF-mutated melanoma. N. Engl. J. Med. 377, 1813–1823 (2017).28891408
3. Robert C Five-year outcomes with Dabrafenib plus Trametinib in metastatic melanoma N. Engl. J. Med. 2019 381 626 636 10.1056/NEJMoa1904059 31166680
Robert, C. et al. Five-year outcomes with Dabrafenib plus Trametinib in metastatic melanoma. N. Engl. J. Med. 381, 626–636 (2019).31166680
4. Ojha R ER translocation of the MAPK pathway drives therapy resistance in BRAF-mutant melanoma Cancer Discov. 2019 9 396 415 10.1158/2159-8290.CD-18-0348 30563872
Ojha, R. et al. ER translocation of the MAPK pathway drives therapy resistance in BRAF-mutant melanoma. Cancer Discov. 9, 396–415 (2019).30563872
5. Johannessen CM A melanocyte lineage program confers resistance to MAP kinase pathway inhibition Nature 2013 504 138 142 10.1038/nature12688 24185007
Johannessen, C. M. et al. A melanocyte lineage program confers resistance to MAP kinase pathway inhibition. Nature 504, 138–142 (2013).24185007
6. Snyder A Genetic basis for clinical response to CTLA-4 blockade in melanoma N. Engl. J. Med. 2014 371 2189 2199 10.1056/NEJMoa1406498 25409260
Snyder, A. et al. Genetic basis for clinical response to CTLA-4 blockade in melanoma. N. Engl. J. Med. 371, 2189–2199 (2014).25409260
7. Kleffel S Melanoma cell-intrinsic PD-1 receptor functions promote tumor growth Cell 2015 162 1242 1256 10.1016/j.cell.2015.08.052 26359984
Kleffel, S. et al. Melanoma cell-intrinsic PD-1 receptor functions promote tumor growth. Cell 162, 1242–1256 (2015).26359984
8. Iwai Y Involvement of PD-L1 on tumor cells in the escape from host immune system and tumor immunotherapy by PD-L1 blockade Proc. Natl. Acad. Sci. U. S. A. 2002 99 12293 12297 10.1073/pnas.192461099 12218188
Iwai, Y. et al. Involvement of PD-L1 on tumor cells in the escape from host immune system and tumor immunotherapy by PD-L1 blockade. Proc. Natl. Acad. Sci. U. S. A. 99, 12293–12297 (2002).12218188
9. Network CGA Genomic classification of cutaneous melanoma Cell 2015 161 1681 1696 10.1016/j.cell.2015.05.044 26091043
Network, C. G. A. Genomic classification of cutaneous melanoma. Cell 161, 1681–1696 (2015).26091043
10. Tirosh I Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq Science 2016 352 189 196 10.1126/science.aad0501 27124452
Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189–196 (2016).27124452
11. Lawrence MS Mutational heterogeneity in cancer and the search for new cancer-associated genes Nature 2013 499 214 218 10.1038/nature12213 23770567
Lawrence, M. S. et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature 499, 214–218 (2013).23770567
12. Hodis E A landscape of driver mutations in melanoma Cell 2012 150 251 263 10.1016/j.cell.2012.06.024 22817889
Hodis, E. et al. A landscape of driver mutations in melanoma. Cell 150, 251–263 (2012).22817889
13. Anna B Mechanism of UV-related carcinogenesis and its contribution to nevi/melanoma Exp. Rev. Dermatol. 2007 2 451 469 10.1586/17469872.2.4.451
Anna, B. et al. Mechanism of UV-related carcinogenesis and its contribution to nevi/melanoma. Exp. Rev. Dermatol. 2, 451–469 (2007).
14. Vincent KM Postovit L-M Investigating the utility of human melanoma cell lines as tumour models Oncotarget 2017 8 10498 10509 10.18632/oncotarget.14443 28060736
Vincent, K. M. & Postovit, L.-M. Investigating the utility of human melanoma cell lines as tumour models. Oncotarget 8, 10498–10509 (2017).28060736
15. Luebker SA Zhang W Koepsell SA Comparing the genomes of cutaneous melanoma tumors to commercially available cell lines Oncotarget 2017 8 114877 114893 10.18632/oncotarget.22928 29383127
Luebker, S. A., Zhang, W. & Koepsell, S. A. Comparing the genomes of cutaneous melanoma tumors to commercially available cell lines. Oncotarget 8, 114877–114893 (2017).29383127
16. Jerby-Arnon L A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade Cell 2018 175 984 997.e24 10.1016/j.cell.2018.09.006 30388455
Jerby-Arnon, L. et al. A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade. Cell 175, 984-997.e24 (2018).30388455
17. Sade-Feldman M Defining T cell states associated with response to checkpoint immunotherapy in melanoma Cell 2019 176 404 10.1016/j.cell.2018.12.034 30633907
Sade-Feldman, M. et al. Defining T cell states associated with response to checkpoint immunotherapy in melanoma. Cell 176, 404 (2019).30633907
18. El Marsafy S Bagot M Bensussan A Mauviel A Dendritic cells in the skin–potential use for melanoma treatment Pigment Cell Melanoma Res. 2009 22 30 41 10.1111/j.1755-148X.2008.00532.x 19040502
El Marsafy, S., Bagot, M., Bensussan, A. & Mauviel, A. Dendritic cells in the skin–potential use for melanoma treatment. Pigment Cell Melanoma Res. 22, 30–41 (2009).19040502
19. North JP Vemula SS Bastian BC Chromosomal copy number analysis in melanoma diagnostics Methods Mol. Biol. 2014 1102 199 226 10.1007/978-1-62703-727-3_12 24258981
North, J. P., Vemula, S. S. & Bastian, B. C. Chromosomal copy number analysis in melanoma diagnostics. Methods Mol. Biol. 1102, 199–226 (2014).24258981
20. McPherson LA Loktev AV Weigel RJ Tumor suppressor activity of AP2alpha mediated through a direct interaction with p53 J. Biol. Chem. 2002 277 45028 45033 10.1074/jbc.M208924200 12226108
McPherson, L. A., Loktev, A. V. & Weigel, R. J. Tumor suppressor activity of AP2alpha mediated through a direct interaction with p53. J. Biol. Chem. 277, 45028–45033 (2002).12226108
21. Sun Y Long J Zhou Y Angiopoietin-like 4 promotes melanoma cell invasion and survival through aldolase A Oncol. Lett. 2014 8 211 217 10.3892/ol.2014.2071 24959248
Sun, Y., Long, J. & Zhou, Y. Angiopoietin-like 4 promotes melanoma cell invasion and survival through aldolase A. Oncol. Lett. 8, 211–217 (2014).24959248
22. Liu R Silencing of PKG1 expression enhances the efficacy of vemurafenib against melanoma cell Zhongguo Ying Yong Sheng Li Xue Za Zhi 2017 33 289 293 29926631
Liu, R. et al. Silencing of PKG1 expression enhances the efficacy of vemurafenib against melanoma cell. Zhongguo Ying Yong Sheng Li Xue Za Zhi 33, 289–293 (2017).29926631
23. Andrzejewski S Siegel PM St-Pierre J Metabolic profiles associated with metformin efficacy in cancer Front. Endocrinol. 2018 9 372 10.3389/fendo.2018.00372
Andrzejewski, S., Siegel, P. M. & St-Pierre, J. Metabolic profiles associated with metformin efficacy in cancer. Front. Endocrinol. 9, 372 (2018).
24. Jaune E Rocchi S Metformin: Focus on melanoma Front. Endocrinol. 2018 9 472 10.3389/fendo.2018.00472
Jaune, E. & Rocchi, S. Metformin: Focus on melanoma. Front. Endocrinol. 9, 472 (2018).
25. Boudhraa Z Annexin A1 in primary tumors promotes melanoma dissemination Clin. Exp. Metastasis 2014 31 749 760 10.1007/s10585-014-9665-2 24997993
Boudhraa, Z. et al. Annexin A1 in primary tumors promotes melanoma dissemination. Clin. Exp. Metastasis 31, 749–760 (2014).24997993
26. Rohwer N Annexin A1 sustains tumor metabolism and cellular proliferation upon stable loss of HIF1A Oncotarget 2016 7 6693 6710 10.18632/oncotarget.6793 26760764
Rohwer, N. et al. Annexin A1 sustains tumor metabolism and cellular proliferation upon stable loss of HIF1A. Oncotarget 7, 6693–6710 (2016).26760764
27. Zhu G-H Dai H-P Shen Q Zhang Q Downregulation of LPXN expression by siRNA decreases the malignant proliferation and transmembrane invasion of SHI-1 cells Oncol. Lett. 2019 17 135 140 30655748
Zhu, G.-H., Dai, H.-P., Shen, Q. & Zhang, Q. Downregulation of LPXN expression by siRNA decreases the malignant proliferation and transmembrane invasion of SHI-1 cells. Oncol. Lett. 17, 135–140 (2019).30655748
28. Lin X C-myc overexpression drives melanoma metastasis by promoting vasculogenic mimicry via c-myc/snail/Bax signaling J. Mol. Med. 2017 95 53 67 10.1007/s00109-016-1452-x 27543492
Lin, X. et al. C-myc overexpression drives melanoma metastasis by promoting vasculogenic mimicry via c-myc/snail/Bax signaling. J. Mol. Med. 95, 53–67 (2017).27543492
29. Rohde M Members of the heat-shock protein 70 family promote cancer cell growth by distinct mechanisms Genes Dev. 2005 19 570 582 10.1101/gad.305405 15741319
Rohde, M. et al. Members of the heat-shock protein 70 family promote cancer cell growth by distinct mechanisms. Genes Dev. 19, 570–582 (2005).15741319
30. Budina-Kolomets A HSP70 inhibition limits FAK-dependent invasion and enhances the response to melanoma treatment with BRAF inhibitors Cancer Res. 2016 76 2720 2730 10.1158/0008-5472.CAN-15-2137 26984758
Budina-Kolomets, A. et al. HSP70 inhibition limits FAK-dependent invasion and enhances the response to melanoma treatment with BRAF inhibitors. Cancer Res. 76, 2720–2730 (2016).26984758
31. Park S-Y DNAJB1 negatively regulates MIG6 to promote epidermal growth factor receptor signaling Biochim. Biophys. Acta 2015 1853 2722 2730 10.1016/j.bbamcr.2015.07.024 26239118
Park, S.-Y. et al. DNAJB1 negatively regulates MIG6 to promote epidermal growth factor receptor signaling. Biochim. Biophys. Acta 1853, 2722–2730 (2015).26239118
32. Yu T The effect of tumor purity on next generation sequencing of colorectal cancer J. Clin. Orthod. 2022 40 e15557 e15557
Yu, T. et al. The effect of tumor purity on next generation sequencing of colorectal cancer. J. Clin. Orthod. 40, e15557–e15557 (2022).
33. Birkeland E Patterns of genomic evolution in advanced melanoma Nat. Commun. 2018 9 2665 10.1038/s41467-018-05063-1 29991680
Birkeland, E. et al. Patterns of genomic evolution in advanced melanoma. Nat. Commun. 9, 2665 (2018).29991680
34. Larkin J Five-year survival with combined Nivolumab and Ipilimumab in advanced melanoma N. Engl. J. Med. 2019 381 1535 1546 10.1056/NEJMoa1910836 31562797
Larkin, J. et al. Five-year survival with combined Nivolumab and Ipilimumab in advanced melanoma. N. Engl. J. Med. 381, 1535–1546 (2019).31562797
35. Rossi A Drug resistance of BRAF-mutant melanoma: Review of up-to-date mechanisms of action and promising targeted agents Eur. J. Pharmacol. 2019 862 172621 10.1016/j.ejphar.2019.172621 31446019
Rossi, A. et al. Drug resistance of BRAF-mutant melanoma: Review of up-to-date mechanisms of action and promising targeted agents. Eur. J. Pharmacol. 862, 172621 (2019).31446019
36. Chen H-L Effect of metformin on proliferation capacity, apoptosis and glycolysis in K562 cells Zhongguo Shi Yan Xue Ye Xue Za Zhi 2019 27 1387 1394 31607288
Chen, H.-L. et al. Effect of metformin on proliferation capacity, apoptosis and glycolysis in K562 cells. Zhongguo Shi Yan Xue Ye Xue Za Zhi 27, 1387–1394 (2019).31607288
37. Morales DR Morris AD Metformin in cancer treatment and prevention Annu. Rev. Med. 2015 66 17 29 10.1146/annurev-med-062613-093128 25386929
Morales, D. R. & Morris, A. D. Metformin in cancer treatment and prevention. Annu. Rev. Med. 66, 17–29 (2015).25386929
38. Niehr F Combination therapy with vemurafenib (PLX4032/RG7204) and metformin in melanoma cell lines with distinct driver mutations J. Transl. Med. 2011 9 76 10.1186/1479-5876-9-76 21609436
Niehr, F. et al. Combination therapy with vemurafenib (PLX4032/RG7204) and metformin in melanoma cell lines with distinct driver mutations. J. Transl. Med. 9, 76 (2011).21609436
39. Ryabaya O Metformin increases antitumor activity of MEK inhibitor binimetinib in 2D and 3D models of human metastatic melanoma cells Biomed. Pharmacother. 2019 109 2548 2560 10.1016/j.biopha.2018.11.109 30551515
Ryabaya, O. et al. Metformin increases antitumor activity of MEK inhibitor binimetinib in 2D and 3D models of human metastatic melanoma cells. Biomed. Pharmacother. 109, 2548–2560 (2019).30551515
40. Vujic I Metformin and trametinib have synergistic effects on cell viability and tumor growth in NRAS mutant cancer Oncotarget 2015 6 969 978 10.18632/oncotarget.2824 25504439
Vujic, I. et al. Metformin and trametinib have synergistic effects on cell viability and tumor growth in NRAS mutant cancer. Oncotarget 6, 969–978 (2015).25504439
41. Martin MJ Hayward R Viros A Marais R Metformin accelerates the growth of BRAF V600E-driven melanoma by upregulating VEGF-A Cancer Discov. 2012 2 344 355 10.1158/2159-8290.CD-11-0280 22576211
Martin, M. J., Hayward, R., Viros, A. & Marais, R. Metformin accelerates the growth of BRAF V600E-driven melanoma by upregulating VEGF-A. Cancer Discov. 2, 344–355 (2012).22576211
42. Alexandrov LB Signatures of mutational processes in human cancer Nature 2013 500 415 421 10.1038/nature12477 23945592
Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).23945592
43. Morita K Clonal evolution of acute myeloid leukemia revealed by high-throughput single-cell genomics Nat. Commun. 2020 11 5327 10.1038/s41467-020-19119-8 33087716
Morita, K. et al. Clonal evolution of acute myeloid leukemia revealed by high-throughput single-cell genomics. Nat. Commun. 11, 5327 (2020).33087716
44. Petti AA A general approach for detecting expressed mutations in AML cells using single cell RNA-sequencing Nat. Commun. 2019 10 3660 10.1038/s41467-019-11591-1 31413257
Petti, A. A. et al. A general approach for detecting expressed mutations in AML cells using single cell RNA-sequencing. Nat. Commun. 10, 3660 (2019).31413257
45. Lagonigro MS CTAB-urea method purifies RNA from melanin for cDNA microarray analysis Pigment Cell Res. 2004 17 312 315 10.1111/j.1600-0749.2004.00155.x 15140079
Lagonigro, M. S. et al. CTAB-urea method purifies RNA from melanin for cDNA microarray analysis. Pigment Cell Res. 17, 312–315 (2004).15140079
46. Butler A Hoffman P Smibert P Papalexi E Satija R Integrating single-cell transcriptomic data across different conditions, technologies, and species Nat. Biotechnol. 2018 36 411 420 10.1038/nbt.4096 29608179
Butler, A., Hoffman, P., Smibert, P., Papalexi, E. & Satija, R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).29608179
47. Hafemeister C Satija R Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression Genome Biol. 2019 20 296 10.1186/s13059-019-1874-1 31870423
Hafemeister, C. & Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296 (2019).31870423
48. Cibulskis K Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples Nat. Biotechnol. 2013 31 213 219 10.1038/nbt.2514 23396013
Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219 (2013).23396013
49. Koboldt DC VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing Genome Res. 2012 22 568 576 10.1101/gr.129684.111 22300766
Koboldt, D. C. et al. VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).22300766
50. Saunders CT Strelka: Accurate somatic small-variant calling from sequenced tumor-normal sample pairs Bioinformatics 2012 28 1811 1817 10.1093/bioinformatics/bts271 22581179
Saunders, C. T. et al. Strelka: Accurate somatic small-variant calling from sequenced tumor-normal sample pairs. Bioinformatics 28, 1811–1817 (2012).22581179
51. Ye K Schulz MH Long Q Apweiler R Ning Z Pindel: A pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads Bioinformatics 2009 25 2865 2871 10.1093/bioinformatics/btp394 19561018
Ye, K., Schulz, M. H., Long, Q., Apweiler, R. & Ning, Z. Pindel: A pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads. Bioinformatics 25, 2865–2871 (2009).19561018
52. Talevich E Shain AH Botton T Bastian BC CNVkit: Genome-wide copy number detection and visualization from targeted DNA sequencing PLoS Comput. Biol. 2016 12 e1004873 10.1371/journal.pcbi.1004873 27100738
Talevich, E., Shain, A. H., Botton, T. & Bastian, B. C. CNVkit: Genome-wide copy number detection and visualization from targeted DNA sequencing. PLoS Comput. Biol. 12, e1004873 (2016).27100738
53. Bray NL Pimentel H Melsted P Pachter L Near-optimal probabilistic RNA-seq quantification Nat. Biotechnol. 2016 34 525 527 10.1038/nbt.3519 27043002
Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34, 525–527 (2016).27043002
54. Soneson C Love MI Robinson MD Differential analyses for RNA-seq: Transcript-level estimates improve gene-level inferences F1000Res 2015 4 1521 10.12688/f1000research.7563.1 26925227
Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: Transcript-level estimates improve gene-level inferences. F1000Res 4, 1521 (2015).26925227
55. Qiu X Reversed graph embedding resolves complex single-cell trajectories Nat. Methods 2017 14 979 982 10.1038/nmeth.4402 28825705
Qiu, X. et al. Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14, 979–982 (2017).28825705
