Clonal dynamics shaped by diverse drug-tolerant persister states in melanoma resistance

Abstract

Most advanced melanomas initially respond to targeted therapy but eventually relapse. Rather than acquiring new mutations, resistance is driven by drug-tolerant persister cells that enter a reversible drug-refractory state. We developed MeRLin, a high-resolution lineage tracing platform integrating cellular barcoding, single-cell transcriptomics, RNA fluorescence in situ hybridization (FISH), and computational analyses to track clonal and transcriptional dynamics in patient-derived melanoma models during prolonged therapy. Clonal dynamics revealed that persister subpopulations first responded to treatment but persisted and expanded during minimal residual disease, ultimately leading to tumor recurrence. Pre-treatment melanoma populations diversified into four conserved persister states characterized by stress-like, lipid metabolism, PI3K signaling, and extracellular matrix remodeling programs associated with adaptive resistance. Spatial transcriptomics showed the organization of these adaptive programs and a complex signaling network of autocrine and paracrine interactions among persister subpopulations. Barcoded RNA-FISH enabled spatial mapping of clonal identity and gene expression, revealing in situ co-localization of a dominant resistant clone with SLC2A1 expression. MeRLin provides a robust framework for dissecting cancer heterogeneity and identifying vulnerabilities in persister populations.

Fig. 1.. Distinct clone fates are associated with variable expression programs.

a, Schematic representation of the MeRLin construct and barcoded WM4237–1 cells exposed to targeted therapy. The vector encoded firefly luciferase and fluorescent protein mNeptune2.5. Lineage tracing is enabled by semi-random barcodes (BC) transcribed within the 3’-untranslated region of mRNAs. WM4237–1 cells were transduced with the MeRLin barcode library. EF1α, elongation factor 1α; T2A, self-cleaving peptide; WPRE, woodchuck hepatitis virus post-transcriptional regulatory element; PolyA, polyadenylation signal. b, Uniform manifold approximation and projection (UMAP) to visualize the separation of BRAFi/MEKi-treated (red) and control (CTRL, blue) WM4237–1 cells using scRNA-seq data. c, Proportion of cells in each cell cycle phase inferred from scRNA-seq. P < 2.2 × 10⁻¹⁶, two-tailed Fisher's exact test. d, Shannon diversity indices (CTRL vs. BRAFi/MEKi-treated cells). P < 2.2 × 10⁻¹⁶, two-tailed Hutcheson’s t-test. e, Lineage tracing of BRAFi/MEKi-treated cells to pre-treatment ancestor cells (CTRL), quantified by the percentage of cells carrying each specific barcode. f, UMAP showing sensitive clones (blue) and resistant clones (red) following treatment. g, Resistant large-sized, resistant small-sized and sensitive clones showed significance in cellular expression programs. And in differentially expressed genes h, *P < 0.05, **P < 0.01, ***P < 0.001; two-tailed Wilcoxon rank-sum test.

Fig. 2.. Clonal dynamics reveal the emergence of adaptive resistance.

a, Experimental outline: (1) Barcoded WM4237–1 cells were injected subcutaneously into NSG mice. (2) BRAFi/MEKi treatment and development of resistance. (3) Tumors were excised and processed at multiple time points. (4) Transcribed barcodes were sequenced and visualized in situ. Pie charts at the bottom illustrate the changes in barcode distribution during treatment. b, Growth curve of barcoded WM4237–1 tumors before and after BRAFi/MEKi treatment. Tumor collection time points: Day 0 (pre-treatment), n = 3; Day 21 (early MRD), n = 3; Day 57 (pre-recurrence), n = 3; and endpoint (resistant), n = 6. c, Normalized abundance of individual barcodes revealed distinct clonal trajectories: (1) persister subpopulations that initially regressed but later expanded under treatment; (2) sensitive subpopulations eliminated by BRAFi/MEKi; and (3) multi-fate subpopulations that fluctuated in abundance without dominating the resistant tumors. d, Stacked plot illustrating clonal dynamics driving tumor growth during BRAFi/MEKi treatment. e, Top panel: Pie chart illustrating the hierarchical clonal composition of endpoint tumors from scRNA-seq data. The 15-bp 3′ end sequences of the five most dominant barcodes are shown on the right, while the proportions of barcodes ranked 6–10, 11–20, minor populations, and singletons were grouped respectively; Bottom panel: normalized abundance of top-ranked barcodes (ranks 1–3) across prolonged treatment. f, UMAP to visualize single-cell transcriptome clusters. Day 0 (pre-treatment), Day 21 (early MRD), and Endpoint (resistant). g, Proportion of cells in each cell cycle phase inferred from scRNA-seq. ***P < 0.001; two-tailed Fisher's exact test. h, UMAP visualization of subpopulations from the endpoint tumor including persister (red), sensitive (green), and multi-fate (blue) cells. i, Proportions of each endpoint subpopulation escaping the non-proliferative MRD state and re-entering a proliferative state resembling the pre-treatment tumor.

Fig. 3.. DTP cells acquire resistance via distinct transcriptional programs.

