Epigenetic modulation reveals differentiation state specificity of oncogene addiction

Hyperactivation of the MAPK signaling pathway motivates the clinical use of MAPK inhibitors for BRAF-mutant melanomas. Heterogeneity in differentiation state due to epigenetic plasticity, however, results in cell-to-cell variability in the state of MAPK dependency, diminishing the efficacy of MAPK inhibitors. To identify key regulators of such variability, we screen 276 epigenetic-modifying compounds, individually or combined with MAPK inhibitors, across genetically diverse and isogenic populations of melanoma cells. Following single-cell analysis and multivariate modeling, we identify three classes of epigenetic inhibitors that target distinct epigenetic states associated with either one of the lysine-specific histone demethylases Kdm1a or Kdm4b, or BET bromodomain proteins. While melanocytes remain insensitive, the anti-tumor efficacy of each inhibitor is predicted based on melanoma cells’ differentiation state and MAPK activity. Our systems pharmacology approach highlights a path toward identifying actionable epigenetic factors that extend the BRAF oncogene addiction paradigm on the basis of tumor cell differentiation state.

T herapeutic inhibition of oncogenic signaling often leads to variable responses due to cell-to-cell heterogeneity in the state of oncogene dependency 1 . Heterogeneity may result from secondary genetic mutations, or may be caused by epigenetic differences associated with a cell's developmental lineage or differentiation state [2][3][4][5][6][7][8][9][10] . An example of such epigenetic heterogeneity is observed in BRAF-mutated melanomas, causing fractional responses to Braf/Mek-targeted therapies [11][12][13][14][15][16][17][18][19][20][21][22][23] . Numerous studies have associated fluctuations in the state of MAPK dependency to melanoma differentiation state plasticity 6,18,[24][25][26][27] . Such plasticity spans a spectrum, ranging from a pigmented melanocytic phenotype associated with transcriptional regulators Sox10 and Mitf 28 , to a neural crest-like state that expresses Ngfr 12,29 , to an undifferentiated state characterized by high expression of receptor tyrosine kinases such as Axl 19,30 . Melanoma tumors consist of a mixture of these phenotypes at variable single-cell frequencies 25 . Although the consequences of such heterogeneities for drug resistance are widely recognized, there is still little known about actionable epigenetic factors that may link heterogeneity in the state of MAPK dependency to plasticity in differentiation state. A system-wide exploration of these factors may reveal novel opportunities for epigenetic treatments that will minimize the emergence of drug resistance. In principle, identifying epigenetic treatments that promote cellular requirement for MAPK signaling may enhance the efficacy of Braf/Mek inhibitors when used in combination 12,[31][32][33] . Alternatively, treatments that target synthetic lethal partners of the BRAF oncogene may serve as a strategy to kill intrinsically drug-resistant cells 34 .
In this paper, we take a systems pharmacology approach to test the hypothesis that heterogeneity in the state of MAPK dependency may result from a subset of key epigenetic variations across tumor cells of heterogeneous differentiation states. To identify the regulators of such variations, we use a library of 276 smallmolecule epigenetic modulators in BRAF-mutant melanoma cell lines that cover a wide spectrum of differentiation states. We evaluate the effect of each epigenetic modulator (individually and in combination with Braf/Mek inhibitors) on cell survival, proliferation, differentiation state, and MAPK activity. Integrating multiplexed single-cell analysis with multivariate modeling and genetic experiments, we identify three classes of small molecules that target seemingly distinct epigenetic states in melanoma cells. These states include: (i) a lysine demethylase 1A (Kdm1a)dependent state, predominantly observed in undifferentiated cells, that is efficiently inhibitable by the reversible Kdm1a inhibitor SP2509, (ii) a lysine demethylase 4B (Kdm4b)-dependent state, associated with neural crest-like cells, that is sensitive to JIB-04 (a pan-inhibitor of Jumonji histone demethylases), and (iii) a state induced by BET bromodomain inhibitors such as OTX015 (Birabresib), which enhances cells' requirement for MAPK signaling. Single-cell analysis shows that these states might co-exist in different combinations and frequencies, highlighting mutual epigenetic vulnerabilities among genetically diverse melanoma cell populations. Importantly, non-transformed primary melanocytes are not sensitive to these inhibitors. These results, therefore, provide a path for identifying actionable epigenetic factors that may extend the BRAF oncogene addiction paradigm on the basis of tumor cell differentiation state.

