Superenhancer Analysis Defines Novel Epigenomic Subtypes of Non-APL AML, Including an RARα Dependency Targetable by SY-1425, a Potent and Selective RARα Agonist

Acute myeloid leukemia (AML) is the most common acute leukemia in adults (1), with an estimated 4 new cases diagnosed per 100,000 individuals in the United States in 2015 (SEER data, 2016, ref. 2). Current treatment involves chemotherapy-based induction and consolidation with potential curative use of hematopoietic cell transplantation (HCT) in those eligible patients who respond to treatment. New therapeutic options are needed for patients with AML given a 5-year overall survival of only 35% to 40%. Myelodysplastic syndrome (MDS) is a heterogeneous hematopoietic disease commonly associated with bone marrow failure. Although early MDS is considered an indolent disease, the prognosis of patients with later-stage MDS is poor and is associated with a low overall survival similar to AML. We sought to characterize the epigenomic circuitry and transcriptional state of leukemic cells to identify potential new treatments beyond those discovered by mutational analysis of protein coding genes.

To date, a large part of cancer research has focused on somatic mutations in protein coding regions to identify putative oncogenic drivers (3). Recently, large, highly active chromatin regions, named “superenhancers” (SE), have been described that define cell identity and cell state in normal and malignant cells, including cancer cells (4–6). Although typical enhancers are short (∼1 kb) stretches of DNA that regulate the expression of genes through the binding of transcription factors (TF) and subsequent recruitment of RNA Pol II and associated proteins, SEs are larger in size (>20 kb on average) and are characterized by a much higher occupancy of transcriptional regulatory proteins (5, 6). Previous work has shown that the histone three lysine-27 acetylation (H3K27ac) mark can be used as an efficient and robust means of SE demarcation (6–13). It has been shown that within a given cell, these few, highly active chromatin regions (<5% of all enhancers) tend to regulate key genes governing cell identity and phenotype. In tumor cells, this includes oncogene and nononcogene drivers of the transformed state (7, 8).

To define contributions of the noncoding genome in AML, we characterized the enhancer, SE, gene expression, and mutational landscapes of blasts from 66 patients with AML and used this information to dissect transcriptional regulatory circuitry and predict therapeutic targets. De novo stratification of these enhancer profiles identified 6 distinct patient subgroups, independent of and orthogonal to mutational status. From this, we demonstrated that retinoic acid receptor alpha (RARα) occupies a critical node in the transcriptional network of a new subgroup of AML, putatively identifying a novel tumor cell liability that could be exploited therapeutically.

The selective and potent RARα agonist SY-1425 (tamibarotene) has been demonstrated to possess improved pharmacologic properties over first-generation pan-retinoids such as all-trans retinoic acid (ATRA; refs. 14–18). We found that SY-1425 potently and selectively inhibits the proliferation of AML cell lines that contain a strong SE at the RARA locus. This observation was confirmed in AML patient-derived xenograft (PDX) models whereby RARA-high AML models were more sensitive to SY-1425 treatment than RARA-low AML models. We also observed RARA SEs and overexpression of RARA mRNA in MDS. These data provide the foundational support for investigating SY-1425 as a treatment option for patients with AML or MDS with a strong RARA SE.

Using a high-throughput H3K27ac chromatin immunoprecipitation (ChIP) platform (Fig. 1A), we profiled the landscape of active enhancers in blasts from 66 AML donors, 28 AML cell lines, 4 FACS-purified hematopoietic stem and progenitor cell (HSPC) samples, and 6 FACS-purified monocyte samples (Supplementary Fig. S1A). Four patient samples were excluded from the analysis due to low ChIP-seq quality. Using methods described previously (7), we identified individual enhancer loci and SEs, revealing a median of 807 SEs per sample with a median length of 22 kb (Fig. 1B and Supplementary Fig. S1B and S1C) with enhancer “score” reflecting both enhancer size and density of reads (see Supplementary Methods). As identified in previous studies, SEs are large, highly active regions of chromatin that regulate a small subset of genes critical to cell identity and state (4–6). This approach provided robust enhancer profiles with high technical replicate consistency (Supplementary Fig. S1D) and enabled the identification of a sample purity threshold (% AML blast cells) of 80% as a requirement for patient inclusion in this study (Supplementary Fig. S1E–S1G).

To link SEs to genes, we utilized the correlation between SE score and H3K27ac abundance at promoters of genes within 500 kb of the enhancer, predicting a link when the P value of the correlation was significant genome-wide, or when an SE lacked any significant correlation, linking it to the best correlated and nominally significant gene (Supplementary Fig. S2A–S2D; see Supplementary Methods). Most SEs are linked to multiple genes with this approach, yielding on average 2.7 links per SE (Supplementary Fig. S2E). We then compared the SE to gene linking from this correlation method in patients with primary AML to published chromatin conformation data from promoter capture Hi-C (PCHi-C) in monocytes (19) and found that 62% of SEs had validated linking.