a, ClonoCluster was used to incorporate clonal barcode information (colored) into the transcriptome (grey). UMAPs show the effect of increasing Warp Factor values (0, 6, and 10) on the structure of endpoint tumor scRNA-seq data (Methods). b, UMAP visualizations of hybrid clusters (left) integrating clonal and transcriptomic information and their clonal fate (right). Barcode groups 1–4 (red, blue, green, brown in the left panel) represent persister subpopulations (red in the right panel), whereas barcode group 5 (gray in the left panel) corresponds to sensitive (green in the right) and multi-fate (blue in the right) subpopulations. n, unique barcode number in each barcode group or clonal fate. c, Proportion of cells in each cell cycle phase inferred from scRNA-seq. *** P < 0.001; two-tailed Fisher's exact test. d, Cells on the UMAP recolored by the expression of the melanoma differentiation markers MITF and TYRP1. e, Functional enrichment terms of barcode groups 1–4 identified in b. P-values were determined using two-tailed Fisher’s exact tests (Table S8). f, The AUCell scores (color scale) of the top functionally enriched programs per state projected on UMAP (Table S7). g, Discriminative marker genes (n = 5–7) for each barcode group (Table S6). Fold change ≥ 2, cell percentage ≥ 50%, and FDR < 0.05. h, UMAPs showing SERPINE2, DUSP4, and FXYD3 expression within the barcoded WM4237–1 persister subpopulations.

Fig. 4.. Regulatory mechanisms underlying adaptive resistance.

a, SCENIC analysis identifying key regulators of each barcode group. b, Regulon activity of ETS1 on its target genes in the ECM remodeling state. c, Survival analysis of TCGA data revealed that high ETS1 expression was associated with worse patient outcomes (P = 0.018). d, Volcano plots showing depleted sgRNAs from the CRISPR screen, highlighting the top negatively selected genes. Log₂FC (fold change) < − 1, and 18 genes were identified with FDR < 0.15. e, Alternative polyadenylation (APA) analysis showing 3′UTR lengthening and shortening events across in vitro treatment and in vivo time points, compared to control cells and pre-treatment tumors. ***P < 0.001; two-tailed Fisher's exact test. f, Copy number variation (CNV) scores for each barcode group (left). ***P < 0.001; two-tailed t-test. Spearman's rank correlation coefficients (rs) between CNV scores and gene expression in each barcode group (right). rs₁ = 0.21, P₁ = 0.047; rs₃ = 0.28, P₃ = 0.012 (median rs, subscripts indicate barcode groups). g, Copy number variation (CNV) scores of BRAF in each barcode group (left). Spearman's rank correlation coefficients between CNV and BRAF expression: rs₁ = 0.44, rs₂ = 0.56, rs₃ = 0.65, rs₄ = 0.23, rs₅ = 0.29, (subscripts indicate barcode groups), ***P < 0.001. UMAP (right) shows BRAF expression. h, Copy number variation of CCND1 in each barcode group (left). Spearman's rank correlation coefficients between gene expression and CNV of CCND1 in each barcode group: rs₁ = 0.23, rs₂ = 0.1, rs₃ = 0.17, rs₄ = 0.13 (subscripts indicate barcode groups). *P < 0.05, **P < 0.01, ***P < 0.001. UMAP (right) shows CCND expression.

Fig. 5.. Persister states occur through treatment and across models.

a, ClonoCluster (top panel) applied to scRNA-seq data from day 21 barcoded WM4237–1 tumors. UMAPs show barcode groups 1–5 (red, blue, green, brown, gray) from endpoint tumors; UMAP (bottom panel) illustratesing upregulated TYRP1 expression in barcode group 1 in the day 21 barcoded WM4237–1 tumor. b, Pathway enrichment scores derived from bulk RNA-seq of barcoded WM4237–1 tumors across BRAFi/MEKi treatment. *P < 0.05, **P < 0.01, ***P < 0.001; two-tailed t-test. c, UMAPs showing persister states from non-barcoded WM4007 scRNA-seq at endpoint and day 13 MRD stage (d). e, UMAPs illustrating persister signatures from non-barcoded WM4380–2 scRNA-seq data at the endpoint. f, Proportions of persister states across BRAF V600E-mutant PDX models WM4007 (day 13 and endpoint), WM4237–1, and WM4380–2. g, Survival analysis of TCGA data revealed that high expression of the ECM remodeling signature was associated with worse patient outcomes (P = 0.011).

Fig. 6.. Spatial profiling of diverse drug-tolerant persister states.

a, Spatial transcriptomics showing persister states overlaid on H&E-stained sections of a recurrent WM4237–1 PDX tumor. Stress-like (red), lipid metabolism (blue), PI3K signaling (green), ECM remodeling (brown), and melanocytic (gray) states. b, Distinct spatial distribution of each persister state with quantification scale. c, Moran's indices showing spatial autocorrelation for each persister state. Stress-like (red), I = 0.27; lipid metabolism (blue), I = 0.31; PI3K signaling (green), I = 0.23; ECM remodeling (brown), I = 0.083; and melanocytic (grey), I = 0.21. **P < 0.01, ***P < 0.001. d, Chord diagrams reveal the strength of cell-cell communication among persister states. Stress-like (yellow) and lipid metabolism (blue) states emitted TGFβ signals (left); lipid metabolism (blue) and PI3K signaling (purple) states exhibited autocrine and paracrine BMP signaling (middle); and ECM remodeling state (red) sent collagen signals (right). e, Specific ligand-receptor pairs involved in TGFβ signaling in d (left). Outgoing TGFβ signaling from the stress-like state (top) and lipid metabolism state (bottom). f, UMAPs illustrating the expression of the dominant barcode GTTGAACGACCACAA-3 and the stress-like marker SLC2A1. g, RNA-FISH (top panel) of a recurrent WM4237–1 barcoded tumor revealed co-localization of a resistant clone (barcode: GTTGAACGACCACAA-3) and the stress-like marker SLC2A1. Statistical analysis (bottom panel) showed significant co-localization between SLC2A1 and the resistant clone across all patch sizes (Fisher’s exact test, P < 0.014), compared to co-localization with a non-specific probe set (patch size of 10 × 10, Fisher's exact test, P = 0.089). The analysis was based on three ROIs (regions of interest) each.