Results
Single-cell analysis uncovers heterogeneities in differentiation, proliferation, and signaling states. To elucidate single-cell heterogeneities in differentiation state and their variations across melanoma cell populations, we utilized high-throughput, multiplex immunofluorescence microscopy. We exposed a panel of 16 BRAF V600E/D melanoma cell lines and a batch of non-transformed human primary epidermal melanocytes to the Braf inhibitor vemurafenib (at 100 nM), alone or in combination with the Mek inhibitor trametinib (at 10 nM). For the purpose of comparison, we also included an NRAS Q61K -mutated variant of the A375 cell line, representing a common mechanism of acquired resistance to Braf/Mek inhibitors 35 . All cells were fixed following 3-5 days of treatment and protein levels of three validated differentiation state markers, Mitf, Ngfr, and Axl (Supplementary Fig. 1), were quantified at a single-cell level ( Supplementary Figs. 2, 3). To visualize baseline and treatmentinduced variations in all three markers, we performed t-distributed stochastic neighbor embedding (t-SNE) analysis on a total population of 6069 randomly selected cells from all of the 102 tested conditions, covering the entire panel of cell lines, drugs and timepoints (Fig. 1a, Supplementary Fig. 4a, b). Single-cell analysis revealed a continuum of differentiation states ranging from Mitf High to Ngfr High to Axl High cells. We then used the t-SNE map as a reference to visualize differentiation state variations in isogenic cell populations from each cell line (Fig. 1b, Supplementary Fig. 4c, d). We also quantified the extent of heterogeneity in each marker by computing the Fano factor, a standardized measure of dispersion of probability distribution ( Supplementary  Fig. 5). To assess the combined effect of heterogeneity in all three differentiation markers, we determined the average cell-to-cell distance in each cell population ( Supplementary Fig. 6). All clonal cell lines exhibited a high degree of heterogeneity in at least one of the differentiation markers and most cell lines expressed substantial plasticity following exposure to Braf/Mek inhibitors. Interestingly, non-transformed melanocytes also exhibited heterogeneity at a level comparable to melanoma cell lines, suggesting that plasticity in differentiation state is not unique to cancer cells.
Because proliferation of BRAF-mutant cells is attributed to their aberrant MAPK signaling, we asked whether heterogeneity in differentiation state was associated with potentially distinct patterns in MAPK signaling. We thus multiplexed single-cell immunofluorescence measurements of p-Erk T202/Y204 , Ki67 (a proliferation marker), and p-S6 S235/S236 (a marker of TORC1 activity that is upregulated in MAPK inhibitor-tolerant cells 36 ) in the same group of cell lines exposed to the same Braf/Mek inhibitors for 3-5 days (Supplementary Figs. 7,8). Fano factor, cellto-cell distance, and t-SNE analysis also revealed substantial variability in MAPK signaling (Fig. 1c,d,Supplementary Figs. 9,10). While vemurafenib led to partial inhibition of MAPK signaling relative to drug-naive cells, the combination of vemurafenib and trametinib suppressed the pathway more strongly and drove a larger proportion of cells toward a fully inhibited (p-Erk Low /p-S6 Low /Ki67 Low ) state (Fig. 1d). The fraction of fully inhibited cells, however, varied among cell lines. We hypothesized that such variability in MAPK signaling would explain differences in the overall MAPK inhibitor sensitivity and might be related to heterogeneity in melanoma differentiation state.
To compare the overall Braf/Mek inhibitor sensitivity among cell lines, we computed drug-induced normalized growth rates (a.k.a. DIP rates 37 ) by normalizing the average net growth rate of each cell line in the presence of each drug to that in DMSOtreated cells (Fig. 1e, f and Supplementary Fig. 11a). By correlating normalized growth rates to the state of MAPK signaling across 16 cell lines, we identified a strong correlation between the fraction of fully inhibited (p-Erk Low /p-S6 Low / Ki67 Low ) cells and the overall sensitivity to Braf/Mek inhibitors (Pearson's r = −0.65, P = 4.3 × 10 -9 ) (Supplementary Fig. 11b and Fig. 1g; top left panel). We next extended the systematic correlation analysis to the diversity of differentiation states. We discovered that upon Braf/Mek inhibition, predominantly undifferentiated (Axl High ) cell lines generated the smallest Fig. 1 Single-cell analysis uncovers heterogeneities in melanoma differentiation, proliferation, and MAPK signaling states across a wide range of Braf/ Mek inhibitor sensitivity. a, b Single-cell protein levels of three melanoma differentiation state markers, Mitf, Ngfr, and Axl, measured by multiplexed immunofluorescence microscopy and overlaid on t-SNE plots. Cells were exposed to either vehicle (DMSO), Braf inhibitor (vemurafenib at 100 nM), or the combination of Braf and Mek inhibitors (vemurafenib at 100 nM and trametinib at 10 nM) for 3 and 5 days. Cells for each experimental condition were randomly selected from the pool of two biologically independent replicates. Single-cell t-SNE maps for all cell lines combined (a) and projections of variation within each individual cell line (b) are shown. c, d Single-cell protein levels of p-Erk T202/Y204 , p-S6 S235/S236 , and Ki67, measured by multiplexed immunofluorescence microscopy and overlaid on t-SNE plots. Treatment conditions are the same as in (a, b). Single-cell t-SNE maps for all cell lines combined (c) and projections of variation within each individual cell line (d) are shown. e Average net growth rates calculated from measurements of live cell count (across at least two biologically independent replicates) at three timepoints (including 0, 3, and 5 days) following exposure to DMSO, vemurafenib (at 100 nM) or vemurafenib (at 100 nM) plus trametinib (at 10 nM). f Drug-induced normalized growth rates (a.k.a. DIP rates) calculated by dividing the average net growth rate for drug-treated cells to that for DMSO-treated cells in each cell line. Normalized growth rates <0 indicate a net cell loss (i.e., drug-induced cytotoxicity), a value of 0 represent no change in viable cell number (i.e., cytostasis), a value >0 indicates a net cell gain, and a value of 1 represents no drug effect as cells grow at the same rate as in the DMSO condition. g Two-sided Pearson's correlation analysis between the fraction of p-Erk Low /p-S6 Low /Ki67 Low cells (referred to as fully inhibited cells) and drug-induced normalized growth rate (top left), the fraction of undifferentiated (Ngfr Low /Axl High ) cells (top right), the fraction of neural crest-like (Ngfr High /Axl Low ) cells (bottom left), and the fraction of melanocytic (Mitf High ) cells (bottom right) across 16 melanoma cell lines treated with Braf/Mek inhibitors for 3 and 5 days. Source data are provided as a Source data file.  Fig. 11b and Fig. 1g; top right panel). Braf/Mek inhibitor resistance in these cells was associated with a high frequency of proliferating (Ki67 High ) cells exhibiting incomplete inhibition of the MAPK pathway. In contrast, populations of neural crest-like (Ngfr High / Axl Low ) cells or differentiated (Mitf High ) cells exhibited substantial co-inhibition of p-Erk, p-S6, and Ki67 (Supplementary Fig. 11b and Fig. 1g; bottom panels). Drug adaptation in these populations was, therefore, associated with an overall reduced requirement for MAPK signaling.
Together, these analyses revealed a spectrum of heterogeneities in melanoma differentiation state that corresponded to two previously described but distinct ways through which cells may tolerate the effect of Braf/Mek inhibition: (1) incomplete inhibition of the MAPK pathway (as seen in Axl High cells), and (2) reduced requirement for MAPK signaling (as seen predominantly in Ngfr High /Axl Low cells).
A chemical screen identifies epigenetic modulators of phenotypic heterogeneity. To systematically search for epigenetic factors that may link heterogeneity in melanoma differentiation state to MAPK dependency, we performed a multi-stage phenotypic screen using a library of 276 epigenetic-modifying compounds. Each compound was used to modulate either a key epigenetic writer, eraser or reader, or a related protein, after which phenotypic responses in the presence or absence of Braf/Mek inhibitors were investigated. For the first stage of the screen, we selected COLO858 and MMACSF cell lines which show distinct patterns of differentiation states (Fig. 2a). Cells from both cell lines were initially exposed to two different doses (0.2 and 1 μM) of each of the 276 epigenetic compounds or vehicle (DMSO) for 24 h. Vemurafenib alone (at 100 nM), vemurafenib in combination with trametinib (at 10 nM), or vehicle (DMSO), was then added and cells were grown for a further 72 or 120 h prior to fixation ( Fig. 2b and Supplementary Tables 1, 2). To differentiate the impact of epigenetic compounds, the growth rates for cells treated with each compound were compared to cells treated without any epigenetic treatment. Statistical analysis identified 58 compounds that led to a significant decrease in normalized growth rate in at least one of the tested conditions (Supplementary Figs. 12,13).
To infer potential variations in the mechanisms of action of the epigenetic compounds, we co-stained cells for p-Erk T202/Y204 , p-Rb S807/S811 , and Mitf, quantifying changes induced by each of the 58 compounds in MAPK signaling, cell cycle progression, and differentiation state, respectively. Unsupervised clustering of the normalized responses revealed remarkable similarities among classes of compounds with common nominal epigenetic targets, suggesting that these responses are most likely the consequence of their on-target effects ( Fig. 2c and Supplementary Fig. 14a). Each class of compounds induced responses that were either cell lineor MAPK inhibitor-specific, or common between both cell lines, or independent of MAPK inhibitor condition. For example, CUDC-907 (Fimepinostat), Quisinostat, Panobinostat, Dacinostat, and Trichostatin A were identified as a cluster of histone deacetylase (HDAC) inhibitors (labeled as group 1 in Fig. 2c) that significantly reduced growth rate in both COLO858 and MMACSF cells to similar degrees and independent of the MAPK inhibitor conditions (Supplementary Fig. 14b). In contrast, Kdm1a inhibitors inhibited net growth rate selectively in COLO858 cells, and they showed a higher efficacy in the absence of Braf/Mek inhibitors.
Correlated patterns of melanoma responses to mechanistically distinct epigenetic inhibitors. To identify potential relationships between the efficacy of each class of epigenetic compounds, we analyzed their effects across a more diverse group of cell lines. Thus, in the second stage of the screen, we focused on seven epigenetic inhibitors representing the most effective classes (including HDAC inhibitors CUDC-907 and Givinostat, pan-Jmj-KDM inhibitor JIB-04, Tankyrase inhibitor AZ6102, BET inhibitors I-BET762 and OTX015, and Kdm1a inhibitor SP2509) and a group of eight cell lines representing a wider spectrum of differentiation states (Fig. 3a, b). By testing each compound in non-transformed human primary melanocytes, we chose effective concentrations of each inhibitor that had little to no effect on healthy cells (Supplementary Fig. 15). To quantify the benefit resulting from combining each of the epigenetic inhibitors with vemurafenib and trametinib, we computed the deviation from Bliss independence (DBI), a metric that compares the observed cellular response to the combination treatment with that expected given independent action for the two individual treatments 38 . We found that the HDAC inhibitor CUDC-907 induced substantial tumor cell killing in all of 8 melanoma cell lines and exhibited additive (independent) to synergistic responses when combined with vemurafenib and trametinib in 6 cell lines (Fig. 3b, c). Treatment with the other 6 epigenetic compounds uncovered more heterogeneous patterns of response across cell lines and MAPK inhibitor conditions. For example, BET inhibitors OTX015 and I-BET762 were not effective in any of the cell lines when used as a single agent, but they induced strong cytotoxic and synergistic responses in combination with vemurafenib and trametinib in the majority of cell lines. SP2509, on the other hand, was effective in 5 cell lines, while exhibiting antagonism in combination with vemurafenib and trametinib in all cell lines.
To systematically explore potential relationships between patterns of responses to epigenetic inhibitors, we computed all pairwise correlations between the efficacy of seven epigenetic compounds across the eight cell lines (Fig. 3d). As expected, the efficacy of two BET inhibitors, OTX015 and I-BET762, was strongly correlated (Pearson's r = 0.99, P < 10 −4 ). In addition, the efficacy of BET inhibitors and the pan-Jmj-KDM inhibitor JIB-04 was positively correlated when cells were co-treated with vemurafenib plus trametinib (r = 0.85, P < 10 −3 ). In contrast, a negative correlation was observed between the efficacies of JIB-04 and SP2509. We thus asked whether such correlations would be maintained if these compounds were assayed across a larger panel of cell lines. By expanding the analysis to a total of 16 cell lines, we confirmed the statistical significance of the negative correlation between responses to SP2509 and JIB-04 (Fig. 4a). In more than half of the cell lines tested, 5 days of treatment with SP2509 induced substantial cell killing ( Fig. 4b and Supplementary  Fig. 16a). SP2509 sensitivity correlated with sensitivity to SP2577 (seclidemstat; a clinical formulation of SP2509), which also suppressed melanoma cell growth (when used at a daily dose of 80 mg kg −1 ) in corresponding melanoma xenografts (Supplementary Fig. 16b-d). Interestingly, however, cell lines with the highest level of resistance to SP2509 and SP2577 were sensitive to JIB-04, showing a range of responses from cytostatic to cytotoxic ( Fig. 4a-c, Supplementary Fig. 16a). The triple combination of JIB-04, vemurafenib, and trametinib led to additive (independent) to synergistic cell killing in most JIB-04-sensitive cell lines, whereas the combination of SP2509 with Braf/Mek inhibitors was antagonistic in all SP2509-sensitive cell lines (Fig. 4c). When combined with vemurafenib and trametinib, both BET inhibitors OTX015 and I-BET762 also induced synergistic cell killing in most of the JIB-04-sensitive cell lines, while having minimal effects when used as a single agent. To determine if the epigenetic inhibitor effects persisted with long-term exposure, we extended the duration of growth inhibition assays. We found that the optimal epigenetic treatments identified for three of the MAPK inhibitor-resistant cell lines were highly efficacious, resulting in reduction of up to 1000-fold in live cell count over a period of 20 days (Fig. 4d).
Together, our multi-stage epigenetic screen and systematic correlation analysis identified a set of seemingly distinct epigenetic states whose inhibition, either by the Kdm1a inhibitor SP2509, or by the pan-Jmj-KDM inhibitor JIB-04 or BET inhibitors (when used in combination with Braf/Mek inhibitors), led to tumor cell killing in a melanoma cell line-specific manner. Next, we asked whether such cell line-specific patterns of response may be linked to heterogeneity in MAPK signaling, proliferation, and differentiation states.
Multivariate modeling identifies predictors of epigenetic inhibitor efficacy. To determine whether molecular markers of differentiation state, MAPK activity, or other phenotypic markers could predict the differential efficacy of epigenetic treatments, we measured both baseline levels and treatment-induced variations in Mitf, Ngfr, Axl, Sox10, p-Erk T202/Y204 , p-S6 S235/S236 , and Ki67 across eight melanoma cell lines. Given the possible contribution of epigenetic histone modifications in DNA damage repair 39 , we also included phosphorylated Histone p-H2A.X S139 (a marker of DNA damage response) in our analysis. Multiplexed immunofluorescence measurements revealed how multiple protein markers were up-or downregulated depending on the cell line, MAPK inhibitor condition, and epigenetic treatment (Fig. 5a). We then used partial least square regression (PLSR) analysis 40 to generate models that linked epigenetic treatment-induced changes in growth rates (response variables) to input vectors that combined baseline and treatment-induced changes in signaling and phenotypic data. Models were evaluated using leave-one-out cross-validation ( Supplementary Fig. 17a). Overall, PLSR models built for JIB-04, SP2509, I-BET762, and OTX015 proved remarkably accurate with an average Pearson's correlation coefficient of 0.85 ± 0.06 between the measured and predicted responses ( Fig. 5b and Supplementary Fig. 17b).
To identify those measurements that are predictive of treatment efficacy, we computed the variable importance in the projection (VIP) scores for each PLSR model 41 ( Fig. 5c and Supplementary  Fig. 17c). Among the most important determinants of treatment efficacy were the baseline differentiation state markers Axl and Ngfr as well as measurements of MAPK activity, including p-Erk, p-S6, and Ki67. Informed by the PLSR results, single-cell analysis across the entire panel of cell lines showed that SP2509 was most effective in inhibiting Ngfr Low /Axl High populations of cells ( Fig. 5d; left panels). In contrast, Ngfr High /Axl Low populations were sensitive to the triple combination of JIB-04, vemurafenib, and trametinib ( Fig. 5d; right panels). In addition, BET inhibitors I-BET762 and OTX015 (when combined with Braf/Mek inhibitors) enhanced tumor cell killing most significantly in populations of Mitf Low /Axl Low cells ( Supplementary Fig. 17d). Interestingly, global changes in differentiation state markers induced by 5 days of treatment with neither of the epigenetic compounds were identified as statistically significant by PLSR models. Single-cell analysis of SP2509-sensitive cells (e.g., A375 and A375-NRAS Q61K ), however, revealed partial downregulation of Axl and upregulation of Mitf following treatment with SP2509, which is in agreement with the possibility that Axl High /Mitf Low cells are being selectively eliminated by SP2509 ( Supplementary  Fig. 18a). In contrast, for cell lines that were sensitive to the combination of JIB-04, vemurafenib, and trametinib (e.g., MMACSF and WM115), we observed partial reduction of the frequency of Ngfr High cells when they were treated with the tripledrug combination in comparison with cells treated with vemurafenib and trametinib only ( Supplementary Fig. 18b, c).
Together, these data reveal how correlated patterns of responses to SP2509, JIB-04, and BET inhibitors are linked to the state of MAPK signaling and differentiation state. Cells in undifferentiated and neural crest-like states represent two different forms of MAPK inhibitor tolerance observed at variable frequencies across most melanoma tumors. Their selective sensitivity to the identified epigenetic inhibitors, therefore, supports the hypothesis that there are epigenetic features associated with melanoma differentiation state that may be linked to their state of Braf/MAPK dependency.
Kdm4b and Znf217 protein levels predict differentiation statespecific sensitivity to JIB-04 and SP2509. We sought to utilize our knowledge of interactions between JIB-04 and SP2509 and their protein targets to better understand the origins of their selective efficacy in melanoma cells. JIB-04 is known as a paninhibitor of Jmj-KDMs, suppressing the activities of Kdm4a, Kdm4b, Kdm5a, and Kdm5b, all with IC 50 's < 0.5 μM 42 . We thus asked which of the Jmj-KDM proteins targeted by JIB-04 might explain its inhibitory effect on melanoma cells. We first depleted each of Kdm4a, Kdm4b, Kdm5a, and Kdm5a proteins in two JIB-04-sensitive cell lines (WM115 and WM902B) and one JIB-04resistant cell line (A2058) using four target-specific siRNAs combined into a single pool ( Supplementary Fig. 19). Only depletion of Kdm4b in WM115 and WM902B cells led to a statistically significant decrease in live cell count and enhanced cells' sensitivity to the combination of vemurafenib and trametinib (Fig. 6a). To rule out the possibility of off-targeting by KDM4B siRNAs, we then examined the effects of three constituent KDM4B siRNAs independently. All individual siRNAs reduced Kdm4b expression and cell viability in both WM115 and WM902B, while having no impact on Kdm4a expression in these cell lines ( Supplementary Fig. 20a, b). A2508 cells, on the other hand, were sensitive to KDM1A knockout by CRISPR using three independent sgRNAs ( Fig. 6b and Supplementary Fig. 20c), which correlated with their sensitivity to SP2509.
We then asked whether the selective efficacy of SP2509 and JIB-04 might be explained by the differentiation state-specific expression of their epigenetic targets. We used mass spectrometry-based proteomics data for BRAF-mutant melanoma lines in the Cancer Cell Line Encyclopedia (CCLE) 43 to evaluate possible correlations between the expression of each epigenetic target and the relative expression of Ngfr and Axl. Since SP2509 acts by blocking Kdm1a interaction with its coactivator Znf217 44 , Fig. 2 A chemical screen identifies epigenetic modulators of phenotypic heterogeneity. a t-SNE maps comparing single-cell heterogeneity in differentiation state (Mitf, Ngfr, and Axl) and MAPK signaling (p-Erk T202/Y204 , p-S6 S235/S236 , and Ki67) within two BRAF V600E melanoma cell lines, MMACSF and COLO858, following exposure to vehicle (DMSO), vemurafenib (at 100 nM), alone or in combination with trametinib (at 10 nM), for 72 and 120 h. b Log 2 -normalized changes in live cell count following exposure of COLO858 and MMACSF cells to either DMSO, vemurafenib (at 100 nM), or vemurafenib (at 100 nM) plus trametinib (at 10 nM), for a period of 120 h. Cells were pretreated for 24 h with either DMSO (top panels) or two different doses (0.2 and 1 μM) of each of the 276 epigenetic-modifying compounds (bottom panels). Data for treatments without epigenetic modifiers (top panels) are presented as mean values ± s.d. calculated across n = 276 biologically independent samples examined over 25 independent experiments. Data for treatments with epigenetic modifiers (bottom panels) are presented as the average of n = 2 biologically independent samples. c A schematic representation of the overall procedure of data collection, processing (normalization), integration, and hierarchical clustering using measurements of cellular grow rate, Mitf, p-Rb S807/S811 , and p-Erk T202/Y204 at indicated treatment conditions and timepoints in MMACSF and COLO858 cells. Unsupervised clustering analysis was performed on data collected for 58 epigenetic compounds that led to a statistically significant decrease in normalized growth rate (when used either as a single agent, or in combination with Braf/Mek inhibitors) in either or both cell lines. Prior to clustering, data collected for cells treated with each epigenetic compound and MAPK inhibitor condition (i.e., DMSO, vemurafenib, or vemurafenib plus trametinib) were normalized to cells treated without any epigenetic compound and the same MAPK inhibitor condition. Groups of compounds with similar nominal epigenetic targets are listed on the right side. Source data are provided as a Source data file.
we also included this protein in our analysis. While Kdm1a protein levels did not significantly correlate with melanoma differentiation state ( Supplementary Fig. 21a), Znf217 was among the top 4% proteins that correlated with the expression of Axl relative to Ngfr ( Fig. 6c; left panels). In addition, levels of Kdm4b (but not other JIB-04 targets) were among the top 0.6% proteins whose expression was significantly greater in melanoma cell lines that expressed higher levels of Ngfr relative to Axl ( Fig. 6c; right panels, and Supplementary Fig. 21b). These data suggest that the relative levels of Kdm4b and Znf217 proteins and differentiation markers Ngfr and Axl may predict the selective sensitivity of melanoma cells to JIB-04 and SP2509.
To independently test this hypothesis, we first profiled the levels of all four proteins across our original panel of 16 melanoma cell lines plus an additional group of 5 cell lines by immunofluorescence microscopy (Supplementary Fig. 22a-d).
Following single-cell analysis, we performed correlation analyses similar to those performed using the mass spectrometry data from CCLE. Our analysis confirmed significant pairwise correlations between the fractions of Ngfr Low /Axl High cells and Znf217 High cells on one hand and the fractions of Ngfr High / Axl Low cells and Kdm4b High cells on the other hand ( Fig. 6d; top panels). As expected, the sensitivity of melanoma cell lines to SP2509 and the combination of JIB-04, vemurafenib, and trametinib were also correlated with the fractions of Znf217 High and Kdm4b High cells, respectively ( Fig. 6d; bottom panels). To further validate the predictivity of these protein markers, we divided the panel of 21 cell lines into a test set of 14 cell lines and . Cells were pretreated with the indicated epigenetic compounds for 24 h and then treated for a period of 3-5 days with either DMSO, vemurafenib (at 100 nM), or the combination of vemurafenib (at 100 nM) and trametinib (at 10 nM). c Deviation from Bliss Independence (DBI) values was computed across diverse epigenetic inhibitor treatments in combination with vemurafenib plus trametinib. DBI = 1 represents an independent (additive) effect equal to what is expected for the combination of drugs that act independently, DBI < 1 represents a combined effect stronger than expected for an independent combination (i.e., synergism), and DBI > 1 represents a combined effect weaker than expected for an independent combination (i.e., antagonism). d Two-sided Pearson's correlation between the effects (i.e., normalized growth rates) of mechanistically distinct epigenetic inhibitors (used individually or in combination with vemurafenib or vemurafenib plus trametinib) evaluated across eight BRAF-mutant melanoma cell lines. Source data are provided as a Source data file. a validation set of 7 cell lines. We then used the training set to develop multi-linear regression (MLR) models of treatment efficacy for SP2509 and JIB-04 using protein levels of Ngfr, Axl, Kdm4b, and Znf217. The trained models predicted SP2509 and JIB-04 efficacy in the validation set with high accuracy (Fig. 6e,  Supplementary Fig. 22e). Together, these data demonstrate that the sensitivity of BRAF-mutant melanoma cells to SP2509 and JIB-04 can be predicted based on their differentiation state (Ngfr versus Axl) and the relative levels of Znf217 and Kdm4b proteins.

Discussion
Lineage-specific epigenetic mechanisms and their reprograming following oncogene inhibition can generate drug-tolerant states that diminish the efficacy of cancer-targeted therapies 45 . In BRAF-mutant melanomas, comparable states tolerant of MAPK inhibition have been associated with variations in differentiation state. Through high-throughput profiling of human melanoma cell lines, we uncovered recurrent patterns of differentiation state heterogeneity that are comparable to those previously seen in patient tumors 25 . We found that a melanoma cell's ability to tolerate Braf/Mek inhibitors was associated with its differentiation state. Responses of undifferentiated cells were associated with incomplete inhibition of the MAPK pathway, whereas MAPKinhibited neural crest-like cells adapted to treatment by reducing their requirement for MAPK signaling.
To identify epigenetic features linked to melanoma differentiation states and MAPK dependency, we performed an epigenetic compound screen and identified three classes of compounds that target distinct melanoma cell states associated with either one of the lysine-specific histone demethylases Kdm1a or Kdm4b, or BET proteins. While the survival of non-transformed melanocytes remained unaffected by these compounds, their cytotoxic efficacy in melanoma cells depended on their state of MAPK activity and their baseline (drug-naive) differentiation state. The Kdm1a inhibitor SP2509 was most effective in inhibiting undifferentiated populations of cells that were intrinsically insensitive to Braf/Mek inhibitors, whereas the combination of SP2509 with Braf/Mek inhibitors led to antagonistic interactions. This is consistent with a recent finding regarding the role that Kdm1a may play in disabling BRAF V600E oncogene-induced senescence 34 . Kdm1a inhibition may, therefore, require MAPK signaling to restore senescence in Ngfr Low /Axl High cells. In contrast, Ngfr High /Axl Low populations of cells exhibited additive to synergistic responses to the combination of Kdm4b and Braf/Mek inhibition. BET inhibitors had minimal effect on melanoma cells when used as a single agent but led to tumor cell killing when combined with Braf/Mek inhibitors. Together, our systematic studies extend previous findings about the potential of each of the epigenetic regulators, Kdm1a, Jmj-KDMs, and BET proteins, as therapeutic vulnerabilities in melanomas 20,34,46 in two complementary ways. First, by associating these vulnerabilities to the patterns of single-cell heterogeneity in differentiation state, we explain how they may vary from one tumor to another. Second, by linking epigenetic vulnerabilities to the state of Braf/MAPK dependency, we determine whether each epigenetic inhibitor reaches its maximal efficacy when used as a single agent or when combined with Braf/Mek inhibitors. Cell-to-cell heterogeneity in the state of oncogene dependency poses a general challenge to the use of cancer-targeted therapies. Our systems pharmacology approach provides a promising avenue toward the identification of actionable epigenetic factors that may extend the oncogene addiction paradigm on the basis of tumor cell differentiation state. We demonstrate the utility of this approach in BRAF-mutant melanomas by using JIB-04 and SP2509, as investigational tools, to block epigenetically diverse populations of MAPK inhibitor-tolerant cells. We show that the relative baseline levels of Kdm4b and Znf217 (a Kdm1a coactivator), determine differentiation state-specific sensitivity of melanoma cells to their corresponding inhibitors. In cell lines that were highly sensitive to SP2509, 3 to 5 days of treatment with SP2509 reduced the fraction of Axl High /Mitf Low cells in the tumor cell population. In contrast, in cell lines that were sensitive to the combination of JIB-04 and Braf/Mek inhibitors, the frequency of Ngfr High cells was only partially reduced following treatment with the triple-drug combination. A possible explanation for not observing drastic changes in differentiation state markers after epigenetic inhibitor treatments may be that the measured Ngfr and Axl proteins are surrogate markers which, although confirmed to be statistically associated with melanoma differentiation states, are not causally linked to melanoma cell survival 12,19 . Therefore, additional studies are required to identify key downstream targets of histone demethylases and their interactions with melanoma differentiation state that make them selectively sensitive to each epigenetic inhibitor. This may be accomplished via a combination of transcriptomic, proteomic, and epigenomic profiling of melanoma cells of diverse differentiation states following treatment with each inhibitor in the absence or presence of MAPK inhibitors. The identification of Kdm4b and Znf217/ Kdm1a, whose knockdown phenocopies the epigenetic inhibitors' effects, will facilitate such mechanistic studies, e.g., through rationally designed chromatin immunoprecipitation coupled with sequencing (ChIP-seq) experiments.
Our data add to a growing body of research 18,[47][48][49][50] showing that even genetically homogenous populations of tumor cells consist of subpopulations at diverse differentiation states that are also different with respect to their state of oncogene dependency. Stochastic fluctuations in the levels of key proteins from pathways that have crosstalk with the oncogenic pathway and differentiation state-specific epigenetic modifications represent possibly related mechanisms that lead to diverse states of drug tolerance 2,51,52 . A quantitative understanding of the key drivers of such phenotypically consequential epigenetic states at a single-cell level and uncovering how they are linked to signaling networks that function in feedback regulation of oncogenic signaling, are likely key steps to improving the effectiveness and durability of response to oncogene-targeted therapies.
Small-molecule inhibitors, including a library of 276 epigenetic-modifying compounds, chemicals used in the follow-up cell-based assays (Givinostat, CUDC-907, JIB-04, AZ6102, I-BET-762, OTX015, and SP2509), as well as vemurafenib and trametinib were all purchased from Selleck Chemicals. SP2577 was purchased from MedChem Express. The complete list of compounds used in this study, their catalog numbers and purity, as evaluated by HPLC and MS analysis, are presented in Supplementary Table 1. Compounds used for cell-based studies were dissolved in the appropriate vehicle (either DMSO or water) at a stock concentration of 10 mM. The vehicle for SP2577 used for in vivo studies is described below.
The following primary monoclonal antibodies (mAb, clone name) and polyclonal antibodies (pAb) with specified animal sources, catalog numbers, research resource identifiers (RRID), and dilution ratios, were used in immunofluorescence staining assays: Mitf (mouse mAb, clone D5, Abcam, Cat# ab3201, AB_303601, 1:800), p-Erk T202/Y204 (rabbit mAb, clone D13.14.4E, Cell Signaling Technology, Cat# 4370, AB_2315112, 1:800), Ki67 (mouse mAb, clone 8D5, Cell Signaling Technology, Cat# 9449, AB_2715512, 1:1200), Axl (goat pAb, R&D Systems, Cat# AF154, AB_354852, 1:400), p-Rb S807/S811 (goat pAb, Santa Cruz Biotechnology, Cat# sc-16670, AB_655250, 1:400), Ngfr (rabbit mAb, clone  (10 nM). Protein data shown for each condition represent mean values across three timepoints and two biologically independent samples and are then z-scored across all cell lines and treatment conditions. b Two-sided Pearson's correlation between responses (normalized growth rates) to SP2509 (left) or JIB-04 in combination with vemurafenib and trametinib (right), measured for each of the eight melanoma cell lines (x-axis) and corresponding responses predicted by partial least square regression (PLSR) modeling following leaveone-out cross-validation (y-axis). c PLSR-derived variable importance in the projection (VIP) scores, highlighting combinations of protein measurements at the baseline (shown in red), following epigenetic inhibitor treatment (shown in blue), or the ratio of change induced by each epigenetic compound (shown in black), that are predictive of efficacy for SP2509 (left) and JIB-04 in combination with vemurafenib and trametinib (right). The sign of VIP score shows whether the change in variable correlated negatively or positively with the treatment-induced response. Only VIP scores of greater than 1 or smaller than −1 with a statistically significant Pearson's correlation (P < 0.05) are highlighted. d Two-sided Pearson's correlation analysis between the baseline fractions of Ngfr Low /Axl High cells and measurements of normalized growth rate in response to SP2509 when used as a single agent (left), and between the baseline fractions of Ngfr High /Axl Low cells and response to JIB-04 in combination with vemurafenib and trametinib (right), across 16 BRAF-mutant melanoma cell lines. The significance of differences between normalized growth rates was also evaluated based on two-sided Mann-Whitney U-test. For this test, 16 cell lines were divided into two groups of eight based on whether the measured variable had a value above or below the median. Data for each group were then presented using box-and-whisker plots, on which the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to 1.5× the interquartile range as a measure of variance, and datapoints outside the range are plotted individually with asterisks. Source data are provided as a Source data file. Immunofluorescence staining, quantitation, and analysis. Cells in 96-well plates were fixed in either 4% paraformaldehyde (PFA) for 20 min at room temperature or 100% ice-cold methanol for 15 min at −20°C. Cells were then washed with PBS, permeabilized in methanol (in cases where they were fixed with 4% PFA) for 10 min at −20°C, rewashed with PBS, and blocked using Odyssey blocking buffer (LI-COR Biosciences) for 1 h at room temperature. Cells were incubated overnight (~16 h) at 4°C with primary antibodies in Odyssey blocking buffer. The following day, cells were washed three times with PBS supplemented with 0.1% Tween-20 (Sigma-Aldrich) (PBS-T) and incubated for 1 h at room temperature with the secondary rabbit, goat, or mouse antibodies. Cells were then washed twice with PBS-T, once in PBS, and then incubated with Hoechst 33342 (Thermo Fisher, Cat# H3570, 1:20,000) for 20 min at room temperature. Cells were washed twice with PBS and imaged with a ×10 objective using the ImageXpress Micro Confocal High-Content Imaging System (Molecular Devices) or the Operetta CLS High-Content Imaging System (Perkin Elmer). A total of nine sites were imaged per well. Background subtraction was performed with ImageJ. Image segmentation and quantification of signal intensities in the whole cell, nucleus, cytoplasm, or the nucleus/cytoplasm (N/C) ratio were performed with CellProfiler 53 . Population- average and single-cell data were analyzed using MATLAB 2019b. By generating histograms of single-cell data across a variety of conditions for each protein (X), including positive and negative controls, we identified an appropriate binary gate, based on which the percentage of X High versus X Low cells in each condition was quantified.
Single-cell analysis and dimensionality reduction. To visualize single-cell heterogeneities in a two-dimensional space, we used a population of 6069 cells assembled following random selection of 60 individual cells (or fewer in the case of the highly drug-sensitive C32 cell line, for which <60 cells survived following treatments with the combination of vemurafenib and trametinib) from each of the 102 tested conditions, covering the entire panel of 17 cell lines, 3 drugs (DMSO, vemurafenib alone, and vemurafenib plus trametinib) and 2 timepoints. Logtransformed single-cell data for multiplexed measurements of differentiation state markers, Mitf, Ngfr and Axl, and for multiplexed measurements of MAPK signaling, p-Erk T202/Y204 , p-S6 S235/S236 , and Ki67, were processed by z-scoring (for each protein measurement) across all 6069 cells, followed by principal component analysis (PCA) using the MATLAB built-in function pca. Eighty percent of the variance in the differentiation state single-cell data and 92.5% of the variance in MAPK single-cell data were captured by the first two principal components (PCs). We thus performed t-distributed stochastic neighbor embedding (t-SNE) analysis on the scores from the first two PCs of each of the single-cell datasets by applying the built-in tsne function in MATLAB 2019b. We used the barneshut algorithm, learning rate of 1000 for the optimization process, a maximum number of optimization iterations of 2000, perplexity of 480, and exaggeration factor of 4. To analyze the baseline state and treatment-induced changes in the single-cell behavior of each cell line, we used the t-SNE map overlaid with the projection of z-scored single-cell measurements (color-coded between blue and red) for that cell line while showing data for other cell lines in gray.
To compare baseline and treatment-induced heterogeneities in the expression of each protein, we evaluated Fano factor (a.k.a. index of dispersion) by computing the variance-to-mean ratio for single-cell measurements across each population of cells. To compare single-cell heterogeneity in multiple differentiation state markers (Mitf, Ngfr, and Axl) simultaneously, we determined the average cell-to-cell distance by computing the pairwise Euclidean distance among individual cells (after z-scoring the data for each protein marker), and the mean of all possible pairwise distances was reported for each population of cells.
Measurements of growth rate, drug sensitivity, and combination effectiveness. Growth rate inhibition assays were performed in 96-well clear bottom black plates (Corning, Cat# 3904). Cells were counted using a TC20 automated cell counter (Bio-Rad Laboratories) and seeded in 200 µl of full growth media at a density of 1200-5000 cells per well, depending on the baseline proliferation rate of each cell line. Using a D300e Digital Dispenser, cells were treated the next day with smallmolecule compounds at reported doses or vehicle (DMSO). Measurements of live cell count across multiple timepoints were then used to calculate the net growth rate for each treatment condition. To measure the number of surviving cells at each timepoint, cells were fixed in 4% paraformaldehyde for 20 min at room temperature. Cells were washed twice with PBS and incubated with Hoechst 33342 (Thermo Fisher, Cat# H3570, 1:20,000) for 20 min at room temperature. Cells were washed again with PBS and imaged with a ×10 objective using an ImageXpress Micro Confocal High-Content Imaging System (Molecular Devices) or an Operetta CLS High-Content Imaging System (Perkin Elmer). A total of nine sites were imaged per well. Nuclear segmentation and cell counting was performed using CellProfiler 53 .
The net growth rate for cells treated with individual compounds (including MAPK inhibitors or epigenetic inhibitors), or their combination, or vehicle, was calculated from time-dependent changes in live cell count according to the following equation: where N t1 and N t2 represent the number of cells measured at timepoints t = t 1 and t = t 2 , respectively, and μ describes the net growth rate of cells during the time period between t 1 and t 2 . Average net growth rates for each treatment condition were calculated as the mean of growth rates measured across multiple consecutive timepoints. To compare drug sensitivity among cell lines while correcting for differences in their baseline proliferation rates, drug-induced normalized growth rates (a.k.a. DIP rates 37 ) were computed as follows: As shown above, normalized growth rate for cells treated with a specific drug is calculated by normalizing the average net growth rate measured for drug-treated cells (μ drug ) to that measured for vehicle (DMSO)-treated cells (μ DMSO ). Normalized growth rates <0 indicate a net cell loss (i.e., cytotoxicity), a value of 0 represents no change in viable cell number (i.e., cytostasis), a value >0 indicates a net cell gain (i.e., population growth), and a value of 1 represents cell growth at the same rate as in vehicle-treated (control) cells.
To quantify the benefit resulting from combining two drugs, e.g., an epigenetic inhibitor (A) with a MAPK inhibitor (B), we evaluated their interactions based on the Bliss independence model 38 . We first calculated the fraction of cells affected by each treatment as follows: where min(DIP) represents the minimum of DIP values reported across all cell lines and drug treatment conditions in this study. Using f a , we overcome the Bliss metric limitation for the analysis of unbounded drug effects such as the normalized growth rates 54 , while highlighting combined interactions that influence drug efficacy, a parameter that is affected (more obviously than potency) by cell-to-cell heterogeneity and the presence of small subpopulations of drug-tolerant cells 55 . We then used measurements of f a for individual and combination treatments to compute the deviation from Bliss Independence (DBI) as follows: where f a (A) and f a (B) represent the fraction of cells affected by each drug separately, and f a (A + B) represents that of the combination treatment (A plus B). Based on this definition, the calculated DBI compares the observed response to that expected given independent action for the two individual treatments. DBI = 1 represents an independent (additive) effect equal to what is expected for the combination of drugs that act independently, DBI < 1 represents a combined effect stronger than expected for an independent combination (i.e., synergism), and DBI > 1 represents a combined effect weaker than expected for an independent combination (i.e., antagonism). Epigenetic compound screen. COLO858 and MMACSF cells were counted using a TC20 automated cell counter (Bio-Rad Laboratories) and seeded in 96-well clear bottom black plates at a density of 2000 and 3000 cells per well (excluding 36 wells at the edges), respectively. Using the HP D300e Digital Dispenser, cells were treated the next day with either vehicle (DMSO or water) or two different doses (0.2 and 1 µM) of each of the 276 compounds in an epigenetic compound library. After 24 h, cells were either fixed or treated with vemurafenib alone (at 100 nM), or vemurafenib (at 100 nM) plus trametinib (at 10 nM), or vehicle (DMSO). Cells were then grown for a further 72 or 120 h prior to fixation. All experimental conditions with an epigenetic compound treatment (at each dose and timepoint) were tested in two replicates. Experimental conditions that did not include an epigenetic compound treatment (i.e., treatments with only vehicle, vemurafenib, or vemurafenib in combination with trametinib) were repeated in all 96-well plates during the entire period of compound screening, creating a total of 276 replicates that were used to evaluate plate-to-plate and day-to-day experimental robustness. To differentiate the impact of epigenetic compounds on each cell line at different states of MAPK signaling, the growth rates for cells treated with each dose of an epigenetic compound were compared with growth rates for cells treated without any epigenetic compound, following 72 h and 120 h of exposure to each of the MAPK inhibitor treatment conditions (i.e., DMSO, vemurafenib, and vemurafenib plus trametinib). Significant epigenetic compounds were defined as those that led to a statistically significant decrease in normalized growth rate (with an effect size of at least 0.5 day −1 ) in at least one of the tested conditions (i.e., cell lines, timepoints, or MAPK inhibitor conditions). Statistical significance was evaluated as P < 0.05 based on two-sample t-test following correction for multiple comparisons using the Dunn-Sidak method. Additional details about the epigenetic compound screen are available in Supplementary Table 2.
Hierarchical clustering. To infer potential variations in the mechanisms of action of the selected group of 58 epigenetic compounds, treatment-induced changes in growth rate, p-Erk T202/Y204 , p-Rb S807/S811 , and Mitf across different cell lines and MAPK inhibitor conditions were integrated into a matrix for unsupervised clustering. Growth rates for epigenetic inhibitor conditions were averaged across three timepoints and two doses of the inhibitor and their differences relative to the timeaveraged growth rate in the absence of any epigenetic compound were used for each MAPK inhibitor treatment condition (i.e., DMSO, vemurafenib, and vemurafenib plus trametinib). For p-Erk T202/Y204 , p-Rb S807/S811 , and Mitf, logtransformed signal intensities (from quantitative immunofluorescence microscopy) were averaged across three timepoints and two doses of the inhibitor and their differences relative to the corresponding MAPK inhibitor condition in the absence of any epigenetic compound were used. Unsupervised hierarchical clustering of the integrated data was then carried out using MATLAB 2018b with the Correlation distance metric and the Complete (farthest distance) algorithm for computing the distance between clusters.
Partial least square regression (PLSR) analysis. We used PLSR analysis 17,40,56 to generate models that linked epigenetic treatment-induced changes in normalized growth rates to baseline and treatment-induced changes in signaling and phenotypic data. By combining data for each epigenetic inhibitor across eight cell lines, we generated one model to predict responses to each epigenetic inhibitor, when used either individually or in combination with Braf/Mek inhibitors. Response variables for each model were defined as normalized growth rates (averaged across 3-5 days of treatment) in the absence or presence of that compound (together with either DMSO, vemurafenib alone, or vemurafenib in combination with trametinib) across the eight cell lines. Input vectors were constructed by combining measurements of a total of eight proteins (including signaling proteins, differentiation state and phenotypic markers) at the baseline (i.e., either DMSO, vemurafenib alone, or vemurafenib in combination with trametinib), following epigenetic inhibitor treatment (together with either DMSO, vemurafenib alone, or vemurafenib in combination with trametinib), and their changes induced by the epigenetic inhibitor which were computed by taking the ratio of the two sets of measurements. All protein measurements were performed by quantitative immunofluorescence microscopy, log-transformed, and averaged cross three timepoints (including 24 h after epigenetic inhibitor or DMSO treatment, followed by 72 and 120 h of MAPK inhibitor or DMSO treatment). The data were then z-scored across all conditions and cell lines prior to the application of PLSR analysis using the built-in MATLAB function plsregress.
To evaluate the predictability of the linear relationship between the input and output variables in each model, we used leave-one-out cross-validation. The goodness of fit for each model was calculated using R 2 . Prediction accuracy was evaluated by Q 2 and the P values generated from pairwise Pearson's correlations between the measured and predicted responses following cross-validation. For the assessment of relative variable importance in each PLSR model, the information content of each variable was assessed by its variable importance in the projection (VIP) 41 .
Multi-linear regression (MLR) analysis. To statistically test the hypothesis that the relative levels of Kdm4b, Znf217, and differentiation state markers Ngfr and Axl can predict the selective sensitivity of melanoma cells to JIB-04 and SP2509, we performed MLR analysis using the entire panel of 21 cell lines, which were divided into a test set of 14 (to tarin models) and a validation set of 7 cell lines (to evaluate the model performance). To group the cell lines into training and validation sets, we first sorted the 21 cell lines based on their sensitivity to each compound. Beginning from the most sensitive or resistant cell line, we picked 2 cell lines for the training set and 1 cell line for the validation set, then another 2 cell lines for the training set and 1 cell line for the validation set, etc. We then used the training set to develop MLR models of treatment efficacy for SP2509 and JIB-04 using protein levels of Ngfr, Axl, Kdm4b, and Znf217 measured by quantitative immunofluorescence microscopy. Response variables were defined as normalized growth rates (averaged across 3-5 days of treatment) in the presence of either SP2509 (at 1 µM) or the combination of JIB-04 (at 0.2 µM), vemurafenib (at 100 nM), and trametinib (at 10 nM). Input vectors were constructed by combining baseline (drug-naive) measurements of the fractions of Ngfr High , Axl High , Kdm4b High , and Znf217 High cells in each cell line. We then used the built-in MATLAB function regress to return a vector of coefficient estimates for an MLR of the responses to each treatment condition. These coefficients were then used to predict the responses of the validation cell lines to each treatment using an input vector generated for those cell lines. To evaluate the predictability of the models, we evaluated the relationship between the measured and predicted responses for the validation set of cell lines using Pearson's correlation analysis. The goodness of fit for each model was calculated using R 2 of measured versus fitted responses for models generated using all 21 cell lines. Gene knockout by CRISPR-Cas9. A2058 cells were seeded in a 24-well plate at a seeding density of 30,000 cells per well for 24 h. The next day, cells were treated at MOI of 0.3 with the pre-designed Edit-R lentiviral Cas9 nuclease (Dharmacon, Cat# VCAS10124). Polybrene (at 8 μg/ml) was added to enhance the efficiency of the viral infection. Following an incubation period of 72 h with the lentivirus, Cas9-lentiviraltreated cells were selected with blasticidin over a period of 1.5 weeks. Surviving A2058 Cas9-positive cells were subsequently infected with either non-targeting control (Dharmacon, Cat# VSGC10216) or three independent KDM1A lentiviral sgRNAs (Dharmacon, Cat# VSGH10142-246522030, VSGH10142-246999191, and VSGH10142-246522028). The efficiency of KDM1A knockout and the consequential changes in cell viability were measured via immunostaining for the Kdm1a protein at 96 h, followed by quantitative immunofluorescence microscopy.
siRNA and sgRNA sequences. All siRNA and sgRNA sequences can be found in Supplementary Table 3.
In vivo xenograft assays. All mouse experiments were carried out in accordance with procedures approved by the Institutional Animal Care and Use Committee (IACUC) at the University of Michigan. Athymic, 5-6 weeks old female nude (NU/J) mice were purchased from The Jackson Laboratory. Mice were housed under standard animal housing protocol at the University of Michigan Unit for Laboratory Animal Care (ULAM) in special containment rooms in individual nonventilated cages at an ambient temperature of 74°F with a humidity of 30-70% and the light-dark cycle of 12 h each. For xenograft tumor injections, mice were first anesthetized using 3% vaporized isoflurane. 2.5 × 10 6 cells from either of the melanoma cell lines A375 or WM2664 suspended in 200 μl of growth factorreduced Matrigel (Thermo Fisher, CB-40230C) in PBS (1:1) were injected subcutaneously in the right flank of each mouse. Tumor xenografts were monitored three times a week with digital calipers. Once the tumors were palpable and the mean volume across all tumors reached a volume of~170 mm 3 , the mice were randomly allocated to two treatment groups of 5 mice per group. SP2577 (at a daily dose of 80 mg kg −1 ) or vehicle (20% DMSO, 20% Cremophor EL, plus 60% sterile water) was then administered to each mouse via I.P. injection of a solution of 200 μl. Mouse body weight, tumor volume, and general health were monitored three times a week for 12 days. At the endpoint of the study, mice were euthanized by CO 2 asphyxiation and bilateral pneumothorax.
Statistics and reproducibility. No statistical method was used to predetermine sample sizes. Sample sizes were chosen based on similar studies in the relevant literature. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment. All data with error bars were presented as mean values ± s.d. or mean values ± s.e.m. as indicated in figure legends using indicated numbers of biologically independent replicates. The significance of pairwise correlations among drug response data were evaluated based on P values associated with the corresponding two-sided Pearson's correlation analysis (r = correlation coefficient). The statistical significance of the effects of epigenetic compounds (in the first stage of the screen) was evaluated based on P < 0.05 generated from two-sided t-test following correction for multiple comparisons using the Dunn-Sidak method. To identify the statistical significance of differences between mean or median of measurements within two different groups, P values from two-sided t test or two-sided Mann-Whitney U-test were used, respectively. Statistical significance of the difference in % change in tumor size in mouse experiments was determined by two-way analysis of variance (ANOVA). Statistical analyses were performed using MATLAB 2018b and 2019b.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All data generated in this study are included in this published article and its supplementary information files. Because of the large number and size of raw immunofluorescence microscopy images associated with this study, all relevant image files will be made available by hard drive upon reasonable request to the corresponding author. The Cancer Cell Line Encyclopedia (CCLE) mass spectrometry-based proteomics data analyzed in this study were downloaded from the depmap project portal (https://depmap. org/portal/download/). Source data are provided with this paper.