Previous work has identified the differentiation state of AML cells as a primary source of epigenetic variation (20–22). This is consistent with the observation that certain key HSPC enhancers, such as the SE linked to CDK6, show strong variability across our AML cohort (Fig. 1C). In an effort to assess how well SE maps can distinguish differentiation states in primary AML samples, we defined the extremes of myeloid differentiation by constructing SE maps from HSPCs and monocytes. We defined signatures of differentiation state based on the top 20 HSPC- and monocyte-specific SEs (Fig. 1D). Assessment of this enhancer-based HSPC signature within the 62 primary AML blasts showed a robust inverse correlation between HSPC-specific enhancers and monocyte-specific enhancers (Supplementary Fig. S3A). To quantify AML blast utilization of an HSPC-like enhancer signature, we derived a predicted HSPC score (Fig. 1D, green bar). This score defines the stem-like (high HSPC score) or differentiated (low HSPC score) state of a given AML sample. To determine if this differentiation state is linked to one of the major genome-wide enhancer programs, we used independent components analysis (ICA) to find the two most pervasive enhancer signatures in the genome (Fig. 1E). The first IC, termed myeloid differentiation, accounted for, on average, 21.4% of the variance of SEs genome-wide and correlated strongly with the predicted HSPC score (R² = 0.85; Supplementary Fig. S3B). The second IC, termed HOX factor activation, accounted for an average of 11.5% of the variance and correlated with enhancers associated with homeobox (HOX) genes, PBX3 and MEIS1, among others (Supplementary Fig. S3C). These TFs have been shown to co-occupy enhancers in murine models of leukemia, acting in concert to drive leukemogenesis and enforce immature differentiation states (23, 24). Together, these two independent components accounted for 33% of the variance in the data, providing motivation for a higher-dimensional analysis.

To this end, we performed nonnegative matrix factorization consensus clustering (NMF-C; see Supplementary Methods) among the patient samples to group the patients by enhancer score. This approach identified 6 distinct subgroups of patients (referred to as clusters) in our AML cohort (Fig. 1F), each characterized by a unique combination of enhancer profiles (Supplementary Fig. S3D and S3E). Although not included in the factorization, the HSPC and monocyte H3K27ac profiles were inserted into the clustering to illustrate their respective enhancer profiles (Fig. 1F, light blue and purple, respectively). This demonstrates that cluster 2 was characterized by enhancers also found in monocytes, whereas cluster 5 utilized enhancers active in HSPCs. Additionally, the 6 clusters showed robust internal association by Pearson correlation of enhancer activity (Supplementary Fig. S3F). Despite clear separation in NMF-C, some clusters also showed strong intercluster correlation, with 2 groups of clusters—1 and 2 and 3–6—showing high intercluster similarity (Supplementary Fig. S3F). This similarity is driven by the HSPC signature score of the samples within the clusters, with clusters 1 and 2 being far more monocyte-like than clusters 3–6 (Fig. 1G; Supplementary Fig. S3G).

Using the clusters defined by overall SE distribution, we then focused on identifying the specific SE-linked genes forming the basis of these clusters. Prominent among these were SEs associated with known AML master TF genes (e.g., HOXA/HOXB/MEIS1), an HSPC-specific gene (PROM1), and a monocyte-specific factor gene (IRF8; Fig. 2A). These core TFs suggest that the clusters arise from shared underlying master regulatory circuitry.

To assess the interplay between somatic mutations in protein coding genes and our SE-based clustering, we performed targeted genotyping on all samples to identify mutations in known leukemia-associated genes (Fig. 2B; Supplementary Table S1). The 130 genes sequenced included any gene found by whole- exome sequencing to be mutated in more than 2% of AML cases (ref. 25; Supplementary Table S1). These patient samples were also profiled for the presence of common prognostic cytogenetic abnormalities as part of routine clinical pathology tests. Comparing targeted genotyping and cytogenetics with our clustering, we observe that all patients with known mixed-lineage leukemia (MLL) translocation are sorted into cluster 1 by their global enhancer profiles (Fig. 2B, P < 1.6E−7). However, we note that 12 of our patients have not been genotyped for MLL translocation. Nevertheless, this suggests that MLL translocations may induce a unique epigenetic state that dominates the overall enhancer landscape of the given AML sample, not surprising given the function of MLL fusion proteins in misregulating DOT1L function, leading to aberrant histone H3K79 methylation (26). Although mutations in NPM1 and FLT3 are enriched in clusters 2, 3, 4, and 6, no cluster-specific associations were observed for mutations in DNMT3A, TET2, IDH1, and IDH2.

To better understand the cluster-specific SEs driving our patient grouping, we visualized the median per-cluster SE score for the most influential SEs in the NMF decomposition (Fig. 2C). Clear cluster-specific SEs were dominated by genes with known roles in AML such as MEIS1, PROM1, SOCS2, and ERG, among others (27, 28). A pattern of cluster-specific SE usage emerged, with each cluster being regulated by a relatively compact number of individual SEs.

To understand whether clustering of AML patient samples by enhancer landscape is biologically meaningful, we asked whether SE-defined clusters could be associated with other biological differences by bridging to a large patient data set (25). Our patient sample cohort showed a bias toward high blast count and FLT3-ITD positivity compared with other surveys of AML (ref. 29; Supplementary Fig. S4A) and showed a median patient survival of 13.1 months from diagnosis (Supplementary Fig. S4B). However, the relatively small numbers of patients in the SE-based clustering make robust correlations with clinical factors difficult. As an exploratory investigation, we attempted to bridge to well-annotated expression-based databases because ours was already the largest SE database in AML available. Thus, we developed a 1,710-gene transcriptional signature as a proxy for our SE profiling clusters and projected this scheme into the patients profiled by The Cancer Genome Atlas (TCGA) Consortium (Supplementary Fig. S4C and S4D; Supplementary Table S2; see Methods). By expression, the TCGA clusters 1 and 2, which already showed some overlapping SEs, were combined in the classification because their transcriptional signatures were too similar to reliably distinguish their members on cross-validation. Although this analysis is not meant to replace established clinical AML prognostic categories or patient stratification, we note that enhancer clusters were associated with significantly different survival outcomes in univariate analysis (P < 0.013; Supplementary Fig. S4E; Supplementary Table S3), which could not be accounted for by somatic mutation differences (Supplementary Fig. S4F). This suggests that enhancer profiles may provide orthogonal information to mutations in the coding region of the genome, motivating future studies into the clinical relevance of enhancer landscapes.

