Single-Cell Profiling Reveals Metabolic Reprogramming as a Resistance Mechanism in BRAF-Mutated Multiple Myeloma

Precision medicine treatments rarely achieve cures in patients with most types of cancers. Patients with malignant melanoma, for instance, may experience an initial reduction in tumor load but experience disease recurrence later on (1, 2). In multiple myeloma, a hematologic malignancy with oncogenic BRAF mutations in approximately 5% to 12% of patients, recurrent mutations in the BRAF oncogene have been exploited with great therapeutic success (3, 4). In the majority of these patients, however, responses have been transient and competitive outgrowth of clones harboring NRAS or KRAS mutations has been reported in relapsed disease (4, 5). Many patients with BRAF-mutated myeloma who are treated with BRAF inhibitors only experience minimal responses or stable disease in the first place, despite initial absence of activating KRAS or NRAS mutations (6–8). Based on these observations we hypothesized that nongenetic mechanisms of resistance to precision medicine may allow myeloma cells to rapidly adapt to therapeutic challenges.

Using BRAF-mutated myeloma as a model disease in this study, we aimed to investigate why “precision treatment” of BRAF-mutated neoplasms with BRAF and/or MEK inhibition often fails to effectively elicit substantial treatment responses. We focused on BRAF-mutated myeloma, as sensitive liquid biopsy technology allows us to follow the cancer-cell fate over time by interrogating circulating tumor cells from patients (9).

Cellular plasticity has recently been identified as an alternative mechanism driving nongenetic resistance as cancer cells acquire transcriptional states which no longer depend on the drug target (10). A better understanding of how these nongenetic programs integrate with genetic resistance may help with designing more effective therapies to overcome resistance to targeted therapies. We hypothesized that adaptive resistance mechanisms allow myeloma cells to rapidly alter their state, thus compensating for the loss of oncogenic stimulation. Utilizing single-cell RNA (scRNA) sequencing (scRNA-seq) of myeloma cells we explored how inhibition of the BRAF/MEK pathway produces transcriptional state changes and exposes putative therapeutic vulnerabilities. We defined differential enhancer usage in cells that have the ability to persist in the presence of BRAF inhibitor treatment. We explored how transcriptional state changes mediate metabolic reprogramming and selective activation of genes involved in metabolic regulation. We defined the kinetics of state changes and determined if susceptibility to pharmacologic inhibition of oxidative phosphorylation (OxPhos) may represent an opportunity to overcome drug resistance. These insights into the mechanisms that allow malignant plasma cells to escape precision treatments will enable the design of more effective targeted treatments in myeloma and other BRAF-mutated neoplasms.

Samples were collected from patients with relapsed/refractory multiple myeloma (RRMM) harboring activating BRAF mutations (Table 1) who provided informed written consent and were treated according to the Institutional Review Board (IRB)–approved Dana-Farber Cancer Institute (DFCI) protocol 16–352 in accordance with the Declaration of Helsinki. Patients were treated with the BRAF inhibitor dabrafenib and the MEK inhibitor trametinib. Samples were retrieved at screening, after 1 week of treatment and at the time of clinical progression. Clinical progression for patients MM1, MM2, and MM3 was measured at cycle 4 day 1 (MM1), week 1 (MM2), and cycle 7 day 1 (MM3). In brief, dabrafenib was administered orally twice a day and trametinib was administered orally every day on a 28-day cycle until progression, unacceptable toxicity, patient refusal, or changes in the patient's condition that prohibited further treatment. Normal donor plasma cells were obtained commercially from Research Blood Components, LLC.

Plasma cells were separated using the EasySep Human CD138 Positive Selection Kit (StemCell Technologies). Samples were next stained with anti–CD38-FITC (multi-epitope), anti–CD138-PE (clone 44F9), anti–CD319-APC (clone 162.1), anti–CD269-PE/Cyanine7 (clone 19F2), and DAPI according to the manufacturers' recommendations (Supplementary Table S1). Single cells were sorted into 96-well plates containing TCL buffer (Qiagen) on a Sony SH800 sorter. Data were analyzed using the FlowJo software version 10.0 (Becton Dickinson).

