Plasticity of Extrachromosomal and Intrachromosomal BRAF Amplifications in Overcoming Targeted Therapy Dosage Challenges

Genomic instability confers cancer cells with a growing list of hallmarks such as enhanced invasion and deregulated cellular energetics (1). Among many types of instability-driven mutations, focal amplifications (FA) of oncogenes in cancer genomes are a major contributor of neoplastic progression and therapeutic resistance (2–5). There are primarily two modes of FA structural topology: double minute (DM) and homogeneously staining region (HSR). DMs are circular extrachromosomal DNAs (ecDNA) that allow copies of oncogenes to exist freely in nuclei and retain intact, altered, or even elevated transcription activity due to high chromatin accessibility, enhancer hijacking, and formation of transcription hubs (6–12). DMs are able to replicate autonomously, but are acentric and therefore segregate into daughter cells randomly (13–16). HSRs are intrachromosomal amplifications resulting in long segments with uniform staining intensities in cytogenetics (17). Several models regarding the generation and interchange of these two kinds of FAs have been proposed, including but not restricted to episomal, chromo-thripsis, breakage–fusion–bridge, and integration mechanisms (18, 19, 20–29). The high prevalence of both kinds of FAs supports their importance in tumorigenesis (15). DMs have been observed in a large number of tumors of different types, especially in glioblastomas [∼55% by whole-genome sequencing (WGS) inferred ecDNA; ref. 30] and neuroblastomas (∼31% by cytogenetics; ref. 31), but rarely in normal tissues. A high occurrence of the HSR form of FAs is found in particular cancer types such as squamous cell carcinoma (12.1% by cytogenetics), but across all cancers, HSRs have a slightly lower frequency compared with DMs (31).

Mutations in BRAF, a serine/threonine RAF family kinase and a key upstream member of the MAPK pathway, have been associated with many cancer types. The frequency of BRAF mutations varies widely across cancer types. For example, BRAF mutations are relatively common in thyroid cancer and skin melanoma (60% and 52%, respectively), but are very rare in kidney cancers (0.3%), based on The Cancer Genome Atlas (TCGA) database. In melanoma therapy, the development of inhibitors targeting the BRAF^V600E mutation, such as vemurafenib (VEM) and dabrafenib as well as combinatorial treatments with other MAPK pathway inhibitors (MAPKi), has greatly improved patient survival (32). However, acquired resistance often compromises the efficacy of these therapies. To date, many resistance mechanisms to BRAF inhibition emerging during clinical treatment have been identified, including reactivation of the MAPK pathway, activation of the PI3K/AKT pathway, or both. This can occur via genomic mutations, genomic rearrangements such as kinase domain duplication, altered splice isoform variant expression, and cellular dedifferentiation (5, 33–41). One mechanism of reactivating the MAPK pathway that is frequently found in melanoma patient tumors is the acquisition of BRAF amplifications (5). As previously noted, these amplifications can be mediated through both DM and HSR FA modes (35, 42, 43). However, the details related to the generation, structure, dynamics, plasticities, and vulnerabilities of MAPK FAs due to acquired drug resistance in melanoma are incomplete, and as such are the focus of our current study.

In this study, through acquired BRAF and MEK inhibitor resistance, we developed a melanoma model system that dynamically harbors mutant BRAF^V600E in the form of DMs, HSRs, or both. Using single cell–derived clones, we found that increasing and/or decreasing kinase inhibitor dosage is a reproducible modulator of the number of DMs, the length of HSRs, the transition between these FA modes, and coupling with additional genomic rearrangements such as kinase domain duplications and alternative splicing. More-over, we observed plasticity of FAs involving other amplified MAPK genes such as RAF1 and NRAS in NRAS^MUT melanoma. Using optical mapping (OM) and WGS, we profiled the BRAF FA structures and found conserved amplicon boundaries between the DM and HSR modes. Furthermore, the observed junction sequences yielded initial insight into the mechanisms of integration and HSR shortening. In investigating the cellular liabilities of BRAF amplification, we identified an increased sensitivity to ferroptosis via GPX4 inhibition, which extends the spectrum of melanoma resistance–derived ferroptosis sensitivity beyond cases of dedifferentiation. Collectively, our findings on BRAF amplicon structure, DM and HSR plasticity, and potential vulnerabilities associated with BRAF FA-driven resistance highlight key therapeutic challenges and opportunities.

In order to generate an FA-positive melanoma model, we treated a BRAF^V600E human melanoma cell line, M249, with VEM [BRAF inhibitor (BRAFi)] and selumetinib [SEL; MEK inhibitor (MEKi)] to develop resistance (abbreviated as M249-VSR for vemurafenib and selumetinib resistant), as previously described (ref. 41; Fig. 1A). Upon the establishment of cells resistant at 2 μmol/L, FISH was performed and showed a high amplification of BRAF, primarily in the DM/ecDNA form. However, over the course of a few months in culture, these cells spontaneously switched their karyotypes to DM-negative and HSR-positive with a small number of exceptions (Fig. 1B; see Supplementary Fig. S1A–S1C and Methods section for categories and images of BRAF FA FISH-based karyotypes). To quantify the extent of FA, we performed quantitative real-time PCR (qPCR) on M249-VSR-DM and M249-VSR-HSR cells and found that there were 30- to 40-fold increases in BRAF copy number compared with M249 parental (M249-P) cells that contained five copies of BRAF (Fig. 1B and C). These amplifications also led to high protein levels of BRAF (Fig. 1D).

As qPCR is limited to investigating a small DNA region, we used WGS and comparative genomic hybridization (CGH) for M249-P and M249-VSR cells to reveal the full copy-number alterations/variations (CNA) across the genome (Fig. 1E; Supplementary Fig. S2A–S2C). Though there were other alterations, the most striking change upon acquisition of resistance was an FA of size ∼1.62 Mb at chr7q34, the region of the BRAF locus, with a fold increase consistent with qPCR results. The amplicon had highly similar start and end points in both the DM and HSR modes of amplification, supporting that the ensuing HSRs were generated through integration of DMs. Genes adjacent to BRAF on the amplicon were amplified to a similar degree (Fig. 1F), and the transcripts of these genes were also elevated as measured by RNA sequencing (RNA-seq; Fig. 1G). Such coamplifications have also been found on amplicons containing other oncogenes, e.g., MYC and EGFR (11, 12, 44). RNA-seq based single-nucleotide variant calling of DM and HSR M249 cell lines indicated that the BRAF 1799T>A (V600E) mutation was selected during FA development, with both DM and HSR cases displaying greater than 99% major allele frequency compared with 71% in the parental line (Fig. 1H).

We next characterized the structure of DM and HSR amplicons, aided by the inclusion of OM data. The observed OM junctions confirmed the circular structure (6, 30, 45) of the DMs/ecDNAs generated during the acquisition of resistance. In contrast, the parental M249 cells with five BRAF copies per cell show a linear arm-level amplification (Fig. 1I; Supplementary Fig. S2C). For HSRs, we investigated the sites of integration. Through cytogenetic G-banding, we found a limited level of heterogeneity of HSR integration sites, with integration on either chromosome 1 or 3, or on one or more marker (unidentifiable by G-banding) chromosomes (Fig. 1J).

To further dissect changes that occurred during the transition from DMs to HSRs, we isolated single cell–derived clones (SC) from the bulk M249-VSR population at an intermediary timepoint between the DM⁺ and HSR⁻ and DM⁻ and HSR⁺ karyotypes (Fig. 2A). Cultures derived from these clones were expanded and characterized for subsequent changes over a three-month timeframe. At the outset, three of the resultant clones had a DM⁺ and HSR⁻ karyotype (SC3, SC4, and SC401), one clone had a DM⁻ and HSR⁺ karyotype (SC2), and one had a DM⁺ and HSR⁺ karyotype (SC5; Fig. 2B and C; Supplementary Fig. S3A–S3D). Over the matching three-month time course, the bulk population began with a small percentage of DM⁻ and HSR⁺ cells, which gradually expanded to dominate the population (Fig. 2B and D, B1–B4). The SC experiments deconvoluted such changes by displaying a range of evolutionary trajectories that implicated de novo DM integration as HSRs, followed by selection of HSR⁺ cells.