To identify targets for putative therapeutic intervention, we inferred a mutual information network of coordinated SEs across the clusters, and hypothesized that TFs in this network are key drivers of the transcriptional circuitry in these patients (Fig. 3A). Some TFs showed concomitant enrichment of their DNA binding motif in cluster-specific enhancers, setting the foundation for the cell's nuclear circuitry (Fig. 3A, red outline). One enhancer that was highly differential between patient samples is associated with the retinoic acid receptor α (RARA) gene (Fig. 3B and C; Supplementary Fig. S5A). In some cases, this enhancer was one of the strongest present genome-wide (Fig. 3D) and the most correlated enhancer with the RARA gene by expression. Overall, the RARA enhancer ranked among the top 100 enhancers in 25% of samples (Fig. 3E). For 59% of the samples, the RARA enhancer was classified as an SE, suggesting that the AML blasts from these individuals have committed a significant amount of transcriptional machinery to this locus. Indeed, RARA expression was significantly higher in the top 25% of samples in which the RARA enhancer ranked in the top 100 enhancers and was correlated with enhancer score overall (Fig. 3F; Supplementary Fig. S5B). Although translocations involving the RARA gene are common in acute promyelocytic leukemia (APL; ref. 30), only 2 of the patients in our cohort harbored RARA fusions, and enhancer-based clustering did not place these samples into cluster 2 where RARA SE samples were enriched (Fig. 3B; Supplementary Fig. S5A). The RARA SE itself is located over the gene body of RARA in AML, MDS, and normal monocytes (Supplementary Fig. S5A and S5C) with acetylation spanning the active promoter and several introns and exons. Moreover, the size of the enhancer correlates with the expression of RARA in our patient samples, further supporting a direct link of this overlapping enhancer and gene. Importantly, published PCHi-C data (19) validate the three-dimensional interaction of the RARA promoter with the RARA SE (Supplementary Fig. S5C).

Given the similarities between AML and MDS, we also investigated the RARA locus in MDS and normal bone marrow CD34⁺ HSPCs. Analysis of a public data set (31) comparing gene expression in MDS CD34⁺ cells with normal CD34⁺ cells revealed that a subset of MDS cells have elevated RARA mRNA levels, similar to those seen in our AML data set (Supplementary Fig. S5D). By analyzing bone marrow aspirates from 2 patients with MDS, we detected RARA SEs with comparable score to those in AML (Supplementary Fig. S5E). In contrast, in normal bone marrow CD34⁺ cells, the RARA locus was not associated with an SE, and H3K27ac levels were similar to AML blasts that do not contain an SE at this locus.

Overall, the fact that the RARA SE is one of the largest SEs detected in a subset of patient samples raises the hypothesis that these cells may be dependent on the function of overexpressed RARα.

We asked whether the presence of a strong SE for RARA in non-APL AML cells could be predictive for sensitivity of such cells to RARα-targeted compounds and identified a series of models to test this (Supplementary Fig. S5F). As in primary patient samples, the score of the RARA enhancer and the corresponding mRNA expression varied across a collection of non-APL AML cell lines (Fig. 4A; Supplementary Table S4). Four cell lines with strong enhancers (MV4-11, OCI-AML3, Sig-M5, and EOL1) and 2 with weak enhancers (HEL and KG-1a), hereafter referred to as RARA-high and RARA-low, were used as initial models for differences in RARA SE score observed in our primary patient samples. RARα targeted agents were chosen for both selective agonism and antagonism. We first identified a potent and selective RARα agonist from the literature, tamibarotene (referred to as SY-1425 hereafter), which has previously been well validated as a selective and potent RARα agonist (14, 16, 17) and was further confirmed by us in biochemical and cellular assays (Supplementary Fig. S6A and S6B).

RARA-high and RARA-low cell lines were treated with a range of concentrations of SY-1425 and assessed for proliferation (Fig. 4B). Although SY-1425 was a potent inhibitor of proliferation in the RARA-high AML cell lines, with EC₅₀s that ranged from 0.14 nmol/L to 0.9 nmol/L, there was no effect on the RARA-low cell lines up to the highest concentration tested (2 μmol/L). A broader panel of non-APL AML cell lines further confirmed the positive correlation of both RARA mRNA (11 cell lines; Supplementary Fig. S7A) and enhancer score (13 cell lines; Supplementary Fig. S7B) with SY-1425 sensitivity, with responders showing significantly higher RARA enhancer score (one-sided t test, P = 0.03), and a correlation of enhancer score and mRNA levels (Supplementary Fig. S7C). To confirm that the effect of SY-1425 was indeed mediated by agonism of RARα, we tested the effect of a RARα-selective antagonist, BMS195614 (18), alone or in combination with SY-1425 on AML cells. The antagonist did not show any antiproliferative effect on RARA-high or RARA-low cell lines by itself (Supplementary Fig. S7D). Furthermore, coadministration of high concentrations of the antagonist with SY-1425 was found to completely rescue the effect of SY-1425 in the RARA-high cell line Sig-M5. Cotreatment with the RARα antagonist had no effect on growth inhibition by idarubicin, a DNA intercalator (Supplementary Fig. S7E). These results support the hypothesis that SY-1425 exerts its antiproliferative effects on AML cells specifically through activation of RARα.