Full-length scRNA-seq libraries were prepared using Smart-seq2 (11). cDNA was fragmented by Nextera XT and amplified with indexed Nextera PCR primers. Products were purified with Agencourt Ampure XP beads (Beckman Coulter) and quantified using a Bioanalyzer High Sensitivity DNA Kit. Pooled sequencing of Nextera libraries was carried out using a NextSeq 500 (Illumina, C>A). All sequenced reads were paired-end and 36-bp read length. A sequencing depth of 1 × 10e6 reads was aimed for in each cell.

U266 was obtained from DMSZ and DP6 was provided by Dr. D. Jelinek (Mayo Clinic). DP6 cells were maintained in Advanced RPMI 1640 with 4% FBS (not heat-inactivated), 1x Glutamax, and IL6 supplementation (2 ng/mL QOD). Effective inhibition of the MAPK pathway in vitro was established at an IC₉₀ dose of 1 μmol/L dabrafenib for U266 and 10 nmol/L dabrafenib for DP6. To account for genetic heterogeneity at baseline, single-cell clones were generated by single-cell sorting using a Sony SH800 sorter. All single-cell clones were tested and found to be negative for mycoplasma contamination. Authentication was performed using DNA fingerprinting with small tandem repeat (STR) profiling.

Single-cell clones were cultured with dabrafenib-containing media replaced twice per week. Aliquots of 5 × 10e5 cells were retrieved at defined time points throughout drug exposure for immunoblotting, low-pass whole-genome sequencing (LPWGS), and bulk RNA-seq. For Western blot analysis, antibodies were used as indicated in Supplementary Table S2. For bulk RNA-seq, cells were resuspended in trizol and RNA was extracted using the Qiagen RNeasy formalin-fixed, paraffin-embedded (FFPE) kit. Libraries were prepared with the TruSeq RNA Access Library Prep Kit (Illumina). Paired-end sequencing, using 75 bp reads was performed on a NextSeq 500 (Illumina). For LPWGS, genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen) and genomic DNA was quantified using the Qubit HS DNA kit (Invitrogen). An input of 5 to 20 ng DNA was used for library construction with custom adapters (IDT and Broad Institute, Cambridge, MA). Libraries were pooled and sequenced to an average sequencing depth of 0.2x.

DMSO-diluted dabrafenib (ThermoFisher), rotenone (Sigma Aldrich), and IACS-010759 (BioVision) were used. Cell viability was estimated using the CellTiter-Glo assay (Promega) on a SpectraMax M5 reader (Molecular Devices). If not stated otherwise, control cells were incubated for the same duration with equal amounts of DMSO and luminescence was normalized to the respective controls.

For chromatin immunoprecipitation sequencing (ChIP-seq), 2 × 10e7 cells were cross-linked using formaldehyde and sonicated on a Covaris E220 focused ultrasonicator (approximately 200–500 bp fragment size). Immunoprecipitation (IP) was performed with anti-H3K27ac antibody (Active Motif) and parts of the lysate were saved as whole cell extract (WCE). Drosophila chromatin and antibody spike-ins were done for normalization. IP pulldown was performed using Protein A/G dynabeads (Invitrogen), followed by washing and elution off the beads. Cross-linking was next reversed by proteinase K. Libraries were constructed with individual barcodes (IP and WCE).

Cells were seeded overnight in 24-well XF microplates (Seahorse Bioscience) at a density of 1.5 × 10⁵ cells per well (200 μL volume). The following day, media was replaced by XF Cell Mito Stress Test Assay Medium (Seahorse Bioscience) and adjusted to a pH of 7.4. Cells were next incubated for one hour at 37°C without CO₂. Injection ports were loaded with oligomycin (1 μmol/L), carbonyl cyanide-4-(trifluoromethoxy) phenylhydrazone (FCCP, 1 μmol/L+1 μmol/L), and rotenone/antimycin A (0.5 μmol/L each). Each drug was added sequentially using three basal rate measurements for oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) prior to the first injection followed by three measurements after each injection. Intervals between these measurements were 3/2/3 minutes.

Processing and quality filtering of scRNA-seq data was carried out using trimmomatic, STAR, HTSeq, and RSEM. Four different parameters were used to filter out low-quality cells: (i) library size, (ii) number of genes detected, (iii) percentage of reads mapping to mitochondrial genes, and (iv) percentage of reads mapping to house-keeping genes (Supplementary Figs. S1 and S2). Cells that fell below three median absolute deviations were excluded from downstream analyses. As an additional approach to remove bad-quality cells without any predefined cutoffs we used the R package mvoutlier. This resulted in the retention of 1,153 cells with an average of 4,836 genes detected per cell.