First, although the majority of SC4 and SC401 cells kept their DM⁺ and HSR⁻ karyotype, some cells began to have the DM⁻ and HSR⁺ karyotype and some cells harbored both DMs and HSRs (Fig. 2B and C; Supplementary Fig. S3A and S3B). Second, HSR⁺ cases of some SC4 cells presented in a format that had three smaller HSR segments in different chromosomes in each cell (Fig. 2B and C), whereas only one long HSR, with or without an additional short HSR segment, was observed in M249-VSR-HSR bulk cells. These single-cell clone dynamics support de novo integration of DMs as HSRs. A previous study reported that nonhomologous end joining (NHEJ) is needed for efficient DM formation (22). We find that NHEJ may also participate in ecDNA integration (see also the integration junction analysis results). Long-term treatment of SC4 cells with the DNA-PK (also known as PRKDC, a key NHEJ kinase) inhibitor NU7026 at a dose that minimally affects growth (22) decreased the frequency of cells displaying HSR integrations (Fig. 2E). Third, a more pronounced FA mode switch was observed in SC5, with 87% of the cells switching from DM⁺ and HSR⁺ to DM⁻ and HSR⁺, and only 13% retaining the mixed karyotype (Fig. 2B and C). Cells with DM⁺ and HSR⁺ karyotypes appear to reflect an intermediate transition stage in the karyotype switch. In contrast to the initial DM⁺ cases, clone SC2 that only contained HSRs on chromosome 3 at the outset (Fig. 2F) maintained the HSR mode for three months in culture (Fig. 2B and C). However, we did detect HSR plasticity in some cells in terms of duplications and/or translocations of shorter versions of the HSR to other chromosomes, with the concurrent retention of the long HSR (Fig. 2B, SC2-A). Furthermore, long-term inhibition of DNA-PK in the HSR subclone SC2 leads to a lower percentage of multiple HSRs, implicating a role for NHEJ in HSR plasticity (Fig. 2G).

A distinction between the DM and HSR modes of FA was observed during SC long-term culture, in that all subclones with DMs had their BRAF copy number decrease, whereas the BRAF copy number in HSR subclones remained unchanged (Supplementary Fig. S4A and S4B). In addition to de novo karyotype changes, we observed heterogeneity and changes in growth rates over three-month expansion among the SCs (Supplementary Fig. S5A–S5D).

The data above support that de novo integration of BRAF DMs as HSRs did occur under the steady dose of dual MAPKi treatment, likely due to HSRs being a more stable and/or fit mode of FA in the melanoma MAPKi BRAF amplification context. These results are in line with prior findings in different cell types and different treatments (21, 25, 28, 29). However, no changes occurred in one DM⁺ and HSR⁻ clone SC3, indicating that the tendency for integration is not absolute during the time scale observed (Fig. 2B and C). In sum, these SC-based findings demonstrate the plasticity of MAPKi-induced BRAF FAs, with a general trend of fitness-based evolution from DMs to HSRs in these conditions.

The observation that in the M249-VSR system DMs will integrate into chromosomes as HSRs upon continuous culture at a constant drug dose suggests that DM⁺ cells have a fitness disadvantage compared with HSR⁺ cells in these conditions. However, DMs are often observed in tumor samples, and thus may have a fitness advantage in other situations. To test this hypothesis, we aimed to identify a scenario in which DMs would have a fitness advantage. DMs are known to segregate asymmetrically during cell division (13–16), so we tested whether an oscillating drug dose would give DM⁺ cells increased fitness, arguably through increased heterogeneity of the population. We designed an experiment in which we turned the double-drug doses on and off in a cycle of eight days (Fig. 3A, EXP1–2). DMs were indeed retained without a switch to an HSR state for a longer period of time compared with the steady dose scenario (Fig. 3A–C, FIX5; Supplementary Fig. S6A and S6B). However, the number of DMs decreased in these cells, suggesting another MAPK inhibitor resistance mechanism had emerged in these cells. We found that these cells express the shorter BRAF splice isoform associated with acquired resistance, whereas the bulk HSR cells do not show this isoform (ref. 34; Fig. 1D; Supplementary Fig. S6C). Hence, in response to the altered fitness challenge of a regularly changing environment, the emerging cells retained DMs longer than cells experiencing constant drug dose, and the nonconstant conditions furthermore resulted in the expression of an additional resistance-associated BRAF isoform that likely reduced the overall BRAF expression requirement, and thus led to lower DM copy numbers.

Next, we focused on studying the plasticity of DMs and HSRs in M249-VSR cells. As a foundation for this analysis, we first examined whether HSRs were the final stable form of amplicons for cells kept under constant drug dose by checking their karyotypes after a few additional months. We found that most cells still harbored HSRs with similar amplicon length and BRAF copy number (Fig. 3A–D, EXP1). This stable result provides our reference control for comparison with other cases with drug-dose manipulation.

To evaluate if the DM-to-HSR trajectory observed under constant inhibitor dose could be affected by changes in dosing, we next either decreased or elevated the double-drug concentration being applied to DM⁺ or HSR⁺ cells. Previous studies have examined the potential of using drug holidays to eliminate drug-addicted cells (41, 43, 46, 47), thus sparking our interest in studying the effect of this approach on DMs and HSRs. To investigate this, we withdrew VEM + SEL treatment from M249-VSR-DM and M249-VSR-HSR cells. In the DM⁺ case, when doses were acutely brought down from 2 to 0 μmol/L, all DMs were eliminated based on FISH analysis, with the fastest change observed in 12 days. qPCR results showed that the copy number of BRAF was reduced drastically (Fig. 3A–D, EXP3; Supplementary Fig. S7A–S7C). We also performed experiments in which only one of the two drugs was withdrawn. Upon single withdrawal of either drug, we saw substantial, but less-complete, reduction in DM copy number (Fig. 3E and F) in comparison with the double-drug removal. This result supports that in the M249 system, the combined effect of BRAFi and MEKi is a notably stronger amplification selective force than the effect of the single inhibitors. With single-drug withdrawal, there were minimal effects on cell viability and growth rates. In contrast to the DM case, HSR cells upon single-drug withdrawal show reduced viability, supporting HSRs as a less plastic FA mode in this context (Supplementary Fig. S8A and S8B).

A prompt reversion of BRAF copy number to the parental state in three weeks also occurred in HSR cells upon full removal of the dual inhibitors (Fig. 3A–D, EXP4). Notably, there was not a substantial difference between the recovery time of DM and HSR cells in these double-drug washout experiments (Supplementary Fig. S8A and S8B). These results motivated additional experiments to test the plasticity of the HSR FA mode. We next repeated the dose-decrease experiment above using the bulk population in its HSR⁺ state but did not perform a complete withdrawal (Fig. 3A–C, EXP5: 2–0.1 μmol/L). In this experiment, the bulk population demonstrated a substantial shortening of the typical HSR length, but HSRs were still detectable. Using this new subpopulation, we further explored the cellular genomic plasticity by subsequently reinstating the 2 μmol/L double-drug dose. The cells regained resistance in less than a month, and most cells again presented with the longer form of HSRs (Fig. 3A–C, EXP5). During the interval of drug reduction and increase, BRAF DNA copy number also decreased following the 2 to 0.1 μmol/L transition and, accordingly, reincreased following the 0.1 to 2 μmol/L transition (Fig. 3D, EXP5).

We also reinstated a 2 μmol/L drug dose on the bulk population of cells that had drug withdrawal (0 μmol/L) occur while they were in the DM⁺ state (Fig. 3A, EXP3). In this case, it took about four months for the cells to redevelop resistance to VEM + SEL, similar to the time required for the initial establishment of resistance in the parental cells. In this experiment, the melanoma cells demonstrated an additional variation in that upon becoming resistant they typically harbored two or three separate, shorter HSRs on different chromosomes (Fig. 3B and C, EXP3). None of the cells presented with a single larger HSR. This treatment course thus further revealed the plasticity of genomic options available for adjusting to changes in selection pressures.

Notably, we could generalize a subset of these copy-number plasticity findings to other melanoma samples and other amplified genes under MAPKi challenge. The BRAF^V600E cell lines A375 and Mel888 harbor BRAF HSRs upon acquiring dabrafenib (BRAFi) plus trametinib (MEKi) resistance (DTR; refs. 35, 43), and harbor shortened HSRs after dose reduction (Fig. 3G and H). We also characterized NRAS^Q61R or BRAF^S365L melanoma patient-derived xenograft (PDX) models. Upon establishment of MEKi trametinib resistance, these models acquire RAF1 (CRAF) amplifications in extra- or intrachromosomal modes, respectively, and the RAF1 amplifications decreased after drug withdrawal (Fig. 3I–L). To our knowledge, this is the first documented example involving RAF1 DMs mediating MAPKi resistance.

The dose decrease and increase results above could be explained by selection for residual BRAF copy number–low or copy number–high cells in the respective populations. To investigate cellular plasticity to dramatic drug reduction in a more homogeneous population, we turned to the single cell–derived DM⁺ or HSR⁺ clones. In these experiments, we lowered the VEM + SEL double dose from 2 to 0.1 μmol/L using clones SC2 (HSR), SC302 (HSR), SC3 (DM), and SC4 (DM). For controls, we kept the dual dose at 2 μmol/L. In the post–drug decrease populations, SC2 and SC302 showed reduced length of HSRs, and SC3 and SC4 showed a reduced number of DMs. All cases were accompanied by a substantial decrease in BRAF copy number (Fig. 4A–C; Supplementary Fig. S9A–S9C). Another characteristic that indicates the plasticity of HSR-harboring cells is in some regards comparable with that of the DM case is that the recovery times upon dose withdrawal for DM⁺ or HSR⁺ cells, either bulk or as SCs, were not substantially different (Fig. 4D).