To further strengthen our hypothesis that an SE at the RARA locus induces high levels of RARA transcripts, resulting in a cell state that is sensitive to a selective RARα agonist, we tested the effects of SY-1425 on 4 patient-derived mouse xenograft models of AML. The 4 PDX models are disseminated AML leukemia models where treatment is initiated when a 2% to 7% AML blast burden is detected in the peripheral blood of inoculated mice. This usually occurs about 20 to 40 days after inoculation with human tumor cells depending on the model. Mice were treated daily with 6 mg/kg SY-1425 orally, a dose that was chosen to match the human exposure when patients with APL are treated clinically with tamibarotene, and similarly studied for the nonselective retinoid ATRA (Supplementary Fig. S8A and S8B; refs. 32, 33). We used 2 PDX models (AM5512 and AM8096) where the tumor cells contained high levels of RARA mRNA and compared the effects of SY-1425 to 2 models where the tumor cells contained low levels of RARA mRNA (AM7577 and AM7440; Supplementary Fig. S9A). We observed significant reductions of AML blasts (measured using antihuman CD45 antibody) in peripheral blood in the 2 RARA-high models (Fig. 4C and D) but no effect of SY-1425 in the 2 RARA-low models (Fig. 4E and F). The reduction in circulating tumor cells in the RARA-high models was also associated with blast reduction in bone marrow and spleen (Supplementary Fig. S9B and S9C) and led to significant prolongation of survival compared with vehicle in these aggressive models of human AML (Supplementary Fig. S9D and S9E). In contrast, in the models with low RARA mRNA expression, no blast reduction in the bone marrow and spleen or prolongation of survival was observed (Supplementary Fig. S9F–S9I). This suggests that elevated levels of RARA mRNA could be used as a correlated surrogate measurement for the RARA SE to detect cancer cells predisposed to respond to SY-1425 in vivo.

The less-potent and less-selective RAR agonist ATRA, used clinically for treating APL, has had limited success in non-APL AML (34, 35). Indeed, ATRA was found to be less effective than SY-1425 in reducing proliferation of AML cell lines (Supplementary Fig. S10A). To explore whether the increased potency, selectivity, and significantly improved pharmacokinetics of SY-1425 versus ATRA are relevant in vivo, we compared the effect of SY-1425 to ATRA in the RARA mRNA-high PDX model AM5512. SY-1425 and ATRA were dosed twice daily with doses calculated to achieve exposure levels in mice similar to the human exposure observed in the clinic (refs. 32, 33, 36, 37; Supplementary Fig. S8A and S8B). SY-1425 greatly reduced tumor burden in the peripheral blood to <5% blasts, whereas ATRA treatment provided only a slight reduction in peripheral blast count compared with vehicle (Fig. 4G). Unlike SY-1425, the tumor burden in ATRA-treated animals failed to decrease by 2 weeks, at which point the mice were sacrificed due to weight loss and progressive disease. SY-1425 also reduced spleen and bone marrow tumor burden compared with vehicle and ATRA (Supplementary Fig. S10B and S10C). Finally, SY-1425, but not ATRA, prolonged survival compared with vehicle (Fig. 4H). These data are consistent with the less potent antiproliferative effect of ATRA on RARA-high cell lines compared with SY-1425 (7–21-fold less potent; Supplementary Fig. S10A) and improved pharmacology of SY-1425, including longer half-life, more stable repeat dosing exposure, and reduced sequestration.

To further investigate the mechanism of action of SY-1425 in the RARA-high model (AM5512), a study was performed in which bone marrow and blood were sampled weekly and analyzed by hematoxylin and eosin (H&E) staining and Ki67 IHC. Peripheral blood and bone marrow were found to closely mirror each other with respect to percent tumor burden (Supplementary Fig. S10D). At day 14, when the tumor burden in vehicle and treatment arms began to diverge, we observed marked morphologic differences between treated and untreated mice in both blood and bone marrow. In SY-1425–treated mice, we observed morphologic changes consistent with promotion of myeloid differentiation, including banded nuclei and increased cytoplasm to the nuclear area (Fig. 5A) compared with a high degree of poorly differentiated blasts in the vehicle-treated mice (Fig. 5B). In addition, Ki67, a marker of cellular proliferation, was greatly reduced in the bone marrow of SY-1425–treated mice compared with those treated with vehicle (Fig. 5A and B, inset).

To further investigate the hypothesis that SY-1425 induces differentiation in RARA-high, but not RARA-low, AML cells, we analyzed the response to SY-1425 in a number of cell line models in more detail. We found that upon SY-1425 treatment, known granulocytic differentiation and retinoic acid response genes (38–41) were upregulated in the RARA-high cell line MV4-11, but not in the RARA-low cell line OCI-M1. In particular, CD38, which is a marker of early myeloid maturation, was found to be strongly induced in RARA-high, but not RARA-low, cell lines at the mRNA (Fig. 5C) and protein level (Fig. 5D). Other myeloid differentiation markers such as ITGAX (CD11c; Supplementary Fig. S11A), ITGAM (CD11b; Supplementary Fig. S11B), and CD66 (Supplementary Fig. S11C) were more strongly induced in MV4-11 (RARA-high) and NB-4 (APL) than OCI-M1 (RARA-low).