Clustering and marker gene analyses were performed using PAGODA2. In order to rule out the possibility of clusters being purely driven by cell cycle, each individual cell was additionally analyzed for expression of G1, G2M, and S-phase markers using SEURAT2. The marker genes for each cluster were determined using the findMarker function in the scran package. Copy-number analysis was carried out with InferCNV (https://github.com/broadinstitute/inferCNV). Untreated normal donor plasma cells served as control to identify malignant plasma cells based on copy-number alterations (CNA) and transcriptional profile and were used for quality control in our analytical pipeline.

Output for bulk RNA-seq was generated in FASTQ format (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) mapped to the human genome GRCh37 (hg19). Gene expression at transcript-level resolution was calculated using RSEM (v1.2.31) and Picard (https://broadinstitute.github.io/picard/). Annotations were derived from the Ensembl database (October 2018).

Differential gene expression analysis and enrichment were performed using DESeq2 (v1.22.2) and gene set enrichment analysis (GSEA; v4.0.3). Detection of CNAs was carried out using HMMcopy (v1.28.0).

H3K27ac chromatin immunoprecipitation (ChIP) data was analyzed using the nf-core/chipseq pipeline as previously described (12). In brief, quality filtering was done by FastQC and TrimGalore. Reads were aligned to UCSC_hg19 using Burrows–Wheeler Aligner (BWA) and SAMtools. IP efficiency was calculated by normalizing peak counts to their WCE using deepTools. Narrow consensus peaks were called using MACS, HOMER, BEDTools, and featureCounts. Enhancer characterization was done by ChIPseeker (v1.20.0; ref. 13). Gene-enhancer relationships were annotated with GREAT. GSEA was performed using GSEA (v4.0.3) and was visualized in IGV (v2.5.0).

Sequencing data have been deposited at Gene Expression Omnibus (GEO) database (reference number GSE168951). Codes used for all analyses will be made available upon request.

We hypothesized that cancer cells are capable of rapidly adjusting to therapeutic pressure when challenged with targeted inhibitors. We previously described a methodology to isolate myeloma cells from blood or bone marrow at multiple time points during the course of the disease with high sensitivity (9). We used this approach to define transcriptional state changes in patients with BRAF-mutated myeloma treated with the BRAF inhibitor dabrafenib and the MEK inhibitor trametinib. A total of 1,495 single CD138+/CD38+ myeloma and normal plasma cells from 3 patients and 1 normal donor were examined by full-length scRNA-seq (Table 1). All 3 myeloma patients had received at least two prior therapeutic regimens. One patient harbored a BRAF V600E and 2 patients had D594N mutations (Table 1). Two of these patients also harbored additional NRAS Q61 mutations, confirmed by next-generation sequencing. All patients showed varying degrees of clinical response to dabrafenib/trametinib treatment (Table 1). We isolated myeloma cells from the bone marrow or peripheral blood of these patients before and at various time points after the treatment was started (Supplementary Fig. S3). After employing quality control filtering steps, 1,153 single cells remained for downstream analyses. We performed PAGODA2 and t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization at screening, for the earliest available time point on treatment (1 week) and the time of clinical relapse. Our analysis revealed 8 distinct clusters (Fig. 1A). One cluster comprised plasma cells from the normal donor, whereas the remaining clusters reflected the heterogeneity of myeloma and segregated according to the individual patients MM1, MM2, and MM3 (Fig. 1A, left). Single myeloma cells that were obtained from patients at the time of clinical relapse clustered distinctly from myeloma cells that were obtained prior to initiation of treatment. Notably, transcriptional state changes were already evident as early as 7 days after patients had been started on treatment with dabrafenib/trametinib (Fig. 1A, right).

We postulated that the rapid transcriptional adaptation that we observed allows myeloma cells to compensate for the loss of oncogenic stimulation and to persist in the presence of BRAF/MEK inhibitor treatment. First, we ensured correct classification of myeloma cells and normal nonmalignant hematopoietic cells, by determining large-segment CNAs and expression of translocation target genes in all single cells (Fig. 1B). Pathogenic CNA in single cells were consistent with aberrations as detected by clinical routine cytogenetics, including del13q in MM1 and MM3, del1p in MM2, and del8p in MM3 (Table 1). As expected, differential gene-expression analysis of MAPK pathway genes showed various degrees of signaling pathway deactivation in MM1, MM2, and MM3 at the time of relapse to BRAF/MEK inhibitor treatment (Fig. 1C).

While MAPK-related genes were downregulated with varying kinetics in the individual patients, heterogeneous programs that have previously been linked to myeloma drug resistance were upregulated at the same time (Fig. 1D). In MM1, persisting cells showed a transient increase in cellular stress response (NEAT1; refs. 14, 15) and glutamine dependency (SLC1A5; refs. 14, 15), which was no longer detectable at the time of clinical relapse. At relapse, this patient showed significant downregulation of BCL2L11 coding for the tumor suppressor Bcl2-interacting mediator of cell death (BIM) (Supplementary Fig. S4), a BCL-S2 family member which has previously been defined as a central mediator for BRAF and MEK inhibitor efficacy in myeloma (16). Persistent cells in MM2 showed upregulation of genes involved in epigenetic modulation (HDAC3), tyrosine kinase signaling (INSR), regulation of cullin-ring ligases (NEDD8), and triggering of apoptosis (LGALS3; refs. 17–20). In addition, persistent cells in MM2 showed increased expression of potential immunotherapeutic targets such as SDC1 and TNFRSF17 coding for the myeloma surface markers CD138 and BCMA, respectively (21, 22).

The transcriptional state changes in persistent myeloma cells of MM3 differed markedly from MM1 and MM2 and included upregulation of the transcription factors IRF1, STAT1, and STAT3 (23–25). These data demonstrate that myeloma cells show rapid adaptation to therapeutic challenges in vivo with detectable transcriptional state changes as early as 1 week after BRAF/MEK inhibition, suggesting that these state changes may occur before genomic evolution. These observations also raised the question if transcriptional state changes generate potential therapeutically exploitable vulnerabilities.

We next asked how transcriptional state changes relate to clonal selection as BRAF inhibitor resistance develops in BRAF-mutated myeloma. Since the characterization of clonal evolution in patients is confounded by extensive spatial heterogeneity and sampling bias, we established an in vitro model of BRAF inhibitor resistance. To this end we generated 3 single-cell clones from U266, a BRAF^K601N-mutated myeloma cell line to optimize genetic homogeneity in the input population of myeloma cells (Fig. 2A). This approach provided a genetically homogeneous input population and allowed for testing of multiple replicates, given that currently only the two myeloma cell lines U266 and DP6 have been described to harbor BRAF mutations (23, 26). Effective inhibition of the MAPK pathway was established at an averaged IC₉₀ dose of 1 μmol/L dabrafenib (Fig. 2B and C).

All 3 clones were next subjected to long-term dabrafenib treatment. The enumeration of viable cells followed a similar pattern in all U266 clones with a sharp decrease in viability, a plateau phase of persistence beginning at approximately day 7, and outgrowth of resistant myeloma cells at approximately day 60 of treatment (Fig. 2A). Serial profiling by immunoblotting over time demonstrated late recovery of MAPK pathway activity (Fig. 2D; Supplementary Fig. S5). As expected, dabrafenib sensitivity was lower in resistant clones as compared with their respective untreated controls after 72 hours of treatment (Fig. 2E; Supplementary Fig. S6). In all clones, resistance was associated with the occurrence of heterogeneous CNA at these late time points as detectable by LPWGS (Fig. 2F; Supplementary Fig. S7), indicating evolving competitive outgrowth of genetically distinct myeloma cells at subclonal level. Notably, although all three replicate U266 clones had identical copy-number profiles initially, the clones that were competitively selected after 2 months as BRAF inhibitor resistance had developed all had distinct CNA. These data demonstrate that the late phase of dabrafenib resistance is characterized by outgrowth of multiple drug-resistant clones, all of which are genetically distinct.

To determine if similar genetic heterogeneity was already present during the early phase of drug resistance to BRAF inhibition, we performed LPWGS of BRAF-mutated U266 myeloma cells after 7 days of BRAF inhibition. In contrast to the late phase of drug resistance we found identical CNA for all 3 U266 single-cell clones at baseline and persistent stage, indicating that this early persistence is not driven by clonal selection but by transcriptional plasticity of myeloma cells (Fig. 3A; Supplementary Fig. S8).

To further characterize early persistent cells, we performed bulk RNA-seq and examined how transcriptional states in U266 and DP6, a second myeloma cell line with known BRAF^V600E mutation (27), at baseline diverge from persistent cells after 7 days of dabrafenib exposure, respectively (Fig. 3B; Supplementary Figs. S9 and S10). Interestingly, even though the genomics of resistant U266 clones ultimately differed substantially at the time of competitive outgrowth, the transcriptional changes after 7 days were homogeneous (Fig. 3C). This was similar for DP6-derived single-cell clones, although differentially regulated genes were distinct between U266 and DP6. To explore if any differentially regulated transcriptional signatures were shared across both cell lines and all single-cell clones, we performed GSEA. As expected, BRAF inhibition led to downregulation of MAPK pathway perturbation genes in U266 and DP6 (Fig. 3D), whereas OxPhos emerged as the most consistently enriched signaling pathway in persistent cells as compared with baseline (P = 0.0001, FDR = 0.0001; Fig. 4A and B; Supplementary Figs. S11 and S12). These data suggest that BRAF inhibition in BRAF-mutated myeloma cells leads to transcriptional reprogramming with induction of OxPhos-related genes within a brief period of time.

We next examined the mechanisms of transcriptional adaptation in persistent cells and asked if transcriptional state changes can be attributed to alterations in enhancer activity. To this end, we performed ChIP-seq for H3K27ac as a key enhancer mark in 3 U266 single-cell clones at baseline versus persistent stage and asked to what extent H3K27ac peaks were distinct between naïve and persistent cells, and which cellular functions might be affected by differential epigenetic regulation.

We observed substantial remodeling illustrated by a greater number of unique H3K27ac peaks in dabrafenib-persistent cells as compared with dabrafenib-naïve controls. Out of a total of 22,756 merged peaks, which were called in at least three samples, 7,645 (34%) were unique to persistent cells compared with 454 (2%) in untreated cells (Fig. 4C, left), illustrating greater epigenetic activity being associated with transcriptional adaptation and persistence to BRAF inhibition. While the majority of overlapping peaks was found in promoter regions, there was a significant change to distal regulatory elements, which were exclusively detectable in clones at the persistent state (Fig. 4C and D). As H3K27 acetylation of histones in these regions is typically associated with enhancer activity, we determined the number of those putative enhancer peaks. Notably, 5,240/11,756 (45%) of these peaks were specific for persistent cells, with the majority of those peaks mapping to intronic and distal intergenic regions, while only 388/11,756 (3%) of the enhancer peaks were specific for cells at baseline (Fig. 4C, right). We performed a gene set analysis for peak intensities of putative enhancers focusing on intergenic peaks and their corresponding genes as predicted by GREAT (Fig. 4E; ref. 28). This analysis demonstrated that the OxPhos pathway signature was enriched in persistent cells as third-most enriched from a total of 186 gene sets for intergenic enhancers (n = 3,236 peaks). For example, H3K27ac signal intensities were increased at the enhancer sites of the OxPhos-related gene SDHB, which has been described to stimulate OxPhos via HIF-1α degradation and stabilization of mitochondrial cytochrome oxidase IV (COX4–1; refs. 29, 30), and at enhancers near ATP6V1E1, coding for a V-ATPase V1 subunit which has been shown to regulate senescence by promoting lysosomal acidification, metabolic reprogramming, and mitochondrial recovery in yeast (Fig. 4F; ref. 31). These data therefore indicate that differential enhancer usage is associated with transcriptional adaptation to selective inhibition of the BRAF oncogenic pathway with induction of genes related to OxPhos.

We next asked if the observed transcriptional shifts toward OxPhos translate into functional dependencies. To assess mitochondrial metabolism, we performed metabolic flux analysis in dabrafenib-persistent U266 and DP6 single-cell clones after 14 days of drug exposure. OCR and ECAR were measured as indicators of OxPhos and glycolysis, respectively. Sequential addition of the ATP synthase inhibitor oligomycin, the mitochondrial uncoupler FCCP, and the complex III and complex I inhibitors antimycin A and rotenone enabled measurement of the largest possible capacity for OxPhos under maximal stress conditions (Fig. 5A). The baseline ratio of OCR:ECAR indicating a preference for OxPhos over glycolysis was significantly elevated in dabrafenib-persistent cells compared with DMSO controls (Fig. 5B and C; Supplementary Fig. S13A and S13C). Correspondingly, the maximal OCR:ECAR ratio under FCCP stimulation was significantly elevated in dabrafenib-persistent cells compared with controls (Supplementary Fig. S13B and S13D), demonstrating that while responding to cellular stress, dabrafenib-persistent myeloma cells rather maintain OxPhos than glycolysis as a preferential resource to recruit energetic metabolites.

Next, we asked if this functionally validated OxPhos dependency could be therapeutically exploited to target dabrafenib persistence in vitro. To this end we generated dabrafenib-persistent cells by exposing U266 and DP6 single-cell clones to 1 μmol/L and 10 nm dabrafenib as previously established for both cell lines. At day 14, persistent cells and corresponding DMSO controls were washed, purified by density gradient centrifugation, and sequentially treated for 7 days with increasing doses of rotenone or IACS-010759, a clinical-grade small-molecule inhibitor of complex I of the mitochondrial electron transport chain. All dabrafenib-persistent U266 and DP6 single-cell clones demonstrated greater sensitivity to rotenone or IACS-010759 than their respective controls (Fig. 5D; Supplementary Figs. S14 and S15)

These data demonstrate that the transcriptional state changes and enhancer rewiring that occur with BRAF inhibitor resistance may create metabolic dependencies that can be targeted subsequently as therapeutic vulnerabilities.

To determine if the induction of OxPhos-related genes can also be found after inhibition of BRAF and MEK in vivo, we defined the transcriptional regulation of OxPhos-related genes in the single-cell dataset from patients with BRAF-mutated myeloma (Fig. 1). Equivalent to the transcriptional regulation in myeloma cell lines, OxPhos-related genes were upregulated in response to inhibition of the BRAF/MEK pathway in patients with myeloma (Fig. 5E), which was inversely correlated with MAPK-associated transcriptional profiles (Fig. 1D). While gene set analysis demonstrated significant activation of OxPhos programs in all patients at single-cell level, the kinetics of this activation differed between individuals (Fig. 5E). Moreover, the magnitude of OxPhos induction varied, suggesting that partial OxPhos induction may also be a consequence of resistance to other myeloma treatments, given that the patients received a number of different medications prior to BRAF/MEK inhibition (Table 1). Taken together these data demonstrate that successful inhibition of the BRAF/MEK pathway in patients with myeloma is accompanied with a state of metabolic resistance characterized by preferential transcriptional activation of OxPhos.

Precision medicine typically refers to targeted pharmacologic inhibition of oncogenes with well-defined mutations. While such approaches have proven to be a successful therapeutic concept in various neoplasms, they are rarely curative (32–34). BRAF mutations occur in 5% to 12% of patients with myeloma and have attracted considerable attention because of their therapeutic potential for targeted inhibition (3, 4, 35). Although remarkably effective in some patients, BRAF inhibitors typically induce only transient responses despite initial absence of established resistance-conferring mutations (6–8, 36). Using BRAF-mutated myeloma as a model for resistance to precision medicine we investigated if BRAF-mutated myeloma cells have the ability to ensure their survival by rapidly adapting to BRAF-inhibitor treatment.

Cellular plasticity and dedifferentiation have recently emerged as a principle of nongenetic resistance in cancer, as cells acquire transcriptional states which no longer depend on the drug target (10). Here, we demonstrate that myeloma cells undergo transcriptional state changes when exposed to BRAF-inhibitor treatment in vitro and in patients. While clonal evolution and outgrowth of genetically distinct clones become apparent at late time points, transcriptional adaptation occurs very rapidly, and can be observed within one week of treatment when clonal selection has not yet commenced. Wide-spread enhancer recruitment appears to facilitate dynamic switching between transcriptional cell states, as emergence of transcriptionally altered myeloma cells coincides with significant changes in chromatin regulation. Epigenetic remodeling has previously been shown to be provoked by a diverse array of stimuli including macrophage activation, inflammation, atherogenesis, and kinase inhibition (37–39). In this study, we demonstrate that BRAF inhibition is associated with differential enhancer usage and greater transcriptional promiscuity, likely in order to permit wider sampling of the transcriptome. These findings are in line with a previous report which indicates that transcriptional plasticity in myeloma is driven by transregulated epigenetic modulation and develops largely independent of predisposing genetics (40). The mechanisms that control permissive chromatin states and transcriptional adaptation remain interesting subjects of future investigations. In this study, pharmacologic inhibition of the oncogenic BRAF/MEK pathway led to induction of alternative survival signals, which may have been conveyed either by other oncogenes, by environmental regulatory pathways through cytokines and surface receptors, or by hardwired survival programs established for normal plasma cells.

Mitochondrial reprogramming has previously been described to promote metabolic recovery and attenuate senescence in fruit flies and yeast (31). OxPhos is normally utilized as the main pathway for energy breakdown in aerobic organisms and has raised interest as an emerging targetable pathway in several types of cancer: studies in melanoma have demonstrated that BRAF-mutant cells, which develop resistance to BRAF inhibitors, display increased activity of oxidative metabolism, increased dependency on mitochondria, and higher levels of reactive oxygen species (41–43). OxPhos activation has also been linked to resistance to the Bruton's tyrosine kinase inhibitor ibrutinib in mantle cell lymphoma (44) and has been observed in therapy-resistant myeloid leukemia stem cells (45, 46).

In myeloma cell lines targeting of glycolysis and OxPhos by repurposing the FDA-approved drugs ritonavir and metformin has cytotoxic activity (47), and OxPhos overexpression has previously been associated with resistance to the proteasome inhibitor bortezomib (48). Pharmacologic inhibition of the transcriptional coactivator PGC-1α in myeloma cell lines has been reported to result in reduced proliferation, G₂–M arrest, and impaired OxPhos (49). Taken together, these data indicate that the OxPhos pathway may represent a putative target in patients with myeloma but its regulation and the patient population that might benefit from OxPhos inhibition have been unclear thus far.

Our data provide understanding why OxPhos inhibitors may not be effective as single agents but require combinatory drugs to induce metabolic vulnerability. By performing metabolic flux analysis we demonstrate that functional dependency on mitochondrial ATP generation arises as myeloma cells become resistant to the BRAF inhibitor dabrafenib. Although OxPhos may be difficult to be targeted therapeutically, agents that exploit metabolic dependencies in cancer have been introduced. IACS-010759, a mitochondrial ETC complex I inhibitor, showed promising in vitro efficacy in this study and is currently being investigated as single agent in phase I clinical trials for advanced tumors, including triple-negative breast cancer, pancreatic adenocarcinoma, recurrent lymphoma, and acute myeloid leukemia (NCT03291938, NCT02882321).

Upfront combination treatments rather than sequential administration may be a more effective path toward early killing of myeloma cells, which would also slow down clonal evolution, simply by reducing the number of remaining cells that can continue to divide and acquire new genomic aberrations. More studies will be needed to determine the best possible sequence or combinations of individual lines of therapy, including novel immunotherapeutic approaches and stem-cell transplantation, which almost all patients with myeloma require. Interestingly, in 1 patient (MM2) BRAF treatment induced upregulation of TNFRSF17 (BCMA), which is a clinically actionable immunotherapeutic target, suggesting that precision-medicine treatment may unexpectedly sustain immunotherapeutic target expression in some patients. More studies are needed to identify which patients under which circumstances might benefit from this approach.

Although we observed that OxPhos activation is shared as a metabolic resistance program across patients with myeloma and cell lines, the majority of activated pathways appeared to be proprietary to the individual patient or cell line. Indeed, our findings indicate that dependency on the BRAF/MEK pathway may be greater in some patients than in others. A deeper understanding on the relative contribution of different oncogenic pathways that drive myeloma growth in individual patients would aid to identify those pathways that should be prioritized for maximum therapeutic benefit.

Overcoming the apparent heterogeneity of alternate cell states and obvious diversity of regulatory programs between patients represents a therapeutic challenge and would argue for repeated molecular characterization of individual patients to identify resistance pathways, for example by liquid biopsy (9, 50). Epigenetic drugs could potentially be used to inhibit epigenetic remodeling and subsequent metabolic reprogramming, especially since clinical efficacy of HDAC inhibitors has already been established in myeloma with panobinostat being FDA approved for the treatment of double-refractory disease.

In summary, we demonstrate that differential enhancer usage is associated with transcriptional plasticity in BRAF-mutated myeloma cells that are exposed to BRAF/MEK inhibition. This transcriptional adaptation is rapid, precedes detectable outgrowth of genetically discernible resistance variants, and promotes metabolic reprogramming before genetic alterations provide additional survival advantages. While the majority of transcriptional state changes that occur with BRAF/MEK inhibition appear to be patient specific, preferential utilization of oxidative phosphorylation appears to represent an inducible vulnerability across patients with drug-resistant myeloma and cell lines. We propose that our findings inform successful implementation of precision medicine in BRAF-mutated neoplasms beyond myeloma. Similar mechanisms may also extend to patients with malignancies that harbor other targetable oncogenes and to a broader population of patients with myeloma irrespective of their BRAF/MEK mutation status.

Therapeutic efficacy will require rational combination therapy as well as frequent molecular characterization, for example by liquid biopsy, to accurately identify resistance pathways and therapeutic escape (9, 50). Targeting the interplay of epigenetic and metabolic adaptation as a basis of nongenetic plasticity may help to overcome refractory disease in myeloma and other BRAF-mutated cancers (6, 51–54).

J.M. Waldschmidt reports personal fees from Sanofi and Janssen outside the submitted work. C. Grassberger reports serving on the Scientific Advisory Board for Nanolive. A.J. Yee reports personal fees and non-financial support from Adaptive, Amgen, BMS, Celgene, Janssen, and Takeda, as well as personal fees from Karyopharm, Regeneron, Sanofi, and GSK outside the submitted work. P.G. Richardson reports personal fees from AbbVie, AstraZeneca, GSK, Janssen, Protocol Intelligence, Regeneron, Sanofi, and Secura Bio, as well as grants and personal fees from Celgene/BMS, Karyopharm, Oncopeptides, and Takeda outside the submitted work. K.C. Anderson reports personal fees from Janssen, Amgen, Pfizer, AstraZeneca, Precision Biosciences, BMS, Mana, Starton, and Window during the conduct of the study, as well as other support from C4 Therapeutics, Oncopep, Raqia, and NextRNA outside the submitted work. N.S. Raje reports personal fees from Amgen, BMS, Takeda, Janssen, Pfizer, Caribou, Immuneel, and GSK, as well as grants and personal fees from Bluebird outside the submitted work. B. Knoechel reports grants from NCI during the conduct of the study. J.G. Lohr reports grants from NCI, V Foundation, and Anna Fuller Fund during the conduct of the study; J.G. Lohr also reports grants from Bristol Myers Squibb, as well as personal fees from T2 Biosystems outside the submitted work. No disclosures were reported by the other authors.

J.M. Waldschmidt: Conceptualization, data curation, software, formal analysis, funding acquisition, validation, investigation, methodology, writing–original draft, writing–review and editing. J.A. Kloeber: Conceptualization, resources, data curation, software, formal analysis, investigation, methodology. P. Anand: Conceptualization, resources, data curation, software, formal analysis, validation. J. Frede: Resources, data curation, investigation, methodology. A. Kokkalis: Resources, validation, investigation, methodology. V. Dimitrova: Resources, methodology. S. Potdar: Resources, methodology. M.S. Nair: Resources, data curation, software, formal analysis. T. Vijaykumar: Resources, methodology. N.G. Im: Validation, methodology. A. Guillaumet-Adkins: Validation, methodology. N. Chopra: Visualization, methodology. H. Stuart: Data curation, software, formal analysis, validation. L. Budano: Resources, data curation. N. Sotudeh: Data curation, software, formal analysis, validation. G. Guo: Resources, data curation, software, formal analysis, validation. C. Grassberger: Visualization, methodology. A.J. Yee: Conceptualization, resources, data curation, project administration. J.P. Laubach: Conceptualization, resources, data curation, project administration. P.G. Richardson: Conceptualization, resources, data curation, project administration. K.C. Anderson: Conceptualization, resources, data curation, project administration. N.S. Raje: Conceptualization, resources, data curation, project administration, writing–review and editing. B. Knoechel: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. J.G. Lohr: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

J.M. Waldschmidt is supported by a postdoctoral fellowship of Deutsche Forschungsgemeinschaft (German Research Foundation, 391926441) and the International Myeloma Society Young Investigator Award (IMW, Boston, 2019). J.G. Lohr is supported by the NCI (K08CA191026), the V Foundation for Cancer Research and the Anna Fuller Fund. B. Knoechel is supported by the NCI (K08CA191091). We are grateful to D. Jelinek for generously providing the DP6 multiple myeloma cell line. I.C. Newbert and M. Young provided clinical information for correlative studies. J. Eismann and O. Sonzogni provided expertise for metabolic assays. W. Oldham and F.P. Bowman provided experimental expertise and kindly assisted with all Seahorse experiments. M. Wal, H. Yun, A.J. Rogers, and J. Poller discussed data and provided methodologic resources throughout the study.

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.