Although DM plasticity can be explained by uneven segregation (13–16), HSR plasticity, especially such rapid change in one month, is uncommon during dose challenging—purportedly due to the stability provided by chromosomal integration (48–50). We thus further analyzed the structural data related to the long-to-short HSR transition upon dose reduction. To reduce heterogeneity, we used the HSR⁺ SC2 clone, with its initial long HSR on chromosome 3 (Fig. 2E). In most cells from clone SC2, the post-dose reduction, short HSR remained located on the same chromosome based on FISH staining. Nevertheless, the aforementioned small HSRs (via HSR duplication) on different chromosomes were either not yet present or not favored by selection upon dose reduction, in comparison with the shortening of the Chr3 HSR. Taken together, these results demonstrate the de novo evolution of the clonal long HSR both during constant drug dose (HSR duplication) as well as upon dose reduction (HSR shortening; Fig. 4E–G).

The longer fragment lengths of OM aided in the investigation of the structure of such plastic HSRs. The OM data indicate that the BRAF HSR structure in SC2 is complex and involves duplications (primarily head-to-tail) and some inversions. The HSR was integrated at the PAK2 gene locus near the telomere of chr3 (Fig. 4H and I; Supplementary Figs. S3D and S10A; Supplementary Table S1), in line with previous findings that telomeres and telomere-proximal sites are more frequent locations for integration (22, 28). The PAK2 locus was duplicated, and the integration occurred between the two PAK2 copies (Fig. 4H). Increases in the copy number based on both OM and WGS data support such PAK2 duplication at the site of integration (Supplementary Fig. S10B and S10C). WGS additionally demonstrates that the breakpoint junction between chromosome 3 and 7 contains a two-nucleotide nontemplated insertion, which supports a potential role for NHEJ (ref. 51; Supplementary Fig. S10D) and is line with the aforementioned finding about ecDNA integration dependency on DNA-PK (Fig. 2E).

After VEM + SEL dose reduction, the number of BRAF amplicon repeats decreased, along with the creation of new breakpoints and the generation of a more heterogeneous population. However, the integration junction next to PAK2 was preserved in a subset of the cell population (Fig. 4H and I; Supplementary Table S1). Overall, the combined OM and WGS data support a model of in situ excision of BRAF amplicon repeats (from within the HSR, not from the ends), potentially through error and repair mechanisms, in the long-to-short HSR transitions that are observed upon dose reduction.

We expanded such finding of DM and HSR plasticity to MEKi-resistant subclones from a human NRAS^MUT melanoma cell line (M245), involving different amplified oncogenes. Clone 3 (C3) of M245 cells harbors RAF1 amplification as DM upon becoming resistance to trametinib, whereas clone 5 (C5) harbors NRAS amplification as an HSR. Drug withdrawal caused copy-number decrease in both cases: reducing RAF1 DM number and shortening the NRAS HSR (Fig. 4J and K). The RAF1 and NRAS amplification events have previously been shown to mediate resistance to MEKi in these cell lines (52).

Although these bulk and single-cell clone experiments demonstrate the plasticity of HSRs, we also identified a melanoma cell line with HSR-based focal amplification that does not show shortened or lost HSR upon BRAFi + MEKi removal (Supplementary Fig. S11A–S11D), which is similar to some previous observations and conclusions about HSR stability (48, 49).

To further investigate HSR plasticity, we increased the double MAPKi doses applied to the bulk M249-VSR cells at a timepoint when they were predominantly HSR⁺. Interestingly, this treatment converted the population of predominantly HSR⁺ cells to predominantly DM⁺ cells (Fig. 3A–C, EXP6). Contrary to the expectation that in the higher drug dose the cells would have higher levels of BRAF DNA copy number, we found that the copy number had decreased (Fig. 3D, EXP6).

To investigate this change further, we repeated the experiment using bulk M249-VSR-HSR cells at various time points over the entire HSR-harboring period, roughly 260 days onward from the beginning of resistance development (Fig. 5A). Four out of five dose-increase experiments resulted in changes of FA types from HSRs to DMs (VS5-1, VS5-2, VS5-3, VS5-4, and VS5-6; five sampling points in total). One out of five resulted in cells that were DM⁻ and HSR⁻ VS5-5 (Fig. 5B and C). Notably, we found that the five DM⁺ 5 μmol/L–resistant samples all expressed a BRAF protein variant with a molecular weight of approximately 140 kDa (Fig. 5D). Four of the five DM⁺ 5 μmol/L–resistant samples also expressed the 62 kD variant of BRAF, the BRAFi-resistant splice variant observed in the oscillating dose experiment above. The 140 kD size matches a previously reported BRAF variant with a kinase domain duplication (KDD) that leads to BRAFi resistance (35, 53, 54). Based on RT-qPCR using primer pair that spans the BRAF exon 18–10 junction, we discovered that all of the HSR-to-DM–transformed samples carried exon 18–10 junctions, whereas other cultures, including M249 parental (VS0), M249-HSR cells prior to dose increase (VS2-1, VS2-2, VS2-3, and VS2-4), and M249-HSR cells that showed DM⁻ and HSR⁻ post dose increase (VS5-5), contained none or only a minimal amount of such junctions (Fig. 5E and F). The close to unity ratio of 18–10 to 9–10 junctions supports that each DM unit contains one KDD region in the KDD-expressing sublines.

We next investigated whether the KDD was developed due to selection of an existing subpopulation or de novo kinase domain duplication after the 2 to 5 μmol/L dose increase. Under constant 2 μmol/L dose VEM + SEL, the M249-VSR–resistant cells were initially primarily DM⁺ and HSR⁻ (circa day 150), turned primarily DM⁻ and HSR⁺ with time (circa day 260), and then with additional time reacquired a small percentage of DM⁺ and HSR⁻ cells (450 days and onward; Fig. 5A–C). Their late timepoint DM⁺ and HSR⁻ fractions were 2/46 (4.4%, VS2-3) and 3/30 (10%, VS2-4). This expanding DM⁺ population could have been the source of KDD that expanded post drug-dose increase to 5 μmol/L.

To further test if rare DM⁺ cells were present at earlier times below the level of detection by bulk FISH analysis, we used both single-cell sorting and a replica-plating approach. First, cells from the earlier-stage M249-VSR-HSR bulk population (322 days) were single-cell sorted. This collection of single-cell clones was then treated either at the original 2 μmol/L dose or at an elevated 5 μmol/L dose of VEM + SEL (Supplementary Fig. S12A). We found that 3.2% of the single-cell clones could grow under 5 μmol/L in a similar manner compared with their counterparts at 2 μmol/L (large colonies, Fig. 5G).

Next, we added a replica-plating step. Forty-one single-cell clones derived at 2 μmol/L were replica plated and then treated in parallel at either the original 2 μmol/L dose or the elevated 5 μmol/L dose (Supplementary Fig. S12A). After two rounds of screening, the clone with the highest relative growth rate (SC101) was revealed to be DM⁺ and HSR⁻ both before and after the dose increase, with no observed cellular heterogeneity of FA modes (Fig. 5H and I). The second fastest clone (SC137) started with a 10% DM⁺ and HSR⁻ population, but finished at 100% DM⁺ and HSR⁻ at the end of the replica plating (Fig. 5H and I). Four other randomly selected SCs from either the near top of the relative growth rate–sorted list and the bottom of the list displayed no DM⁺ and HSR⁻ karyotype (SC122, SC124, SC111, SC106; Fig. 5H and I; Supplementary Fig. S12B and S12C). The two fastest SCs, SC101 and SC137, harbored BRAF KDD on their DMs based on immunoblot analysis (Fig. 5J). In a second quantitative viability assay, the SC101 and SC137 KDD⁺ SCs again demonstrated the best ability to tolerate drug-dose increases (Supplementary Fig. S13A). This replica screening result supports that the cells harboring the BRAF KDD–containing DMs preexisted in the bulk population prior to increases in the dual MAPKi dose. These cells were starting to expand in the 2 μmol/L drug condition with a relative fitness slightly higher than other cells, but the increase in drug dose sharply increased such fitness advantage. Using a barcode-based clone tracing system (ClonTracer; ref. 55) to keep track of the subpopulations in the bulk M249-VSR-HSR cells (from day 318), we observed that even under the constant 2 μmol/L drug dose and at such later timepoints, certain cells expanded faster than others (Fig. 5K and L), indicating the bulk population of cells continues to evolve in regard to its subpopulation distributions.