We then sought to expand these findings beyond cell lines and analyzed 5 primary patient samples, 3 of which displayed high RARA mRNA levels (RARA-high) and 2 of which displayed low RARA mRNA levels (RARA-low). When treated ex vivo with 50 nmol/L SY-1425, a similar trend of CD38 maturation marker induction in RARA-high, but not RARA-low, cells was observed (Fig. 5E). Although some AML is classified as CD38-positive on the basis of CD38 dim expression, these were gated on a high mean fluorescent intensity level, which was achieved only after exposure to SY-1425. Comparing the induction in RARA-high versus RARA-low finds a significantly greater percentage of CD38-high phenotype cells in the RARA-high patients (Supplementary Fig. S11D).

To investigate the similarity of SY-1425 response in RARA-high AML to retinoid response in APL, we assayed the gene expression of 3 RARA-high cell lines and the NB-4 APL cell line following treatment with either SY-1425 or vehicle. At a global level, the transcriptional changes observed after SY-1425 were broadly consistent between the RARA-high cell lines and NB-4, an APL cell line (Fig. 6A). We next performed gene set enrichment analysis (GSEA) to query whether the transcriptional response of RARA-high cell lines to SY-1425 resembles gene sets of chemical or genetic perturbations found in the literature. Two of the most highly associated gene sets (FDR < 10⁻⁴) were relevant to APL pathogenesis or response to ATRA (tretinoin; Fig. 6B; Supplementary Table S5). We next compared the transcriptional response of RARA-high cell lines with changes previously described in APL cells treated with ATRA ex vivo (42). We found that the genes differentially regulated in RARA-high AML cell lines have largely consistent directionality with those in ATRA-treated APL patient samples (t test P < 10⁻¹⁶; Fig. 6C). To further characterize the transcriptional response to SY-1425, we performed GSEA on pathway gene sets and found that immune cell pathways were induced and MYC target genes were downregulated (FDR < 0.01; Supplementary Fig. S12A and Supplementary Table S6), consistent with the observed growth arrest and prodifferentiation effect. We conclude that the transcriptional response of APL cells to retinoids (ATRA or SY-1425) is similar to the response of RARA-high cell lines to SY-1425, suggesting an overlap in the mechanisms of action of these two agents in these two different forms of AML.

To further study the changes in epigenetic cell states induced by SY-1425, we performed ChIP-seq for RARα and H3K27ac in cell lines representing APL, RARA-high AML, and RARA-low AML. Focusing on the TGM2 locus, a well-described PML-RARα target gene (38), whose expression was induced in RARA-high AML and APL (Fig. 6D), we found an increase in H3K27ac in both NB-4 and the RARA-high line OCI-AML3 (Fig. 6E). There was a RARα binding site in both cell lines downstream of the promoter around which acetylation was induced (Fig. 6E). Globally, the majority of genes bound by RARα in NB-4 were also bound by RARα in the RARA-high cells and to a lesser degree in the RARA-low cells (Supplementary Fig. S12B). Because SE analysis has been hypothesized to reveal cell states (6, 7), we also performed a clustering analysis of the differential SEs following SY-1425 treatment in all 3 cell types. This analysis revealed that the RARA-high cells clustered more closely with NB-4 than with the RARA-low cells (Supplementary Fig. S12C and S12D). To examine whether the H3K27ac peaks induced by SY-1425 regulate known sets of differentiation genes, we used the GREAT method (43), which detects enrichment of genes near genomic regions of interest. We found strong overlap between the genes near enhancers induced in the RARA-high cells and genes induced by ATRA and bound by the PML–RARα fusion protein in APL cell lines (38; Fig. 6F; Supplementary Table S7). These results show that the response of enhancers in RARA-high AML cells to SY-1425 is similar to that associated with treatment of APL cells with ATRA.

In addition to the similarities with APL, we sought to better understand the transcriptional and epigenomic consequences of SY-1425 treatment. Therefore, the global transcriptional changes and chromatin alterations were examined to study the relationship between RARα target genes, chromatin activation, and mRNA regulation. Of particular interest was whether RARα required redistribution to activate genes upon treatment or showed switch-like behavior at its original positions. We observed large changes after treatment with SY-1425 in both gene expression and enhancer score in the RARA-high lines, but little to no changes in the RARA-low cell lines, suggesting that the transcriptional and epigenomic changes are due to RARα cell state and correlate with cell line sensitivity to SY-1425 (Fig. 7A; Supplementary Fig. S13A and S13B). Further supporting the role of RARα in the SY-1425 response, we found that genes and enhancers bound by RARα in the vehicle condition were more likely to be induced by SY-1425 treatment (Fig. 7B and C; Supplementary Fig. S13C). To further test whether part of the mechanism of action of SY-1425 could be through redistribution of RARα among potential binding sites, we investigated the difference in RARα binding between SY-1425 and vehicle treatment conditions in the RARA-high AML cell lines OCI-AML3 and MV4-11. Large changes in RARα binding locations with SY-1425 treatment were not observed in either RARA-high AML cell line (only 1/9,691 RARα peaks in OCI-AML3 and 480/14,919 in MV4-11 changed). Examining the genes that were near the RARα binding sites that changed in MV4-11, only 3 genes bound by RARα specifically upon treatment with SY-1425 were differentially expressed and only 2 genes bound by RARα specifically in the vehicle condition were differentially expressed (Fig. 7D). This supports the hypothesis that RARα undergoes a functional switch from repression to activation upon SY-1425 treatment rather than inducing a redistribution of RARα on the chromosome.

In order to find a pharmacodynamic marker for SY-1425 treatment, we investigated the gene expression changes in RARA-high cell lines upon treatment with SY-1425. DHRS3 showed the largest expression increase following SY-1425 treatment across RARA-high cell lines (Fig. 7A). This gene is implicated in endogenous retinoic acid metabolism and is responsive to retinoic acid levels (44–46). Furthermore, DHRS3 had multiple RARα binding sites both within and distal to the gene around which strong H3K27ac signals were observed following SY-1425 treatment (Fig. 7E). The positions of these RARα binding sites were largely unchanged upon SY-1425 treatment, supporting the hypothesis that new RARα binding sites are not required for the large induction of DHRS3 expression.

Here, we report the enhancer and transcriptional landscapes in 66 primary human AML patient samples. Using these patient-derived data sets, we identify SEs that are active in AML and define the transcriptional regulatory circuits that govern this aggressive malignancy. As a comparison, we mapped the enhancer landscapes of healthy hematopoietic stem, progenitor cells, and monocytes to provide a framework for hematopoietic differentiation states. SE analysis allows us to map the differentiation state of primary AML blast samples into the continuum of myeloid development. Applying nonnegative matrix factorization to our SE patient data set, we identified 6 distinct patient groups that were defined by the activity of enhancers driving expression of important TF genes such as IRF8, PROM1, MEIS1, and HOXA/HOXB.

To assess a potential clinical relevance of our SE-defined patient clusters, we made use of the well-curated transcriptomic and clinical data provided by TCGA. We projected our patient subgroup classification onto this large cohort to understand how enhancer landscapes might correlate with patient survival. This projection analysis demonstrated divergent overall survival among individual patient clusters. We note, however, that our AML cohort was biased toward high blast count AML which correlates with FLT3-ITD (29) mutation, and this analysis should be validated in additional patient cohorts in the future. The FLT3-ITD, although more frequent overall, was not correlated with the RARA SE.

Because we also analyzed the somatic mutation spectrum in our patient samples, we asked how mutations are distributed throughout our clusters and whether a similar clustering could have been achieved through exonic mutation analysis. Although some mutations were enriched or depleted in various clusters (for example, no NPM1 mutations were found in cluster 5 or 3, and MLL translocations were enriched in cluster 1), other mutations, such as those occurring in TET2, DNMT3A, and RAS, did not appear to be concentrated in a specific cluster. From these data, we conclude that clustering using SEs can provide orthogonal information that cannot be obtained by somatic mutational analysis alone.

SEs can also provide an avenue to identify novel druggable dependencies by finding those which are cancer specific or specific to certain tumor subtypes. By identifying SE-associated, cluster-specific TFs, we can begin to discern transcriptional regulatory circuitry. One SE-driven gene that appeared to be a key node in cluster 2 was the ligand-dependent TF RARA. The RARA SE was exceptionally strong in many cluster 2 patients, and RARα binding sites were enriched in cluster 2 SEs. The fact that this SE ranked as one of the largest enhancers in about 25% of the patient samples analyzed, while not present in normal, nonmalignant CD34-positive blasts, raises the possibility that this RARα circuitry represents a novel tumor liability that could be pharmacologically exploited.

Indeed, we found that a potent and selective RARα agonist (SY-1425) inhibited the proliferation of cell lines that contained the RARA SE, but not of cell lines that did not. This correlation was also true for RARA mRNA levels, presumably because SEs can drive high transcriptional activity for the genes they regulate (47). We verified this unique sensitivity of RARA-high AML cells in vivo in 4 AML PDX models where 2 RARA-low models showed no response to SY-1425, whereas SY-1425 treatment led to a significant reduction of tumor burden and prolongation of survival in the 2 RARA-high models. This observation does not extend to the natural retinoid ATRA, presumably due to its lower potency, selectivity, and PK properties (36, 48, 49). Compared with SY-1425, ATRA has a significantly less favorable PK profile in both mice and humans due to the fact that, as a natural derivative of vitamin A, ATRA is degraded by CYP26A1 and has been reported to induce the expression of this metabolizing enzyme, resulting in lower exposures upon repeat dosing (36, 48). SY-1425, as a synthetic retinoid, is resistant to CYP26A1 degradation (49) and has reduced binding to other retinoic acid sequestration proteins.