Interestingly, we did not observe de novo generation of DMs, either KDD-bearing or not, from BRAF DM⁻ and HSR⁺ SCs upon performing dual drug (VEM + SEL) escalation (Supplementary Fig. S14A and S14B). Successful HSR-to-DM transitions have been demonstrated in different cell types, involving different genes, with corresponding different drug regimens. These reported cases typically involve the creation of fragile sites or chromothripsis on HSRs (22, 26).

After learning the plasticities of BRAF-containing DMs and HSRs in response to dose perturbations, we next investigated whether amplicon boundaries and junctions changed during these processes. We performed structural variant (SV) analysis on the M249 samples collected after various dose challenges (Fig. 6A) and found that genomic amplification boundaries do not differ substantially regardless of their FA mode, amplicon number, subcloning status, nor presence or absence of the KDD selection, supporting a single initial amplicon origin (Fig. 6B). This conclusion is further corroborated by the conserved junctions connecting amplicons in bulk and SC DM⁺ and HSR⁺ sublines cultured at full drug dose (Supplementary Table S2). New junctions were generated during dose decreases and KDD formation, but these alterations did not alter the overall amplification coordinates (Fig. 6B).

High BRAF amplification upon MAPKi-acquired resistance is observed in the clinic with at least 16 patient-based cases reported in the literature, in PDX models (35) and in other cell line models (Supplementary Tables S3 and S4). In amplicon boundary analysis of the cohort, we found the M249 FA boundary to be consistent with the range of all observed boundaries. Furthermore, we did not find evidence for inclusion of coamplification loci adjacent to BRAF (Fig. 6C). This observation is in line with some frequently amplified oncogenes, around which the ecDNA breakpoints distribute randomly (30), and is distinct from MYCN and EGFR amplification in other tumor types that have been shown to involve coamplification of adjacent enhancers (11, 12). Looking across the cohort, we did not observe a relationship between pretreatment BRAF copy number and the likelihood for BRAF focal amplification post acquisition of resistance (Supplementary Fig. S15A).

Given the high plasticity of BRAF amplifications in the M249-VSR series, we next investigated cellular vulnerabilities affiliated with amplification in this dual MAPKi context. Our previous study revealed a correlation between melanoma differentiation stages and sensitivity to proferro-ptotic drugs (39). We thus tested if BRAF-amplified cells have altered sensitivity to disruption of the repair of oxidized lipids. We tested the sensitivity of the proferroptotic drug RSL3, which targets glutathione peroxidase 4 (GPX4), in the M249-VSR series. We found that both M249-VSR-DM and M249-VSR-HSR are substantially more sensitive to RSL3 compared with M249-P (Fig. 7A). We further found that after drug withdrawal, when the cells have reduced BRAF copy number, M249 cells lose sensitivity, reverting to levels closer to the parental cells. Consistently, two additional cases of BRAFi + MEKi resistance mediated by BRAF amplification, A375-DTR and Mel888-DTR (43), also showed increased RSL3 sensitivity compared with their parentals. We next confirmed that the RSL3 sensitivity in M249-VSR sublines has the expected characteristics of ferroptosis, namely that the RSL3 sensitivity is reactive oxidative species (ROS)–, lipid ROS–, and iron-dependent, as cell death can be rescued by adding reduced glutathione (GSH), the lipophilic antioxidant Trolox, and the iron chelator deferoxamine (DFO; Fig. 7B and C). We confirmed an increase in lipid ROS levels in M249-P and M249-VSR cells upon RSL3 treatment, which was protected by the presence of a lipophilic antioxidant (Supplementary Fig. S16A).

The BRAF-amplified M249 dual MAPKi-resistant cells that are sensitive to the GPX4 inhibitor RSL3 were also sensitive to the proferroptotic drug ferroptocide that targets thioredoxin (56), but were not more sensitive to inhibition of the system x_c⁻ cystine/glutamate antiporter by erastin (Supplementary Fig. S16B). Such differential sensitivity to different upstream components of the glutathione synthesis and ferroptosis pathway has been previously observed (e.g., SKMEL28R; ref. 39).

We next investigated why these three cell lines with BRAF amplification–mediated and thus MAPK reactivation–mediated resistance demonstrated higher ferroptosis sensitivity compared with their parental sublines. In previous studies, ferroptosis sensitivity in melanoma is associated with innate or acquired treatment-induced dedifferentiation (39). However, in our past work, resistance mediated by reactivation of the MAPK pathway through genomic changes (e.g., via NRAS mutation: M249P/R) does not lead to dedifferentiation and changes in ferroptosis sensitivity (39). We thus analyzed whether the three BRAF-amplified dual MAPKi-resistant cells studied here demonstrated signs of dedifferentiation. Gene-expression profiles of RSL3-sensitive M249-VSR (DM and HSR amplification mode), A375-DTR and Mel888-DTR cells (both with HSR mode) do not demonstrate dedifferentiation compared with their parental sublines when their gene-expression profiles are projected onto a panel of melanoma lines spanning the full spectrum of differentiation states (39). By contrast, the cell lines M229P/R, M238P/R, and SKMEL28P/R, which became resistant through upregulation of receptor tyrosine kinases (RTK; ref. 33), demonstrated dedifferentiation and increased sensitivity to RSL3 (as tested previously; ref. 39; Supplementary Fig. S16C and S16D). Such findings are also supported by the combination of increases in melanocyte differentiation and pigmentation, decreases in mesenchymal gene set scores, and increases in the melanoma differentiation master regulator MITF in the BRAF amplification samples, and the reverse patterns in the RTK upregulation/dedifferentiation cases (Supplementary Fig. S16E–S16H).

Upon determination that increased RSL3 sensitivity in these three BRAF amplification cases is not due to dedifferentiation, we then turned our focus to mitochondrial pathways, as previous studies have reported that MAPKi resistance can cause melanoma cells to shift their major energy generation program from glycolysis to mitochondrial pathways. This shift then leads to elevated production of ROS and more dependence on ROS detoxifying mechanisms (57–61). Although ROS were implicated, these prior studies did not assess the change in ferroptosis sensitivity of the resistant sublines. We found that upon acquisition of BRAFi/MEKi resistance, M249, Mel888, and A375 all upregulate PPARGC1A (PGC1α; ref. 57), have distinct but overlapping patterns of upregulation of mitochondrial respiration programs [tricarboxylic acid cycle (TCA), electron transport chain, oxidative phosphorylation, and mitochondrial biogenesis], and all upregulate lipid oxidation pathways (featured by PPARα and ACOX1; ref. 62). In sum, these changes may cause higher dependence on glutathione metabolism for lipid detoxification via GPX4, whereas all cases do not equally upregulate their ROS detoxification pathways (Supplementary Fig. S16E and S16F).

In accordance with this last observation, (i) expression of the ROS detoxification pathway gene glutathione synthetase (GSS) and inferred activity of the ROS detoxification pathway are downregulated in both dedifferentiation and BRAF amplification cases of MAPKi resistance (Supplementary Fig. S16G and S16H), and (ii) the levels of reduced GSH are decreased upon both dedifferentiation (39) and BRAF amplification (Supplementary Fig. S16I). We furthermore found that the NCOA4 gene, which mediates the selective autophagic degradation of ferritin (63), is upregulated in both dedifferentiation- and BRAF amplification–mediated MAPKi-resistant cells, but is downregulated in NRAS mutation–mediated resistance (where both parental and resistant sublines have similar ferroptosis sensitivity; ref. 39; Supplementary Fig. S16E). This observation is in line with a previous study finding that NCOA4 promotes accumulation of cellular labile iron, leading to higher susceptibility to proferroptotic drugs (64).

Taken together, although dedifferentiation-mediated MAPKi resistance has distinctions from BRAF amplification–mediated resistance, they both demonstrate increased GPX4 inhibition (RSL3) sensitivity, have a common pattern of downregulated ROS detoxification genes such as GSS, and demonstrate upregulation of the iron homeostasis regulator NCOA4.

Focal amplifications of oncogenes in either DM (ecDNA) or HSR mode are clinically observed both as a resistance mechanism for inhibitors targeting oncogenes (e.g., MET in EGFRi-treated lung cancer; ref. 65) and in the targeted therapy–naïve setting (e.g., MYCN in neuroblastoma; refs. 27, 66). The disappearance of oncogene-containing DMs has also been reported upon modeling of oncogene-targeted therapy (67). Although a few-fold amplification of BRAF is sometimes observed in treatment-naïve melanoma tumors (68), higher-fold DM or HSR focal amplifications are typically seen only following MAPK inhibitor therapy (Supplementary Tables S3 and S4). To further elucidate the genomic plasticity enabled by focal amplifications, we developed an expanded version of a BRAF + MEK inhibition and BRAF locus amplification model. This system demonstrated a high degree and broad range of evolutionary plasticity of BRAF amplicon in response to changing drug-dose regimens, which can in part be generalized to other amplified MAPK genes mediating resistance (i.e., RAF1 and NRAS). BRAF plasticity was in cases coupled to multiple genomic rearrangement and related mechanisms such as kinase domain duplications and alternative splicing.