Investigating the molecular mechanism of SY-1425 in non-APL AML cell lines, we found that SY-1425 induced profound transcriptional changes in RARA-high, but not in RARA-low, cell lines (Fig. 7A). When we examined the distribution of RARα on chromatin and the alterations in H3K27ac before and after SY-1425 treatment, it became apparent that the chromatin location of RARα was not affected by treatment with SY-1425; rather, treatment tended to induce acetylation proximal to existing RARα ChIP-seq peaks. This is consistent with previously described activities of RARα, including acting in a transcriptional repressive complex while in the unliganded form and in a transcriptional activator complex when bound by an agonist (50). This supports a model in which chromatin-bound RARα nucleates enhancer formation and gene expression in response to the ligand SY-1425. Applying GSEA to the changed expression profile following SY-1425 treatment of RARA-high cells, we found an induction of immune signaling pathways (integrin, protein secretion, complement, MHC pathways) and a reduction in MYC target genes. These findings are consistent with the decreased proliferation and prodifferentiation effects observed in cell lines and xenograft models following SY-1425 treatment. Specifically, we found myeloid differentiation RARα target genes, like ITGAX (encoding CD11c), ITGAM (encoding CD11b), and CD38 to be induced in RARA-high, but not in RARA-low, AML cell lines. Moreover, histology of the bone marrow from RARA-high PDX models treated with SY-1425 revealed morphologic changes consistent with differentiation and a reduction in the proliferation marker Ki67 (Fig. 5A and B).

Because AML blasts that either contain a PML–RARA fusion (APL) or an SE at the RARA locus seem to both have a dependency on RARα, we compared the transcriptional response and epigenetic state in representative APL or RARA-high cell lines following RARα agonist treatment. We found the gene expression and enhancer responses of RARA-high AML cell lines to SY-1425 to be very similar to the responses of APL cells to ATRA or SY-1425. Gene sets of APL gene expression response to either treatment or genetic perturbation matched the response of RARA-high AML to SY-1425 at the gene expression and enhancer acetylation level. Across the genome, RARα binding was highly conserved between the APL cell line NB-4 and RARA-high AML cell lines and showed less similarity to the RARA-low AML cell lines. We conclude that many aspects of the well-studied action of retinoids on APL cells (42, 51–54) also apply to the treatment of RARA-high cells (defined here) with SY-1425. This includes the initial oncogenic cell state and the induction of differentiation and loss of proliferation after treatment with retinoids.

These findings support a model (Fig. 7F) where, just as in APL where the PML–RARA fusion protein acts as a transcriptional repressor driving the undifferentiated state of AML blasts, SE-induced overexpression of unliganded RARα acts as a transcriptional repressor of genes under the control of RARα. This assumes that the RARA SE and cell state are linked to super-physiologic evels of RARα or local depletion of endogenous retinoic acid such that the balance is distorted. The SE may also favor a more repressive isoform of RARα by altering the transcriptional start site of the RARA gene as noted with RARα2 (55). When APL cells are treated with ATRA, the PML–RARA fusion protein is degraded and the remaining wild-type RARα is induced by ATRA to function as a transcriptional activator (34, 56). Similarly, in the case of RARA-high AML cells, SY-1425 functions as an agonist turning RARα into a transcriptional activator. RARα binding to the genes and regulatory regions nucleates enhancer formation to drive upregulated expression of target genes in the presence of SY-1425. This potently reactivates differentiation pathways and drives the cells to terminal differentiation and growth arrest. Further validation of this model could be explored in future studies to demonstrate increased repressive complex association and activity with RARα in RARA-high cells. Alternatively, it could be conjectured that the RARA-high state merely potentiates the cells to gene activation without an active role in repression.

In summary, through the analysis of SEs in non-APL AML patient blasts, we have discovered a novel patient subgroup that is characterized by an SE at the RARA locus, predisposing cells to exquisite sensitivity to a potent and selective RARα agonist, SY-1425. SY-1425 acts through the induction of differentiation and arrest of proliferation in these cells using a mechanism of action similar to that of retinoids in APL. Although the discovery of this new cancer vulnerability was made possible by genome-wide SE analysis, there are multiple transcriptional readouts associated with the RARA SE that now can be used clinically to identify this patient population.

Previous clinical studies of retinoids in non-APL AML used a variety of patient selection methods, patient background types, and chemotherapy combinations (57–59). Individual responses in these trials and case reports are consistent with the possibility of a therapeutic benefit of retinoids in a subset of patients with AML and MDS distinct from APL. However, prior attempts may have suffered from the lack of an effective patient-selection strategy. In the case of ATRA, it is also possible that insufficient potency, nonselectivity over other RAR subtypes (in addition to other nuclear hormone receptors and sequestration proteins), short half-life, and diminished exposure after repeat dosing may have contributed to limited responses. Moreover, the increased reduction in blast count and survival seen with SY-1425 versus ATRA in our models is consistent with the results from clinical studies demonstrating increased complete-response rates of tamibarotene over ATRA in APL (32, 37, 60). The preclinical data presented here have led to the development of a clinical patient selection strategy utilizing an RARA SE-based biomarker for patient selection in a phase II study of SY-1425 in patients with non-APL AML and MDS (NCT02807558).

Normal donor human blood cells were obtained fresh from AllCells and FACS-enriched. Human AML samples were obtained from patients at the Stanford Medical Center with informed consent, according to Institutional Review Board (IRB)–approved protocols (Stanford IRB nos. 18329 and 6453). Mononuclear cells were isolated from each AML sample via Ficoll separation (GE Healthcare Life Sciences) and cryopreserved in liquid nitrogen in 90% FBS + 10% DMSO. All experiments presented here were conducted on freshly thawed cells. Criteria for inclusion of AML samples in this study were preestablished. Samples were included based on blast percentage at sample collection. For normal donors, no exclusion criteria were used. All AML samples were thawed in Iscove's Modified Dulbecco's Medium + 10% FBS + 200 U/mL DNase (Worthington Biochemical).