In the initial phase of drug resistance to dual BRAF and MEK inhibition, the BRAF amplification appeared via DMs. Under conditions of a stable double-drug dose, the population gradually became dominated by an HSR form of BRAF amplification. Such DM-to-HSR conversion was also observed in single cell–derived clones of the M249-VSR cells, supporting that de novo integrations of (potentially agglomerated; ref. 45) DMs occurred in addition to selection of an existing HSR⁺ population. Such FA mode switch is also supported by the conserved genomic contents and shared breakpoints between DM and HSR amplicons based on WGS and OM data. This mode switch result adds to reports in the literature for other focally amplified oncogenes in different cancer types. For example, one study conducted a long-term observation on a non–drug-treated leukemia cell line and saw formation of MYC-carrying HSRs from DMs (25). Another study proposed a common origin of MYCN DM and HSR in neuroblastoma based on their shared structures (69).

The reproducible observation of the DM-to-HSR transition led us to hypothesize that DMs carry a higher fitness disadvantage than HSRs during stable conditions. In support of this, we found that an oscillating drug dose could prevent or prolong autosomal integration of the amplicon. This difference in fitness is arguably linked to uneven segregation of DMs (13–16) and the resulting uneven BRAF gene copy numbers in daughter cells. During nonstable conditions, cellular heterogeneity provided by uneven segregation can provide a reservoir of cells more adept to grow well in the new conditions. In contrast, during stable drug-dose conditions, a reduction in cellular heterogeneity would produce a fitness advantage, with all daughter cells maintaining the optimal BRAF gene copy number. We also observed that DM numbers tended to decrease during long-term stable culture, probably suggesting that a secondary (undetermined) resistance mechanism allowed these cells to depend less on the DMs (Supplementary Fig. S4A and S4B). Taken together, our results, along with other findings on the evolution of DMs harboring different oncogenes in different cancer types (25, 28, 29), support that in some cell types in nonchanging contexts, DMs are not a fitness-optimized form of amplification and thus tend to be replaced by other mechanisms such as less heterogeneous chromosomally integrated HSRs. However, such fitness considerations are likely affected by cell type and by the characteristics of the oncogene driving the focal amplification.

In contrast, in nonconstant conditions, such as the tumor microenvironment or tumors targeted by therapeutics, the uneven segregation of DMs provides an evidence-supported model for tumor heterogeneity that in turn provides tumors the diversity to withstand changes in conditions that affect fitness (13–16, 70). In the single- and double-drug withdrawal experiments involving DM⁺ cells, we saw rapid decreases in the DM copy number (e.g., BRAF or RAF1). It is possible that the rapid changes in DM copy number were due to selection of a preexisting DM-negative subpopulation or that post-mitosis cells with fewer DMs due to uneven segregation could have been selected for upon drug withdrawal (13–16). It is also possible that DMs were exported out of cells through previously observed micronuclei exclusions (71), especially in the single-drug withdrawal cases where decreases in DM copy number occurred without appreciable changes in cell viability or growth rates. The single-drug withdrawal results also support that dual BRAF and MEK inhibition is required to sustain pressure for high copies of the BRAF gene.

Beyond DM plasticity, our study revealed that “HSR plasticity” can also be a mode of tumor evolution in response to drug challenge. Dose reduction experiments demonstrate that HSRs can offer somewhat comparable levels of plasticity as DMs. Due to the inherent differences between DM and HSR modes of amplification, this is almost undoubtedly through distinct molecular mechanisms. In more detail, we observed single cell–derived HSR-containing cell populations that demonstrated dose-tunable BRAF and RAF1 HSR lengths. OM, WGS, and FISH data reveal that such length shorting involves reducing the number of amplicon repeats rather than changing integration junctions (Fig. 4H and I). Future work will investigate whether errors and repairs made while replicating and segregating intrachromosomal long HSRs may be generating heterogeneity and thus contributing to this plasticity.

In sum, the single-cell clone results support that de novo genetic alterations occur during expansion from a single cell and/or during the stress of drug withdrawal, thus creating population heterogeneity and enabling population plasticity. In these cases, selection alone cannot explain the outcome, and clearly genomic instability, in the HSR case potentially mediated by the challenge of replicating adjacent homogeneous regions, is diversifying the population.

The tumor evolutionary and drug resistance plasticity enabled by focal amplifications extended beyond changes in amplicon copy numbers and DM versus HSR modes. In particular, we observed two additional parallel mechanisms: (i) KDD, representing an additional genomic rearrangement mechanism (35, 53, 54), and (ii) activation of an alternative splicing mechanism (34). Our results indicate that cells harboring the BRAF KDD–encoded DM mechanism can be reproducibly selected from an HSR-predominant population upon dual MAPKi escalation treatment due to an accompanying gain in relative fitness advantage. Further research on KDD formation and KDD-mediated resistance could offer therapeutic insights for pan-cancer therapy, as this alteration occurs to many other kinases, such as EGFR and FGFR1 in glioma and lung cancer (72–75). We also observed the alternative splicing mechanism as a potential method to escape reliance on high DM copy number during an oscillating dose regimen. The drug resistance provided by the splice variant arguably lowers the number of DMs required but maintains the DM-mediated unequal segregation-based heterogeneity.

Therapeutic approaches to target the vulnerabilities of FA-harboring cells are in academic and industry development. Our study demonstrates important challenges, such as mode switching and acquisition of additional genomic rearrangements, that must be coaddressed in these pursuits. Here we report that BRAF-amplified melanomas relapsed from dual MAPKi treatment show increased ferroptosis sensitivity, which extends the spectrum of ferroptosis sensitivity in melanoma therapy resistance. We found that an additional mechanism, distinct from treatment-induced dedifferentiation and mesenchymal transition, can generate sensitivity to GPX4 inhibition (39). This finding links to studies of MAPKi-induced oxidative stress in melanoma. In some melanomas, BRAF^V600E activation leads to enhanced glycolysis and reduced oxidative phosphorylation and mitochondrial respiration (57). However, BRAF inhibition, including acquired resistance to BRAF inhibitors, can switch the energy generation dependency back to oxidative phosphorylation pathway by induction of PPARGC1A and overexpression of other mitochondrial genes (57–60, 76, 77). Reactive oxygen species (ROS) productively mediate redox-based energy production in mitochondrial respiration, but they can also damage lipid, protein, and DNA (78). Hence, respiring cells need to upregulate detoxification programs to compensate for elevated oxidative stress (61, 79). The imbalance of cellular pro-oxidative and antioxidative mechanism can lead to cell death (80). Ferroptosis is one form of cell death that can result from such compromised redox homeostasis, mediated by iron-dependent accumulation of lipid peroxides (81).

In our studies, BRAF amplification–mediated MAPKi-resistant melanoma cells did not exhibit dedifferentiation. However, they did downregulate GSS and had limited reduced glutathione levels, which would limit their capacity to detoxify lipid ROS (Supplementary Fig. S16E–S16I). They also upregulated the iron homeostasis regulator NCOA4 (refs. 64, 82; Supplementary Fig. S16E), similar to other MAPKi parental/resistant melanoma pairs with differential RSL3 sensitivity, consistent with higher vulnerability to ferroptosis induction. Relatedly, one previous report found that MAPKi-acquired resistance through most MAPK reactivation mechanisms, such as RTK overexpression, NRAS^Q61H/Q61K mutation, KRAS^G12C mutation, BRAF splice variant, and BRAF amplification, are all more vulnerable to undergo apoptosis via inhibition of the system x_c⁻ cystine/glutamate antiporter in cells treated with a histone deacetylase inhibitor (83). Our finding complements this by uncovering a different form of cell death, ferroptosis, occurring under a similar MAPKi-resistance context and extends the role of ferroptosis in MAPKi resistance beyond cases of dedifferentiation. Taken together, the melanoma dedifferentiation-independent synthetic lethality between BRAF amplification and ferroptosis identified here provides therapeutic insight for treating BRAF-amplified melanomas relapsed from MAPKi treatment. Although ferroptosis sensitivity was observed in both the HSR- and DM-harboring cells, previous reports have revealed that DM and HSR can have specific targetable vulnerabilities linked to their distinct mechanisms for generation and maintenance (84–90). Future work is needed to confirm such FA mode–specific vulnerabilities in BRAF amplification systems.

Collectively, we observed a high degree and broad range of tumor evolution and drug-resistance plasticity enabled by or coupled to focal amplifications. Through perturbations by a panel of drug regimen challenges, we observed (i) de novo generation of extrachromosomal DMs, (ii) de novo integration of DMs into chromosomal HSRs, (iii) context-dependent HSR-mediated fitness advantage over DMs, (iv) context-dependent DM-mediated fitness advantage over HSRs, (v) coevolution of DMs and a de novo genomic rearrangement creating a kinase domain duplication, (vi) coevolution of DMs and activation of BRAF alternative splicing, (vii) propensity to couple secondary resistance mechanisms (KDD and/or alternative splicing) to DMs to reduce the total number of DMs required, and (viii) a plasticity of HSRs that compares in some kinetic aspects to the known plasticity of DMs. Appreciation of the interplay of focal amplification modes with drug regimens and other resistance mechanisms is central to our understanding of tumor evolution and drug resistance, and to developing therapeutic approaches to overcome the resulting plasticity.

The M249 (RRID: CVCL_D755), M395 (RRID: CVCL_XJ99), and M245 NRAS^Q61K (RRID: CVCL_D754) cell lines are part of the M series melanoma lines established from patient biopsies at UCLA under UCLA Institutional Review Board approval #02-08-06 and were obtained from Dr. Antoni Ribas (91). PDX1 (NRAS^Q61R) and PDX13 (BRAF^S365L) cell lines were derived from PDXs with the same names (47). M245 C3 and C5 sublines were reported previously (52). Mel888 (RRID: CVCL_4632) and A375 (RRID: CVCL_0132) cell lines and their variants were described previously (35, 43). All cell lines have been tested for Mycoplasma. All cells were cultured in RPMI 1640 with L-glutamine (Gibco), 10% (v/v) FBS (Omega Scientific), and 1% (v/v) streptomycin (Gibco). All cells were maintained in a humidified 5% CO₂ incubator. Resistant M249 cell lines were generated by exposing cells to step-wise increasing doses of vemurafenib and selumetinib, similar to the previously described approach (41). Briefly, the doses for both drugs were sequentially increased by roughly 2-fold, with each dose escalation taking place when cells resumed growth rates with doubling in four days or less. The initial and final doses were 0.05 and 2 μmol/L, respectively. Growth and viability were assayed by staining cells with trypan blue (Sigma-Aldrich) followed by cell counting using Vi-cell XR Cell Viability Analyzer (Beckman Coulter) or by CellTiter-Glo luminescence assay. Doubling times for M249 SCs and bulk cells were calculated by fitting exponential growth curves, and their error bars were derived based on a previously published method (92). Cells were sampled only for experiments when they show reasonable growth rate at corresponding dose.

BRAF inhibitors vemurafenib and dabrafenib as well as MEK inhibitors selumetinib and trametinib were obtained from Selleckchem or LC Laboratories. Proferroptotic drugs RSL3 and erastin were obtained from Cayman Chemical and Selleckchem, respectively. Ferroptocide was described previously (56). The DNA-PK inhibitor NU7026 was purchased from Selleckchem. Inhibitors were all dissolved in DMSO.

Resistant subclones were derived by seeding single cells from the bulk population into 96-well plates using a FACSAria cell sorter. Doublets were removed by circling the right area in the FSC height versus area plot. Seeded single cells were then cultured using aforementioned medium or a modified medium with 20% FBS for two weeks. Culture medium was not changed until clear colonies were observed in some wells. If certain treatments were needed, i.e., double-drug dose changes, they were initiated upon seeding the cells. M245-resistant subclones were derived by ring selection (47).

Cells were blocked at metaphases by adding colcemid (KaryoMax, Thermo Fisher Scientific) at a final concentration of 0.05 μg/mL followed by incubation at 37°C for six to eight hours. Cells were then fixed using methanol:acetic acid (3:1). FISH slides were prepared by dropping fixed cells in a humid environment following the manufacture's protocol provided by Cytotest and Empire Genomics. Formalin-fixed, paraffin-embedded (FFPE) xenograft tumor FISH slides were prepared by pepsin digestion followed by similar procedures of cell line FISH. Colored FISH images were taken and processed using confocal microscope Leica TCS SP8 X. Karyotype categorizations were based on the guidelines in Supplementary Fig. S1. The fractions under certain images represent the number of cases for corresponding karyotype divided by the total number of cases analyzed. If not otherwise mentioned, scale bars in FISH images represent 10 μm. Centromere probe names are abbreviated as CEN-x. DM numbers were quantified by directly counting the number of features in the FISH images or by using ecDNA quantification tool EcSeg (93). HSR lengths were quantified by dividing the probe area by chromosomal DAPI area in metaphases. The staining areas were calculated using ImageJ v1.53a. Cells fixed by the same procedure were also used for G-banding. G-banded metaphase spreads were photographed using the 80i Nikon Microscope and Applied Spectral Imaging Karyotyping system. A minimum of 10 metaphases were karyotyped.

qPCRs for BRAF genomic DNA (gDNA) copy-number measurement were performed by combining samples with PowerUp SYBR Green Master Mix (Applied Biosystems) in Optical 96-Well Reaction Plates (Applied Biosystems) with three technical replicates for each sample. Plates were then read by the 7500 Real-Time PCR System (Applied Biosystems) using the standard cycling mode. Input templates for all samples were gDNAs extracted using DNeasy Blood and Tissue Kits (Qiagen). Unless specified, all qPCR runs used M249 parental as the reference sample and GAPDH as the endogenous control. Error bars represent t-distribution–based 95% confidence intervals from triplicates: |${\rm{R}}{{\rm{Q}}_{{\rm{max}}/{\rm{min}}}}{\rm{\ }} = {\rm{\ }}{2^{ - \Delta \Delta {\rm{Ct}} \pm {{\rm{t}}_{0.05,{\rm{df}}}}{\rm{*SE}}}}$|⁠. RQ, relative quantity; Ct, threshold cycle; df, degree of freedom; SE, sample standard error. All primers were ordered from Eurofins Scientific and their sequences are shown below.

• BRAF Forward: 5′-TTTAGAACCTCACGCACCCC-3′ (intron 2)

• BRAF Reverse: 5′-TGTTGTAGTTGTGAGCCGCA-3′ (intron 2)

• GAPDH Forward: 5′-CTGGCATTGCCCTCAACG-3′

• GAPDH Reverse: 5′-AGAAGATGAAAAGAGTTGTCAGGGC-3′

gDNA of M249-P and M249-VSR cells was isolated by using DNeasy Blood and Tissue Kits (Qiagen). Samples were run on Agilent 6 × 80K array. The raw data were then processed by Cytogenomics software v5.2 (Agilent Technologies). Nested genomic regions were flattened, and .seg files were generated, followed by data visualization in Integrative Genomics Viewer (IGV) v2.10.0 (94). Regions with large copy-number changes were identified by comparing every segment in M249-VSR with the corresponding segment in M249-P. The same gDNAs were sent to PacGenomics for low-pass WGS with coverage of 0.04. The library was prepared using the KPA DNA Library Preparation Kit. Sequencing was performed on Illumina NextSeq 500 using 75-bp paired-end reads (2 × 75 bp). CNA was inferred using Ginkgo v3.0.0 (95), which contains a step that used bowtie v1.2.1 (96) to align raw reads to hg19 genome.

gDNA of M249-P and M249-VSR sublines was extracted by DNeasy Blood and Tissue Kits. The samples underwent WGS library preparation and were then sequenced on Illumina Novaseq S1 at 2 × 150 and 10–15× coverage. Raw reads in fastq files were aligned to hg38 using BWA-MEM v0.7.1 (97). The duplicated reads were marked by MarkDuplicates tool from GATK v4.1.2 (98). Next, CNA calls were performed using CNVkit v0.9.7 (99) with flat normal as the control. Segmentation was performed using the HMM (hidden Markov model) tumor method. CNVkit results were used as the input for AmpliconArchitect v1.2 (100). The same genomic region chr7:139410000–141180000, which corresponds to the BRAF amplicon, was used as the seed interval for all M249 samples when running AmpliconArchitect. SVs were also called using SvABA v1.1.3 (101) for analyzing break points and integration junctions.