MDS patient samples were acquired from DVbiologics and processed as described in ChIP-seq methods (below). MDS sample 1 was from the bone marrow of a Caucasian female age 67. MDS sample 2 was from the bone marrow of a Caucasian female age 92.

We obtained the human AML cell lines OCI-AML3, Sig-M5, NOMO-1, and OCI-M1 from Deutsche Sammlung von Mikroorganismen und Zellkulturen GmbH (DSMZ). We obtained other human AML cell lines (HL-60, MV-4-11, THP-1, KG-1a, HEL, and Kasumi-1) from the ATCC. All cell lines were obtained during 2014–2015. Cell lines were incubated at 37°C under 5% CO₂, in RPMI-1640 plus 10% FBS and additional growth factors as necessary according to the tissue bank guidelines for each cell line. Mycoplasma testing was performed at least annually using the MycoAlert (Lonza) Detection Kit. Cell line identity was confirmed at the time of study using the ATCC Human STR Profiling Cell Authentication Service. Cell passage number was maintained below 20 culture splits.

Preparation cell viability/proliferation studies were assayed using ATPlite (Perkin Elmer) following the manufacturer's instructions and analyzed using GraphPad Prism version 6.0. Extended description of study protocol is included in the Supplementary Methods.

Preparation and analysis of ChIP-seq are described in full detail in Supplementary Methods. Raw data have been submitted to the Sequence Read Archive (SRA) under accession numbers SRP103200 (primary samples) and SRP103029 (cell lines).

AML samples were homogenized in 1 mL TRIzol and RNA isolated using RNAeasy (Ambion, ThermoFisher) according to the manufacturer's instructions. RNA libraries were prepared and sequenced as described previously (5). Analysis of RNA-sequencing (RNA-seq) data is described in detail in Supplementary Methods. Raw data have been submitted to the SRA under accession number SRP103200.

Expression studies were performed using TRIzol (Thermo Fisher Scientific) for extraction and mirVana (Thermo Fisher Scientific) purification kits for RNA isolation, and were analyzed on Affymetrix PrimeView arrays. Full details of preparation and analysis are provided in the Supplementary Methods. Raw data have been submitted to GEO under the accession GSE97331.

The PDX studies were performed at Crown Bioscience using their HuKemia systemic models. Raw RNA-seq data on the models were provided by Crown Bioscience and analyzed as per the RNA-seq methods above. FACS analysis was performed by Crown Biosciences. The study protocols and procedures involving the care and use of animals were reviewed and approved by the Institutional Animal Care and Use Committee of Crownbio prior to conduct. During the studies, the care and use of animals was conducted in accordance with the regulations of the Association for Assessment and Accreditation of Laboratory Animal Care. Slides for H&E and Ki67 were made from formalin-fixed, paraffin-embedded blocks from Crown Biosciences that were prepped and stained at the Specialized Histopathology Services (Dana-Farber/Harvard Cancer Center, Brigham and Women's). Additional methods for in vivo studies are provided in Supplementary Methods.

M.R. McKeown has ownership interest (including patents) in Syros Pharmaceuticals. M.L. Eaton has ownership interest (including patents) in Syros Pharmaceuticals. C. Fiore has ownership interest (including patents) in Syros Pharmaceuticals. E. Lee has ownership interest (including patents) Syros Pharmaceuticals. M.W. Chen has ownership interest (including patents) in Syros Pharmaceuticals. M.G. Guenther is a stockholder and patent inventor at Syros Pharmaceuticals. D.A. Orlando has ownership interest (including patents) in Syros Pharmaceuticals. C.C. Fritz has ownership interest (including patents) in Syros Pharmaceuticals. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.R. McKeown, M.R. Corces, M.L. Eaton, C. Fiore, K. Austgen, M.G. Guenther, J. Lovén, C.C. Fritz, R. Majeti

Development of methodology: M.R. McKeown, M.R. Corces, M.L. Eaton, C. Fiore, E. Lee, M.W. Chen, M.G. Guenther, D.A. Orlando, J. Lovén

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.R. McKeown, M.R. Corces, E. Lee, M.W. Chen, J.L. Koenig, K. Austgen, J. Lovén, R. Majeti

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.R. McKeown, M.R. Corces, M.L. Eaton, C. Fiore, D. Smith, S.M. Chan, K. Austgen, D.A. Orlando, J. Lovén, C.C. Fritz, R. Majeti

Writing, review, and/or revision of the manuscript: M.R. McKeown, M.R. Corces, M.L. Eaton, C. Fiore, M.G. Guenther, J. Lovén, C.C. Fritz, R. Majeti

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.L. Eaton, J.T. Lopez, M.W. Chen, D. Smith, J. Lovén

Study supervision: M.G. Guenther, J. Lovén, C.C. Fritz, R. Majeti

Thank you to Dr. Jason Marineau and Dr. Kevin Sprott for help on chemistry review and sourcing, and to Dr. Eric Olson for guidance. Thank you to Dr. Abner Louissaint for advice and interpretation of histopathology.

This study was supported by Stinehart-Reed Foundation (R. Majeti), Ludwig Center for Cancer Stem Cell Research (R. Majeti), and the NIH (R01CA18805 to R. Majeti). R. Majeti is a New York Stem Cell Foundation Robertson Investigator and Leukemia and Lymphoma Society Scholar.