AmpliconArchitect-generated breakpoint graphs were first converted to in silico–digested optical map segments. AmpliconReconstructor v1.01 (ref. 45; https://github.com/jluebeck/AmpliconReconstructor) was then run with default settings on the breakpoint graph segments and the assembled Bionano contigs from the Bionano Genomics optical genome map de novo assembly pipeline. From the collection of reconstructed breakpoint graph paths present, we identified circular or noncircular paths representing the ecDNA or HSR structures. Resulting structures were visualized with CycleViz v0.1.1 (https://github.com/jluebeck/CycleViz).

Ultra-high-molecular-weight (UHMW) DNA was extracted from frozen cells preserved in DMSO following the manufacturer's protocols (Bionano Genomics). Cells were digested with Proteinase K and RNAse A. DNA was precipitated with isopropanol and bound with nanobind magnetic disks. Bound UHMW DNA was resuspended in the elution buffer and quantified with Qubit dsDNA assay kits (Thermo Fisher Scientific). DNA labeling was performed following the manufacturer's protocols (Bionano Genomics). Standard Direct Labeling Enzyme 1 (DLE-1) reactions were carried out using 750 ng of purified UHMW DNA. The fluorescent-labeled DNA molecules were imaged sequentially across nanochannels on a Saphyr instrument.

De novo assemblies of the samples were performed with Bionano's de novo assembly pipeline (Bionano Solve v3.6) using standard haplotype-aware arguments. With the overlap–layout–consensus paradigm, pairwise comparison of filtered DNA molecules (average length approximately 350 kbp) of >200× coverage was used to create a layout overlap graph, which was then used to generate the initial consensus genome maps. By realigning molecules to the genome maps (P value cutoff < 10⁻¹²) and by using only the best-matched molecules, a refinement step was done to refine the label positions on the genome maps and to remove chimeric joins. Next, during an extension step, the software aligned molecules to genome maps (P < 10⁻¹²), and extended the maps based on the molecules aligning past the map ends. Overlapping genome maps were then merged (P < 10⁻¹⁶). These extension and merge steps were repeated five times before a final refinement (P < 10⁻¹²) was applied to “finish” all genome maps.

Optical map alignment of Bionano contigs to the reference genome and Bionano raw molecules to the reference genome was performed with FaNDOM v0.2 (ref. 102; https://github.com/jluebeck/FaNDOM). Alignment and SV detection were done by calling modules named “wrapper_contigs.py” and “wrapper_individual.py” with default settings. Alignments were visualized with the MapOptics v1.0.0 (103) software.

Total RNA was isolated from M249-P and M249-VSR cells using the RNeasy Plus Mini Kit (Qiagen). Samples were sequenced on HiSeq3000 at 150-bp paired-end or on Novaseq SP at 50-bp paired-end. Raw data were then processed using the Toil v3.9.1 pipeline to output transcripts per million, including STAR v2.7.1a, that aligned raw reads to the GRCh38 genome (104, 105). Following data trimming and log transformation, visualizations were done in R. For calculating allele frequencies of BRAF^V600E, all RNA-seq fastq files were aligned using STAR, and the resultant bam files were processed according to GATK RNA-seq short variant discovery best practices until the step of haplotype calling (106). For visualization, we loaded base quality score recalibrated bam files to IGV.

Cell lysates were prepared by using mRIPA buffer supplemented with PMSF, leupeptin, and aprotinin. Western blots were performed using following antibodies: β-actin (AC-15, Sigma-Aldrich), β-actin (13E5, Cell Signaling Technology), BRAF (F-7, Santa Cruz Biotechnology), BRAF (C-19, Santa Cruz Biotechnology), goat anti-rabbit secondary antibodies (IRDye 680RD, LI-COR), and goat anti-mouse secondary antibodies (IRDye 800CW, LI-COR). Images were directly output by the Odyssey CLx Imaging System (LI-COR).

Total RNA was extracted from fresh cells using the RNeasy Plus Mini Kit (Qiagen). Reverse transcriptions were then performed by using the SuperScript VILO cDNA Synthesis Kit (Invitrogen). cDNA was then used for PCR and qPCR. Primers for detecting BRAF exon 18–10 and exon 9–10 junctions were the same as those previously published (35). The regular PCR was performed using Phusion High-Fidelity PCR Master Mix with HF buffer (New England Biolabs). The PCR products that targeted exon 18–10 and exon 9–10 were then combined for each sample and run on 2% agarose gel. For qPCR, each sample was combined separately with PowerUp SYBR Green Master Mix (Applied Biosystems) and loaded on Optical 96-Well Reaction Plates (Applied Biosystems) in triplicate. Plates were then read by the 7500 Real-Time PCR System (Applied Biosystems) using the standard cycling mode.

Each of 41 SCs of M249-VSR-HSR cells (cultured at 2 μmol/L VEM + SEL) was seeded in 6 wells of 96-well plates with the same cell number per well. Three wells of each clone were treated with 5 μmol/L VEM + SEL, whereas the other three stayed at 2 μmol/L. After six days, cell viabilities were measured by CellTiter-Glo Luminescent Cell Viability Assay. Thirteen of 41 SCs were picked for a second round of the dose-increase screen to confirm the findings. The viability of SCs was visualized by heat maps using the R package ComplexHeatmap v2.6.2 (107).

ClonTracer Barcoding Library (55) was purchased from Addgene. The plasmid pool was expanded by electroporation transformation. Lentivirus was made by transfecting 293T cells. M249-VSR-HSR cells were tested for their puromycin dose–response and multiplicity of infection curves. For the actual infection, 54 million M249 HSR cells were spin-infected in 12-well plate with 8 μg/mL polybrene, followed by a six-day puromycin (0.3 μg/mL) selection. Day 0 refers to the end of the selection. Next, cells underwent a standard culture growth period with kinase inhibitors present until the gDNA collection points on days 14 and 35. The sequencing library was prepared by PCR amplification of barcoded regions using the primer sequence provided by the manufacturer. The libraries were paired-end sequenced on Illumina NextSeq 500 at 75 bp read length.

MAPKi-treated melanoma copy-number profiles from multiple previous studies were downloaded and compiled. The list of studies (5, 22, 33, 35, 41, 108–114) can be found in Supplementary Table S3. The software packages used for CNA callings include CNVkit v0.9.7 (99), penncnv v1.0.5 (115), CopywriteR v1.3 (116), rCGH v1.20.0 (117), and the circular binary segmentation algorithm (118). Supplementary Table S4 contains a subset of samples in Supplementary Table S3 that have paired pretreatment and postprogression time points. When the actual normal samples are not available for certain patients, flat normals were used to call CNA. Gene-level copy numbers of BRAF were determined by averaging all length-normalized segments in the BRAF genomic region after removing the gaps.

For Fig. 6C, the expectation of amplicon boundaries around BRAF was simulated using a method from a previous study (11). Briefly, to construct the solid line, real copy-number (CN) profiles of treated melanoma samples used include those that have both pretreatment and postprogression time points with BRAF CN log₂(post/pre) >0.75 and BRAF CN log₂(post/normal) >1.3 as well as those do not have pretreatment data available and BRAF CN log₂(post/normal) >1.7. We have confirmed all selected samples have focal BRAF amplification instead of the arm level. For the dashed line, random amplicons were generated by shifting each of real BRAF amplicon boundaries multiple times but still encompass the BRAF gene. The boundaries are sometimes defined after merging nearby CNA segments with log₂ differences within 1 and gaps smaller than 1 Mb. The genome was binned at 10 kb size. For each bin, amplification frequency is defined by the percentage of samples that have BRAF CN log₂(post/normal) >1.3.

The dose–response curves of ferroptosis-inducing agents RSL3, erastin, and ferroptocide were performed by seeding an appropriate number of cells on day 0 in 96-well plates, treating cells on day 1 with corresponding drugs and reading the plates on day 4 for viability using the CellTiter-Glo luminescent assay. If not otherwise mentioned, resistant cells were maintained in the full dose of MAPK inhibitors throughout the dose–response experiments to keep the BRAF amplifications. Seeding density for each cell line was determined by using the same assay and the same experimental length with multiple cell number titrations. The dose series were generated by serial dilutions. All drugs used for dose–response curves were dissolved in DMSO. DMSO toxicity was performed on the cell lines to determine the appropriate DMSO concentration (0.5%), which was used in all doses. The resulting values from viability assays were normalized to the zero-dose condition after subtracting background (wells with no cells). The curve fittings were performed by using a three-parameter model in the drc v3.0-1 R package (119).

Cells (6,000 per well) were seeded in 96-well plates and treated with ferroptosis-inducing agents in combination with vehicle, GSH (Sigma, G4376), Trolox (Acros Organics, 218940010), or DFO (Sigma, D9533) the next day. CellTiter-Glo luminescence was assessed 24 hours after treatment. For quantification, all values were normalized to vehicle conditions. Resistant cells were maintained in the full dose of MAPK inhibitors throughout the treatments to keep the BRAF amplifications.

In 12-well plates, 80,000 M249-P and M249-VSR-DM cells were seeded per well and treated the next day with RSL3 in combination with Trolox or vehicle. Resistant cells were maintained in the full dose of MAPK inhibitors to keep the BRAF amplifications. After 24 hours, CM-H2DCFDA dye (Invitrogen C6827) was added to each well and incubated for another 20 minutes at 37°C. Cells were then washed with PBS, harvested by trypsinization, suspended in 250 mL PBS, and filtered through cell strainers. The samples were analyzed with BD LSRII Analytic Flow Cytometer at the excitation wavelength of 488 nm.

An appropriate number of cells were seeded in 10-cm dishes for 72-hour growth to reach 80% confluency. MAPKi-resistant cells were maintained in the full dose of VEM + SEL to keep the BRAF amplifications. On the day of collection, cells were rinsed with ice-cold 150 mmol/L NH₄AcO at pH 7.3, incubated with 80% MeOH at −80°C for 20 minutes, scraped off from the plates and transferred to Eppendorf tubes. Cells were then vortexed for 10 seconds and centrifuged at 16,000 × g for 15 minutes at 4°C. The supernatants were transferred to glass vials and dried in a Genevac EZ-2 Elite evaporator at 30°C to obtain metabolite extracts.

Dried metabolites were resuspended in 50% acetonitrile (ACN):water, and one tenth was loaded onto a Luna 3 μm NH2 100 Å (150 × 2.0 mm) column (Phenomenex). The chromatographic separation was performed on a Vanquish Flex (Thermo Scientific) with mobile phases A (5 mmol/L NH4AcO pH 9.9) and B (ACN) and a flow rate of 200 μL/minute. A linear gradient from 15% A to 95% A over 18 minutes was followed by a 9-minute isocratic flow at 95% A and re-equilibration to 15% A. Metabolites were detected with a Thermo Scientific Q Exactive mass spectrometer run with polarity switching (+3.5/−3.5 kV) in full scan mode with an m/z range of 70–975 and 70.000 resolution. TraceFinder 4.1 (Thermo Scientific) was used to quantify the targeted metabolites by area under the curve using expected retention time and accurate mass measurements (<5 ppm). Values were normalized to protein content of extracted material. Data analysis was performed using in-house R scripts.

Raw RNA-seq data of Mel888, Mel888-DTR, A375, A375-DTR, SKMEL28P, SKMEL28R, and M series cell lines were downloaded from corresponding Gene Expression Omnibus (GEO) accessions (39, 43, 120–124). The data were processed through Toil v3.9.1 (105) to obtain RSEM (125) expected counts and normalized by log-transformed counts per million (logCPM) approach. Principal component analysis was performed using mean-centered logCPM values of M series cell lines and served as the framework on which RNA-seq data of other samples were projected onto for determining their dedifferentiation stages. The score of selected gene sets for the parental/resistance-paired cell lines was calculated using the single-sample gene set enrichment analysis (ssGSEA) method in GSVA v1.38.2 R package (126). Nearly all gene sets were taken from MSigDB v7.4 (127, 128) except for ROS detoxifying gene sets which were made by combining (i) a subset of detoxifying genes (a combination of multiple detoxifying gene sets in MSigDB) that correlate well (Pearson correlation >0.4) with the dedifferentiation trajectory scores of M series samples and (ii) the top eight genes that downregulate upon knocking down PGC1α in A375 cells (61).

Most statistical analysis and data visualizations were performed using R v4.0.3 in RStudio v1.3.1093.

The data reported in this paper are deposited in the GEO database under accession number GSE153526 and the Sequence Read Archive database under PRJNA642304.

K. Song reports a patent for methods and compositions for treating melanoma pending. S.R. Dehkordi reports grants from U24CA264379 and R01GM114362 during the conduct of the study. G. Moriceau reports grants from NIH-P01 and the Melanoma Research Alliance during the conduct of the study. P.J. Hergenrother reports that the University of Illinois has a pending patent on compounds described herein, with P.J. Hergenrother listed an inventor. V. Bafna is a cofounder, consultant, and scientific advisory board member of and has equity interest in Boundless Bio, inc. and Abterra, Inc. The terms of this arrangement have been reviewed and approved by the University of California, San Diego, in accordance with its conflict of interest policies. R.S. Lo reports grants from Merck, Pfizer, Bristol Myers Squibb, and OncoSec outside the submitted work. T.G. Graeber reports grants from the NIH, the Melanoma Research Alliance, the W.M. Keck Foundation, the UCLA Eli and Edythe Broad Center of Regenerative Medicine, and Stem Cell Research Hal Gaba Director's Fund for Cancer Stem Cell Research during the conduct of the study; personal fees from Amgen, Auron Therapeutics, Boundless Bio, Coherus BioSciences, and Trethera Corporation, grants from ImmunoActiva, and nonfinancial support from BridgeBio Pharma outside the submitted work; and a patent for methods and compositions for treating melanoma pending. No disclosures were reported by the other authors.

K. Song: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. J.K. Minami: Data curation, validation, investigation, visualization, writing–review and editing. A. Huang: Data curation, formal analysis, validation, investigation, visualization. S.R. Dehkordi: Software, formal analysis, investigation, visualization, methodology, writing–review and editing. S.H. Lomeli: Resources, investigation, writing–review and editing. J. Luebeck: Data curation, software, investigation, visualization, methodology, writing–review and editing. M.H. Goodman: Data curation, formal analysis, and investigation. G. Moriceau: Resources and methodology. O. Krijgsman: Resources, data curation, and methodology. P. Dharanipragada: Data curation, validation, and methodology. T. Ridgley: Data curation, formal analysis, investigation, visualization, and methodology. W.P. Crosson: Resources and investigation. J. Salazar: Data curation, investigation, and visualization. E. Pazol: Data curation, formal analysis, and investigation. G. Karin: Data curation, investigation, and methodology. R. Jayaraman: Data curation and visualization. N.G. Balanis: Resources. S. Alhani: Data curation, investigation, and visualization. K. Sheu: Resources and investigation. J. ten Hoeve: Resources and methodology. A. Palermo: Conceptualization and methodology. S.E. Motika: Resources. T. Senaratne: Resources and methodology. K.H. Paraiso: Writing–review and editing. P.J. Hergenrother: Resources. P. Rao: Resources and methodology. A.S. Multani: Investigation and methodology. D.S. Peeper: Resources. V. Bafna: Software, methodology, writing–review and editing. R.S. Lo: Resources, methodology, writing–review and editing. T.G. Graeber: Conceptualization, resources, supervision, funding acquisition, methodology, project administration, writing–review and editing.

T.G. Graeber is supported by the NIH NCI P01 CA168585 and R01 CA222877, the Melanoma Research Alliance (MRA; Team Science Award 691165), the UCLA SPORE in Prostate Cancer (NIH NCI P50 CA092131), the W.M. Keck Foundation, and the UCLA Eli and Edythe Broad Center of Regenerative Medicine, and Stem Cell Research Hal Gaba Director's Fund for Cancer Stem Cell Research. R.S. Lo is supported by the NIH (1R01CA176111A1, 1R21CA21591001, R21CA25583701, and 1P01CA168585), the MRA (Team Science Award 691165), the V Foundation for Cancer Research (Translational Award), Mary Tanner and Maurizio Grimaldi, and the Ressler Family Foundation. V. Bafna, J. Luebeck, and S.R. Dehkordi were supported in part by grants from the NIH (U24CA264379 and R01GM114362). K.H. Paraiso is a postdoctoral fellow supported by the UCLA Tumor Cell Biology Training Program (USHHS Ruth L. Kirschstein Institutional National Research Service Award # T32 CA009056). G. Moriceau is supported by the NIH (1P01CA168585) and the MRA (Young Investigator Award). P. Dharanipragada and A. Palermo are supported by Jonsson Comprehensive Cancer Center (JCCC) postdoctoral fellowships. FISH microscopy was performed at the Advanced Light Microscopy/Spectroscopy Laboratory and the Leica Microsystems Center of Excellence at the California NanoSystems Institute at UCLA with funding support from NIH Shared Instrumentation Grant S10OD025017 and NSF Major Research Instrumentation grant CHE-0722519. Flow cytometry was performed in the UCLA JCCC and Center for AIDS Research Flow Cytometry Core Facility that is supported by NIH awards P30 CA016042 and 5P30 AI028697, and by the JCCC, the UCLA AIDS Institute, the David Geffen School of Medicine at UCLA, the UCLA Chancellor's Office, and the UCLA Vice Chancellor's Office of Research. Next-generation sequencing was performed at the Technology Center for Genomics and Bioinformatics at UCLA. CGH and low-pass WGS were performed at PacGenomics. Bionano optical mapping assemblies were provided by Andy Pang from Bionano Genomics, Inc. Metabolomics was performed at the UCLA Metabolomics Center supported by the NIH Shared Instrumentation grant S10 OD016387 (T.G. Graeber). A subset of the results shown here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.