Multi-omics profiling reveals a distinctive epigenome signature for high-risk acute promyelocytic leukemia

INTRODUCTION

Acute promyelocytic leukemia (APL) is a subtype of acute myeloid leukemia (AML) and accounts for ~10% of all AML cases. APL is most often characterized by the chromosomal translocation t(15;17) that fuses the PML gene on chromosome 15 to the retinoic acid receptor α (RARA) gene on chromosome 17 resulting in the transcription of the PML-RARA oncofusion protein [1, 2]. The PML-RARA fusion has been reported to be present in >95% of APL cases [1, 3]. In normal cells, PML regulates a wide range of proteins including transcription factors and ensures normal functioning of biological processes such as DNA damage response and microorganism resistance [4], while RARA-involved pathways are critical for terminal differentiation of myeloid cells [5]. The formation of the PML-RARA oncofusion protein disrupts normal functioning of both fusion partners expressed by wild-type alleles [6] and induces a maturation block at the promyelocytic stage [7].

The PML-RARA oncofusion protein has been proposed to recruit epigenetic modifiers, like histone deacetylase complexes such as SMRT [8] and N-CoR [9], DNA methyltransferases [10], and the repressive histone methyltransferases SUV39HI [11] and PRC2 [3], to provoke alterations in chromatin structure. Genome-wide analyses revealed that the majority of PML-RARA binding sites have low histone H3 acetylation [12–14], suggesting a main role for histone deacetylases in PML-RARA leukemogenesis.

Although the PML-RARA fusion event is a major driving force in APL, the treatment of APL patients is clinically different according to risk classification. In general, APL can be subdivided into three risk categories based on white blood cell (WBC) and platelet count. Low risk APL patients have a WBC count lower than 10,000/μL and a platelet count over 40,000/μL, intermediate-risk APL patients carry a WBC count under 10,000/μL but a platelet count lower than 40,000/μL, while high-risk APL patients have a WBC count >10,000/μL. Despite significant advances in the treatment of APL, the category of high-risk APL still poses pressing challenges.

Pharmacological doses of All-trans retinoic acid (ATRA) and/or arsenic trioxide (ATO) are widely used by clinicians in treating low- and intermediate-risk APL patients [15, 16]. ATO and ATRA act on PML and RARA, respectively, and promote degradation of the PML-RARA oncofusion protein. ATRA adopts the proteasome-mediated pathway [17] and caspases [18] to degrade PML-RARA, and ATO induces sumoylation of the PML moiety [19]. Depletion of the PML-RARA oncofusion protein leads to loss of recruitment of epigenetic repressors and remodels the repressive chromatin environment at PML-RARA binding sites, creating a more accessible chromatin state [13]. In contrast to low/intermediate-risk APL patients, high-risk APL patients generally show a higher relapse rate, and are currently often treated with cytotoxic chemotherapy in combination with ATRA/ATO [20].

By performing an in-depth transcriptomic and epigenomic analysis, we aimed to gain further insights into these high-risk APLs (hrAPLs) that might ultimately help in improving their prognosis. We used ChIP-seq to identify H3K4me1, H3K4me3, H3K27ac, H3K9/14ac, H3K27me3, and H3K36me3 enriched regions in patient blasts from two high-risk APLs, and also conducted whole-genome-bisulfite sequencing (WGBS) and RNA-seq to probe DNA methylation and gene expression profiles, respectively. We compared the epigenomes of high-risk APL samples with non-high-risk APL samples as well as with other progenitor cells from the normal myeloid compartment, revealing a hyperacetylation signature at a defined set of enhancer regions in hrAPLs. As we also uncovered hrAPL-specific alterations for H3K27me3, our analysis suggested that a combination of epigenetic drugs that affects acetyltransferase and H3K27 methyltransferase activity might be beneficial. Using the novel HAT/EZH2 inhibitor MC2884 confirmed that altering these hrAPL-specific epigenomic features could indeed significantly inhibit cell growth.

RESULTS

The genomic spectrum of two high-risk APLs

To investigate the molecular determinants of high-risk APLs, we analyzed the genome and epigenome of two rare primary APL samples of patients that were resistant to standard ATRA plus chemo treatment and presented extremely high white blood cell counts (>10,000/μL) (Figure 1A). We confirmed transcription of the PML-RARA fusion in both samples, which revealed expression of the BCR2 isoform in patient 7 (pt#7) and BCR1 in patient 8 (pt#8) (Figure 1B). In addition, we used SNP calls inferred from WGBS data to look for common gene variants in the RARA moiety of the two hrAPLs, that could potentially explain the absence of an ATRA response [21–24]. We found one rare gene variant signature with minor allele frequency (MAF) <0.007 in our two hrAPL samples, but as this variant was synonymous it was not expected to have functional effects on the RARA moiety.

To examine the full mutational landscape of the two hrAPL patients, we used variant calls from WGBS as our primary calls and applied additional thresholds of read depth (DP) ≥20 and Phred score ≥40. To have more confidence in our variant calls, we confirmed these calls using RNA-seq and/or ChIP-seq profiles of these two samples. Next, we looked for only missense variants and evaluated them using the CHASM program (http://www.cravat.us/), on a cut-off of MAF ≤0.05, which resulted in the identification of a range of mutations. Specifically, the analysis identified several candidate driving mutations such as TET2 and BRCA2, but none of them was recurrent in the two patients (Figure 1C).

A hyperacetylation signature typifies PML-RARA binding sites in high-risk APLs

Given the absence of common secondary mutations, we focused on the epigenomic patterns observed in the two hrAPL patients and explored whether an ATRA response was still observed at the epigenetic level as compared to low/intermediate-risk APLs. The APL-associated oncofusion protein PML-RARA is able to recruit HDACs and establish a hypoacetylated chromatin signature, which can be alleviated by ATRA-induced PML-RARA degradation [13]. Hence, we first surveyed binding of PML-RARA in hrAPLs by ChIP-seq using an antibody against PML (the non-DNA-binding moiety of the PML-RARA fusion). Our analysis revealed occupancy of PML-RARA at previously identified high-confidence PML-RARA binding sites (n = 2,723) (Figure 2A) from the NB4 model cell line and low/intermediate-risk APLs [13], such as the TGM2 and ID1 genomic regions (Figure 2B). We generated additional epigenetic ChIP-seq profiles using ex vivo ATRA treated patient cells and examined epigenetic alterations at the same PML-RARA binding sites. This confirmed elevated histone H3 acetylation levels after PML-RARA depletion by ATRA treatment at the TGM2, RARB, ASB2 and ID1 loci (Figure 2C). Next, we further delved into global ATRA response at all 2,723 PML-RARA occupancy loci and found that generally H3K9/14ac levels increased in an ATRA-treatment-dependent manner in hrAPLs (pt#7 and pt#8) (Figure 2D). Interestingly, we noted significantly elevated steady-state levels of histone acetylation at PML-RARA binding sites in hrAPLs as compared to low/intermediate-risk APL (pt#9) (Figure 2D) in both control and ATRA-treated samples, suggesting that abnormal hyperacetylation at PML-RARA binding sites might be a crucial property of high-risk APLs.

Hyperacetylation at enhancer and promoter regions in high-risk APLs

To analyze whether increased acetylation (H3K9/14ac and H3K27ac) was a genome-wide feature, we identified all active hrAPL promoter and enhancer regions based on levels of the promoter mark H3K4me3 or the enhancer mark H3K4me1 in association with H3K27ac (Figure 3A). A total of 9,349 common active promoters (high H3K4me3 and H3K27ac, low H3K4me1 and H3K27me3) and 10,208 common active enhancers (high H3K4me1 and H3K27ac, low H3K4me3 and H3K27me3) were detected. In line with the positive correlation between histone H3 acetylation levels and gene expression, we observed that associated genes of active promoters and enhancers are expressed, but also that these target genes seem higher expressed in hrAPLs as compared to APL (Figure 3B). To further explore whether elevated acetylation levels are common, we compared H3K9/14ac intensity for all active promoter and enhancer regions in two hrAPLs and two APLs. Our findings revealed increased H3K9/14ac levels at the two types of regions in hrAPLs as compared to APLs (Figure 3C), suggesting the hyperacetylation signature is not only present at PML-RARA binding sites, but even is a ubiquitous mark at other regulatory elements.

As the difference in histone H3 acetylation was most obvious at putative enhancer regions, we examined whether these regions might represent an activation signature of a unique cellular population. For this, we compared the H3K27ac occupancy landscape of hrAPLs with those of five progenitor cells from healthy individuals including normal CD34⁺ cells, promyelocytes, metamyelocytes, band neutrophils and segmented neutrophils. This revealed distinct global acetylation levels at the hrAPL enhancer sites in different progenitor population (Figure 3D), but no similar enrichments as observed in hrAPLs, suggesting that the hyperacetylation signature present in hrAPLs is not reminiscent of a specific progenitor population.

To determine enriched consensus sequences at all active enhancers with the hyperacetylation signature in hrAPLs (Figure 3C), we conducted motif screening. This discerned significant enrichment for RUNX1 and ETS consensus motifs (Figure 3E), implying transcription factors belonging to these families might assist in epigenetic marking of these regions.

Finally, we asked whether the hyperacetylation state at identified enhancers is specific for hrAPLs. For this, we compared the hrAPL H3K27ac pattern with levels of 27 AMLs with varying prognosis and cytogenetic background [25]. In order to get unbiased and precise results, we included 4 additional high-risk APLs (pt#A to D; number of leukocytes >10,000/μL) and 1 low-risk APL (pt#E). Indeed, acetylation levels at the set of hrAPL enhancers are increased as compared to acetylation at these regions in almost any other subtype of AML (p < 0.01 for pairwise comparisons) (Figure 3F), implying that the acetylation signature might be used to identify this particular subclass of APLs.

To allow clinical discrimination of this subset using assays that are sequencing-independent, we defined the most differential acetylated regions. Using an FDR <0.01 and p-value < 0.0002, we identified 22 regions specifically acetylated in hrAPLs that could potentially be used for clinical subtype stratification (Supplementary Table 1), although these results will need further confirmation in a larger cohort of hrAPLs.

Altogether these results reveal that the histone acetylation landscape of hrAPLs is different as compared to low/intermediate-risk APLs and myeloid progenitors, and generally represents a hyperacetylated histone signature.

Deregulation in repressive histone methylation impacts on APLs with low overall survival

APL leukemogenesis has also been linked to H3K27me3 alterations [3, 13], suggesting an important role for those genomic regions occupied by the PRC2 complex and H3K27me3 in APL. We speculated that in hrAPLs, the hyperacetylated signature might affect the presence of repressive marks, especially H3K27me3 given the mutual exclusivity of acetylation and methylation on the same lysine. Inspecting H3K27me3 patterns in hrAPLs and APL revealed differential occupancy at a number of loci, for example, H3K27me3 hypermethylation at the ZSCAN12/GPX6 genomic region (Figure 4A, 4B), but mainly hypomethylation at other (n = 214) differential H3K27me3 regions in hrAPLs. These results indicate that besides H3 acetylation also other epigenetic features such as H3K27me3 are altered in hrAPL and likely affect the transcriptional program.

High-risk APLs harbor an altered DNA methylation signature

Next, we assessed whether hrAPL and APLs presented differential global CpG methylation patterns. To extend our dataset, we included DNA methylation data of a cohort of 13 APL patients (450 K, from TCGA; [26]), of which 3 had an exceptionally high white blood cell count (>10,000/μL) and survival below 1 month that could be integrated as hrAPL cases (Figure 5A, 5B).

Supervised clustering to identify differential CpG methylation sites between hrAPLs and APLs revealed a set of 1,817 CpGs that could classify APLs into the high- or low/intermediate-risk category (Figure 5C). These included 399 hypomethylated CpG regions (CpG coverage >4 and p-value < 0.005) which displayed enriched H3K4me3 and H3K27ac levels in hrAPLs (Figure 5D), and 1,418 hypermethylated CpG regions. To further validate this finding, we integrated all DNA methylation data from TCGA AML samples [26]. This indicated that the unique hyper/hypo methylation signature discovered for hrAPLs is specific (Figure 5E, top). Importantly, this signature could be further filtered to 23 CpGs (CpG coverage > 4, p-value < 0.0007) (Supplementary Table 2; Figure 5E, bottom), which could distinguish hrAPLs from all APLs/AMLs.

To validate this signature, we examined methylation levels reported for high-, intermediate- and low-risk APLs in a separate study [27]. This showed that 3 out of 4 regions with sufficient coverage were hypomethylated in high-risk cases as compared to low/intermediate-risk APLs (Figure 5F). These results suggest that the 23 CpG hrAPL signature potentially enables PCR-based approaches to identify the hrAPL subtype at diagnosis. As for the identified Histone H3 acetylation signature above, further confirmation in a larger cohort of hrAPLs will be needed.

MC2884 induces apoptosis in hrAPLs and remodels the epigenome

Since hrAPLs display multiple epigenome deregulations, we wondered whether a multi epi-drug approach targeting HATs and PRC2 would be beneficial. We identified MC2884, a novel HAT/EZH2 inhibitor (Petraglia et al., in press), as a promising candidate. We evaluated MC2884-induced apoptosis in the hrAPLs (pt#7 and pt#8) as well as in an ATRA-sensitive APL (pt#9) after 24 hours ex vivo treatment. Our results show that the hrAPLs did not undergo efficient apoptosis (as measured by the % of cells in the pre-G1 phase) ex vivo in response to HDACi’s such as SAHA or MS275 [28], when compared to ATRA-sensitive blasts (pt#9). In contrast, ex vivo MC2884 treatment in hrAPLs was accompanied by remarkable apoptosis induction (Figure 6A).

As the hrAPLs harbor an aberrant H3K27ac, H3K9/14ac and H3K27me3 signature, we inquired whether MC2884 would reset acetylation and H3K27me3 patterns and thereby trigger the onset of apoptosis. We therefore treated pt#8 blasts with MC2884 ex vivo for 4 and 24 h, and examined H3K27ac and H3K27me3 levels. We identified 3,636 H3K27ac regions (Figure 6B) for which occupancy was changed (FC > 3×STDEV, % gain, % loss), such as at the BCL2, SEC23 and AK2 genomic regions (Figure 6C). Amongst the changed regions, 1,387 represented enhancer regions decreased in H3K27ac (Figure 6D). Functional annotation of these decreased enhancer regions revealed enrichment for genes associated with AML, PML-RARA binding, JAK/STAT signaling, Toll-like receptor signaling, and the RAS pathway (Figure 6E).

In addition to H3K27ac changes, we observed 799 H3K27me3 regions for which occupancy was changed (FC > 3×STDEV, % gain, % loss) (Figure 6F), suggesting that MC2884 treatment might alter the function of histone H3 lysine 27 both by reducing acetylation at this residue as well as methylation. Interestingly, reduced H3K27ac was also observed at a putative enhancer region of EZH2 (Figure 6G), suggesting that MC2884-chromatin action also affects EZH2 gene transcription and further changes polycomb complex localization. Indeed, expression of EZH2 was reduced, similar as for other genes associated with reduced enhancer acetylation (Figure 6H, left). In contrast, genes associated with reduction in H3K27me3 increased in expression (Figure 6H, right). Together, these results reveal that MC2884 treatment reshapes the epigenome of hrAPLs and increases apoptosis.

DISCUSSION

Treatment of APL has significantly advanced over the last decade due to increased knowledge on the molecular mechanisms underlying APL pathogenesis and the widespread application of ATRA and arsenic trioxide (ATO) in treatment protocols. While low/intermediate-risk APL patients have benefited from these developments, especially the reduced usage of chemotherapy, clinical therapy for high-risk APL patients often still include cytotoxic chemotherapy [20]. Here, by performing an in-depth transcriptomic and epigenomic analysis, we aimed to point out major causes of severe drug resistance in these high-risk APLs (hrAPLs) and identify entry points for the application of epi-drugs.

As no common genetic signature could be identified by our genomic scan, we focused on the epigenomic comparison between high-risk and low/intermediate-risk APL samples. Our findings revealed hrAPL-specific alterations in histone acetylation, H3K27me3 and DNA methylation. Given these specific epigenomic features, we tried to establish hrAPL-specific epigenetic signatures that could potentially be used as biomarker. For histone acetylation, we identified 22 hyperacetylated regions that seemed to be hrAPL-specific, and for DNA methylation we proposed a signature of 23 CpGs that distinguishes hrAPLs from other APLs/AMLs. Additional samples and confirmation in different cohorts will be needed to support these regions as real biomarkers. Nevertheless, it suggests that histone acetylation level has similar potential to classify samples as DNA methylation signatures.

The prominent signature of hyperacetylation of histones at a set of defined promoters and enhancers in high-risk APLs suggests a more open chromatin structure in hrAPLs. The hyperacetylation signature was also observed at PML-RARA binding sites, which seemed counterintuitive given the previously reported recruitment of histone deacetylases by us and others [8, 9, 13, 29]. Moreover, it is surprising that despite a hyperacetylation signature at PML-RARA binding sites, histone acetylation is still increased upon ATRA treatment. Previously, we have shown that for another AML associated fusion protein (AML1-ETO) both acetyltransferase and deacetylase enzymes are present at its binding site [30]. This could well be the case in this situation as we already previously reported the acetyltransferase p300 might be present at PML-RARA binding sites [31] to mediate histone acetylation. Maybe through stabilization of its binding or increased activity, hyperacetylation at the PML-RARA binding regions is established in hrAPLs. However, it can still be expected that also HDACs recruited by PML-RARA will be present. These will be removed upon ATRA treatment and subsequent PML-RARA degradation, resulting in an even further increase in acetylation.

Given the hyperacetylation signature and the absence of a clear ATRA response (which would even further increase acetylation), inhibition of HDAC activity which has been previously proposed to treat APL [9, 32–36] would likely be ineffective. As we also observed other epigenetic alterations in hrAPLs as compared to APLs, especially for H3K27me3, our analysis suggested that a combination of epigenetic drugs that affect histone acetyltransferase and H3K27 methyltransferase activity might be beneficial. Using the novel HAT/EZH2 inhibitor MC2884 confirmed that altering these hrAPL-specific epigenomic features could indeed inhibit cell growth. These results suggest that epigenome-based patient stratification, in this case hrAPLs may provide a novel strategy to be coupled with genome analyses to “personalize” precision medicine.

METHODS

Next generation sequencing (NGS) experiments

Samples for NGS experiments were obtained from different laboratories but have been processed similarly, with standard operating procedures set by BLUEPRINT at the start of the project using extensively validated and documented antibody batches and other reagents. Moreover, all samples underwent strict quality control at different levels (e.g. sequence quality, read depth, visual inspection) before they were included in the downstream analysis pipeline. The exact procedures and quality filters can be found at the BLUEPRINT website (www.blueprint-epigenome.eu).

Chromatin immunoprecipitation (ChIP)

ChIP extracts preparation and procedures have been carried out following IHEC procedures using Diagenode antibodies and as reported [13, 29]. The fold enrichment of H3K27ac and H3K27me3 ChIP’ed DNA was evaluated. Primer sequences were as follows: BCL2 promoter region (at -410 and -282 from +1), forward (5ʹ- GTG TTC CGC GTG ATT GAA GAC-3ʹ) and reverse (5ʹ- CAG AGA AAG AAG AGG AGT TAT AA-3ʹ); for chr21:44,496,353-44,496,918/CBS region, forward (5ʹ- CGC AGA ACA GTC GCC TTG-3ʹ) and reverse (5ʹ- GTC CAG AGC ACG ATG TTT GG-3ʹ); for chr17:4,938,577-4,939,058/SLC52A1 region, forward (5ʹ- CGA GTT GGA GAG GGG AGT G-3ʹ) and reverse (5ʹ- AAC AAA ACC CCA GCT GTG TG-3ʹ).

RNA extraction and RT-PCR

Total RNA was extracted with Trizol (Invitrogen) and converted into cDNA using VILO (Invitrogen). For amplification, the following primers were used: EZH2, forward (5ʹ-CATCATAGCTCCAGCTCCCG-3ʹ) and reverse (5ʹ-CATCCCGGAAAGCGGTTTTG-3ʹ); EED, forward (5ʹ-CGATTTGCGACAGTGGG-3ʹ) and reverse (5ʹ-CAGGTGCATTTGGCGTG-3ʹ); SUZ12, forward (5ʹ-GTCCTGCTTGTGAAAGTTTGC-3ʹ) and reverse (5ʹ-CAAATGTCTTTTCCCCATCCT-3ʹ); BCL2, forward (5ʹ-GAACTGGGGGAGGATTGTGG-3ʹ) and reverse (5ʹ-CAGCCTCCGTTATCCTGGAT-3ʹ); GAPDH, forward (5ʹ-TCAACGGGAAGCCCATCACCA-3ʹ) and reverse (5ʹ-ACGGAAGGCCATGCCAGTGA-3ʹ).

Whole-genome bisulfite sequencing

Genomic DNA (1–2 μg) was spiked with unmethylated λ DNA (5 ng of λ DNA per μg of genomic DNA) (Promega). The DNA was sheared by sonication to 50–500 bp using a Covaris E220 and fragments of size 150–300 bp were selected using AMPure XP beads (Agencourt Bioscience Corp.). Genomic DNA libraries were constructed using the Illumina TruSeq Sample Preparation kit (Illumina Inc.) following the Illumina standard protocol: end repair was performed on the DNA fragments, an adenine was added to the 3ʹ extremities of the fragments and Illumina TruSeq adapters were ligated at each extremity. After adaptor ligation, the DNA was treated with sodium bisulfite using the EpiTexy Bisulfite kit (Qiagen) following the manufacturer’s instructions for formalin-fixed and paraffin-embedded (FFPE) tissue samples. Two rounds of bisulfite conversion were performed to assure a high conversion rate. An enrichment for adaptor-ligated DNA was carried out through 7 PCR cycles using the PfuTurboCx Hotstart DNA polymerase (Stratagene). Library quality was monitored using the Agilent 2100 BioAnalyzer (Agilent), and the concentration of viable sequencing fragments (molecules carrying adaptors at both extremities) estimated using quantitative PCR with the library quantification kit from KAPA Biosystem. Paired-end DNA sequencing (2 × 100 nucleotides) was then performed using the Illumina Hi-Seq 2000.

The data are publicly available and can be accessed via the Blueprint DCC: http://dcc.blueprint-epigenome.eu/#/experiments.

Peak calling and identification of promoter and enhancer regions

For peak calling the BAM files were first filtered to remove the reads with mapping quality less than 15, followed by fragment size modeling (http://code.google.com/p/phantompeakqual-tools/). The peak-calling algorithm MACS2 (http://github.com/taoliu/MACS/) was used to detect the binding sites for the three studied histone marks at default q-value (5.00e-02). H3K4me1 peaks were called using the broad setting of MACS2 while peaks for marks H3K27ac, H3K4me3 and transcription factor PML-RARA were called using the default (narrow) setting.

In order to identify promoter and enhancer regions, non-overlapping H3K4me3 and H3K4me1 bound regions were taken into account, respectively. The common promoter and enhancer regions across the hrAPL samples were computed using the in-house script intersectbed. pl. Further, normalized tag enrichment for H3K27ac and H3K27me3 was used to mark active and inactive promoter and enhancer regions. Also, normalized tag enrichment for H3K27ac and H3K9/14ac in common enhancer regions was used to compute differentially active regions across different samples (pt#7, pt#8 and pt#9) and cell types (CD34⁺, promyelocytes, metamyelocytes, band neutrophils and segmented neutrophils). The intensity graphs were generated using the in-house script makeColorprofiles.pl. Discriminating hyperacetylated regions were identified using DESeq package (http://bioconductor.org/packages/release/bioc/html/DESeq.html) and selection for regions with an FDR < 0.01 and p-value < 0.0002 between hrAPL and all other samples (APL and AML). Motif analysis was performed as before [30].

The data are publicly available and can be accessed via the Blueprint DCC: http://dcc.blueprint-epigenome.eu/#/experiments.

RNA-seq alignment and expression analysis

RNA-seq reads were aligned using GSNAP [37] using non-default parameters -m 1 - N 1 -n 1 -Q -s Ensembl_splice_68. RNA-seq library data was initially subjected to a quality control step, where, based on read distribution over the annotated genome, libraries that are outliers were identified and discarded from further analysis. For expression analyses, reads were aligned to the Ensembl v70 (GRCh37.70) human transcriptome using Bowtie. Quantification of gene expression was performed using MMSEQ [38].

The data are publicly available and can be accessed via the Blueprint DCC: http://dcc.blueprint-epigenome.eu/#/experiments.

SNP calling and identification of potential functional SNPs

SNP calls were identified using the CNAG in-house pipeline based on WGBS data. The list of variant calls was reduced by applying a number of filters, starting from the number of tags (≥20) and Phred score (≥40) on the SNP calls from WGBS. The reduced list was then given as an input to the CHASM web server (http://www.cravat.us/ version 3.1) for functional annotation and only the missense variants present in driver genes (as identified by Vogelstein) [39] at MAF ≤ 5% (1000 G) were selected. Further, only the variants reported in ChIP-seq and RNA-seq tracks by applying computational tools samtools (mpileup option) and GATK along with the variations confirmed through visual inspection of the tracks using IGV were selected.

The data are publicly available and can be accessed via the Blueprint DCC: http://dcc.blueprint-epigenome.eu/#/experiments.

Identification of signature methylation regions

Cytosine methylation level was called for every cytosine site in different sequence contexts such as CHH, CXG and CG (H denotes A, T, or C and X denotes A or T). Due to the fact that symmetric methylation level is the common case for CG methylation, we focused on CG methylation and excluded non-CpG methylation for further analysis. Within a symmetric CG context, CG pairs on both forward and reverse strands were combined and a CpG coverage threshold of ≥5 was used. Further, the methylation samples pt#7 and pt#8 were grouped together with the 3 APL samples of low OS from TCGA [26] and compared with the APL samples. We selected the regions for which CG methylation was determined across all samples and applied the Wilcoxon test to determine statistically significant different CpGs between the two groups. A cut-off of the p-value of 0.005 was applied as the selection criteria and regions were plotted. To identify the most discriminating CpGs, we applied the same test at a more stringent cut-off p-value of 0.0007. Moreover, only the CpG regions that had beta values in all 181 AML samples were selected.

The metadata of the APL and AML samples can be found in Supplementary Table 3.

Cell cycle analysis

For cell cycle analyses, cells were plated (2 × 10⁵ cells/mL) and after stimulation (performed as indicated in the text) were harvested, centrifuged at 1,200 rpm for 5ʹ and resuspended in 500 μL of a hypotonic solution containing 1X PBS, Sodium Citrate 0.1%, 0.1% NP-40, RNAase A and 50 mg/mL Propidium Iodide (PI). After 30’ at room temperature (RT) in the dark, samples were acquired by FACS-Calibur (BD Biosciences, San Jose, CA, USA) using CellQuest software (BD Biosciences). The percentage in different phases of the cell cycle (and pre-G1 as a measure of apoptosis) was determined by ModFit LT V3 software (Verity). All experiments were performed in triplicate.

Author contributions

FP, AAS: data generation and analysis of the majority of the experiments; AN, MC, SV, LS: data interpretation and analyses; AM, RL, HK. EMJM, FP, EH, KB, BK, GY, CL, SH, LC, PF, TK, MLY, VDV, OB, IG, HGS, AM: genome/epigenome data generation and analyses; JHJ, EV, ATJW: AML collection and clinical information incorporation; JHAM & LA: conceived the study and wrote the manuscript.

ACKNOWLEDGMENTS AND FUNDING

A special thanks to the generosity of the patients whose samples allowed the present study and discoveries. This work was supported by BLUEPRINT (282510), the Italian-Flag Project-EPIGEN, PRIN-20152TE5PK, the Italian association against cancer (AIRC-17217), the Dutch Cancer Foundation (KWF kun 2011–4937), KIKA (project 311) and Netherlands organization for Scientific Research (NWO-VIDI to JM).

CONFLICTS OF INTEREST

The authors declare that they have no conflicts of interest.

REFERENCES

1. de The H, Chomienne C, Lanotte M, Degos L, Dejean A. The t(15;17) translocation of acute promyelocytic leukaemia fuses the retinoic acid receptor alpha gene to a novel transcribed locus. Nature. 1990; 347:558–61. https://doi.org/10.1038/347558a0.

2. Kakizuka A, Miller WH Jr, Umesono K, Warrell RP Jr, Frankel SR, Murty VV, Dmitrovsky E, Evans RM. Chromosomal translocation t(15;17) in human acute promyelocytic leukemia fuses RAR alpha with a novel putative transcription factor, PML. Cell. 1991; 66:663–74.

3. Villa R, Pasini D, Gutierrez A, Morey L, Occhionorelli M, Vire E, Nomdedeu JF, Jenuwein T, Pelicci PG, Minucci S, Fuks F, Helin K, Di Croce L. Role of the polycomb repressive complex 2 in acute promyelocytic leukemia. Cancer Cell. 2007; 11:513–25. https://doi.org/10.1016/j.ccr.2007.04.009.

4. Lallemand-Breitenbach V, de The H. PML nuclear bodies. Cold Spring Harb Perspect Biol. 2010; 2:a000661. https://doi.org/10.1101/cshperspect.a000661.

5. Gratas C, Menot ML, Dresch C, Chomienne C. Retinoid acid supports granulocytic but not erythroid differentiation of myeloid progenitors in normal bone marrow cells. Leukemia. 1993; 7:1156–62.

6. Brown NJ, Ramalho M, Pedersen EW, Moravcsik E, Solomon E, Grimwade D. PML nuclear bodies in the pathogenesis of acute promyelocytic leukemia: active players or innocent bystanders? Front Biosci (Landmark Ed). 2009; 14:1684–707.

7. Wang ZY, Chen Z. Acute promyelocytic leukemia: from highly fatal to highly curable. Blood. 2008; 111:2505–15. https://doi.org/10.1182/blood-2007-07-102798.

8. Lin RJ, Nagy L, Inoue S, Shao W, Miller WH Jr, Evans RM. Role of the histone deacetylase complex in acute promyelocytic leukaemia. Nature. 1998; 391:811–4. https://doi.org/10.1038/35895.

9. Grignani F, De Matteis S, Nervi C, Tomassoni L, Gelmetti V, Cioce M, Fanelli M, Ruthardt M, Ferrara FF, Zamir I, Seiser C, Grignani F, Lazar MA, et al. Fusion proteins of the retinoic acid receptor-alpha recruit histone deacetylase in promyelocytic leukaemia. Nature. 1998; 391:815–8. https://doi.org/10.1038/35901.

10. Di Croce L, Raker VA, Corsaro M, Fazi F, Fanelli M, Faretta M, Fuks F, Lo Coco F, Kouzarides T, Nervi C, Minucci S, Pelicci PG. Methyltransferase recruitment and DNA hypermethylation of target promoters by an oncogenic transcription factor. Science. 2002; 295:1079–82. https://doi.org/10.1126/science.1065173.

11. Carbone R, Botrugno OA, Ronzoni S, Insinga A, Di Croce L, Pelicci PG, Minucci S. Recruitment of the histone methyltransferase SUV39H1 and its role in the oncogenic properties of the leukemia-associated PML-retinoic acid receptor fusion protein. Mol Cell Biol. 2006; 26:1288–96. https://doi.org/10.1128/MCB.26.4.1288-1296.2006.

12. Hoemme C, Peerzada A, Behre G, Wang Y, McClelland M, Nieselt K, Zschunke M, Disselhoff C, Agrawal S, Isken F, Tidow N, Berdel WE, Serve H, et al. Chromatin modifications induced by PML-RARalpha repress critical targets in leukemogenesis as analyzed by ChIP-Chip. Blood. 2008; 111:2887–95. https://doi.org/10.1182/blood-2007-03-079921.

13. Martens JH, Brinkman AB, Simmer F, Francoijs KJ, Nebbioso A, Ferrara F, Altucci L, Stunnenberg HG. PML-RARalpha/RXR Alters the Epigenetic Landscape in Acute Promyelocytic Leukemia. Cancer Cell. 2010; 17:173–85. https://doi.org/10.1016/j.ccr.2009.12.042.

14. Singh AA, Mandoli A, Prange KH, Laakso M, Martens JH. AML associated oncofusion proteins PML-RARA, AML1-ETO and CBFB-MYH11 target RUNX/ETS-factor binding sites to modulate H3ac levels and drive leukemogenesis. Oncotarget. 2017; 8:12855–65. https://doi.org/10.18632/oncotarget.14150.

15. Breccia M, Lo-Coco F. Arsenic trioxide for management of acute promyelocytic leukemia: current evidence on its role in front-line therapy and recurrent disease. Expert Opin Pharmacother. 2012; 13:1031–43. https://doi.org/10.1517/14656566.2012.677436.

16. Lo-Coco F, Orlando SM, Platzbecker U. Treatment of acute promyelocytic leukemia. N Engl J Med. 2013; 369: 1472. https://doi.org/10.1056/NEJMc1309895.

17. Zhu J, Gianni M, Kopf E, Honore N, Chelbi-Alix M, Koken M, Quignon F, Rochette-Egly C, de The H. Retinoic acid induces proteasome-dependent degradation of retinoic acid receptor alpha (RARalpha) and oncogenic RARalpha fusion proteins. Proc Natl Acad Sci U S A. 1999; 96:14807–12.

18. Nervi C, Ferrara FF, Fanelli M, Rippo MR, Tomassini B, Ferrucci PF, Ruthardt M, Gelmetti V, Gambacorti-Passerini C, Diverio D, Grignani F, Pelicci PG, Testi R. Caspases mediate retinoic acid-induced degradation of the acute promyelocytic leukemia PML/RARalpha fusion protein. Blood. 1998; 92:2244–51.

19. Jeanne M, Lallemand-Breitenbach V, Ferhi O, Koken M, Le Bras M, Duffort S, Peres L, Berthier C, Soilihi H, Raught B, de The H. PML/RARA oxidation and arsenic binding initiate the antileukemia response of As2O3. Cancer Cell. 2010; 18:88–98. https://doi.org/10.1016/j.ccr.2010.06.003.

20. Norsworthy KJ, Altman JK. Optimal treatment strategies for high-risk acute promyelocytic leukemia. Curr Opin Hematol. 2016; 23:127–36. https://doi.org/10.1097/MOH.0000000000000215.

21. Raelson JV, Nervi C, Rosenauer A, Benedetti L, Monczak Y, Pearson M, Pelicci PG, Miller WH Jr. The PML/RAR alpha oncoprotein is a direct molecular target of retinoic acid in acute promyelocytic leukemia cells. Blood. 1996; 88:2826–32.

22. Rosenauer A, Raelson JV, Nervi C, Eydoux P, DeBlasio A, Miller WH Jr. Alterations in expression, binding to ligand and DNA, and transcriptional activity of rearranged and wild-type retinoid receptors in retinoid-resistant acute promyelocytic leukemia cell lines. Blood. 1996; 88:2671–82.

23. Marasca R, Zucchini P, Galimberti S, Leonardi G, Vaccari P, Donelli A, Luppi M, Petrini M, Torelli G. Missense mutations in the PML/RARalpha ligand binding domain in ATRA-resistant As(2)O(3) sensitive relapsed acute promyelocytic leukemia. Haematologica. 1999; 84:963–8.

24. Cote S, Zhou D, Bianchini A, Nervi C, Gallagher RE, Miller WH Jr. Altered ligand binding and transcriptional regulation by mutations in the PML/RARalpha ligand-binding domain arising in retinoic acid-resistant patients with acute promyelocytic leukemia. Blood. 2000; 96:3200–8.

25. Adams D, Altucci L, Antonarakis SE, Ballesteros J, Beck S, Bird A, Bock C, Boehm B, Campo E, Caricasole A, Dahl F, Dermitzakis ET, Enver T, et al. BLUEPRINT to decode the epigenetic signature written in blood. Nat Biotechnol. 2012; 30:224–6. https://doi.org/10.1038/nbt.2153.

26. Cancer Genome Atlas Research N, Ley TJ, Miller C, Ding L, Raphael BJ, Mungall AJ, Robertson A, Hoadley K, Triche TJ Jr, Laird PW, Baty JD, Fulton LL, Fulton R, et al. Genomic and epigenomic landscapes of adult de novo acute myeloid leukemia. N Engl J Med. 2013; 368:2059–74. https://doi.org/10.1056/NEJMoa1301689.

27. Schoofs T, Rohde C, Hebestreit K, Klein HU, Gollner S, Schulze I, Lerdrup M, Dietrich N, Agrawal-Singh S, Witten A, Stoll M, Lengfelder E, Hofmann WK, et al. DNA methylation changes are a late event in acute promyelocytic leukemia and coincide with loss of transcription factor binding. Blood. 2013; 121:178–87. https://doi.org/10.1182/blood-2012-08-448860.

28. Nebbioso A, Clarke N, Voltz E, Germain E, Ambrosino C, Bontempo P, Alvarez R, Schiavone EM, Ferrara F, Bresciani F, Weisz A, de Lera AR, Gronemeyer H, et al. Tumor-selective action of HDAC inhibitors involves TRAIL induction in acute myeloid leukemia cells. Nat Med. 2005; 11:77–84.

29. Martens JH, Stunnenberg HG. The molecular signature of oncofusion proteins in acute myeloid leukemia. FEBS Lett. 2010; 584:2662–9. https://doi.org/10.1016/j.febslet.2010.04.002.

30. Mandoli A, Singh AA, Prange KHM, Tijchon E, Oerlemans M, Dirks R, Ter Huurne M, Wierenga ATJ, Janssen-Megens EM, Berentsen K, Sharifi N, Kim B, Matarese F, et al. The Hematopoietic Transcription Factors RUNX1 and ERG Prevent AML1-ETO Oncogene Overexpression and Onset of the Apoptosis Program in t(8;21) AMLs. Cell Rep. 2016; 17:2087–100. https://doi.org/10.1016/j.celrep.2016.08.082.

31. Saeed S, Logie C, Francoijs KJ, Frigè G, Romanenghi M, Nielsen FG, Raats L, Shahhoseini M, Huynen M, Altucci L, Minucci S, Martens JH, Stunnenberg HG. Chromatin accessibility, p300, and histone acetylation define PML-RARα and AML1-ETO binding sites in acute myeloid leukemia. Blood. 2012; 120:3058–68. https://doi.org/10.1182/blood-2011-10-386086.

32. Valiulienė G, Stirblytė I, Jasnauskaitė M, Borutinskaitė V, Navakauskienė R. Anti-leukemic effects of HDACi Belinostat and HMTi 3-Deazaneplanocin A on human acute promyelocytic leukemia cells. Eur J Pharmacol. 2017; 799:143–153.

33. De Bellis F, Carafa V, Conte M, Rotili D, Petraglia F, Matarese F, Françoijs KJ, Ablain J, Valente S, Castellano R, Goubard A, Collette Y, Mandoli A, et al. Context-selective death of acute myeloid leukemia cells triggered by the novel hybrid retinoid-HDAC inhibitor MC2392. Cancer Res. 2014; 74:2328–39. https://doi.org/10.1158/0008-5472.

34. Amin HM, Saeed S, Alkan S. Histone deacetylase inhibitors induce caspase-dependent apoptosis and downregulation of daxx in acute promyelocytic leukaemia with t(15;17). Br J Haematol. 2001; 115:287–97.

35. Guidez F, Ivins S, Zhu J, Soderstrom M, Waxman S, Zelent A. Reduced retinoic acid-sensitivities of nuclear receptor corepressor binding to PML- and PLZF-RARalpha underlie molecular pathogenesis and treatment of acute promyelocytic leukemia. Blood. 1998; 91:2634–42.

36. He LZ, Guidez F, Tribioli C, Peruzzi D, Ruthardt M, Zelent A, Pandolfi PP. Distinct interactions of PML-RARalpha and PLZF-RARalpha with co-repressors determine differential responses to RA in APL. Nat Genet. 1998; 18:126–35. https://doi.org/10.1038/ng0298-126.

37. Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics 2010; 26:873–881.

38. Kanitz A, Gypas F, Gruber AJ, Gruber AR, Martin G, Zavolan M. Comparative assessment of methods for the computational inference of transcript isoform abundance from RNA-seq data. Genome Biol 2015; 16:150.

39. Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz LA Jr, Kinzler KW. Cancer genome landscapes. Science 2013; 339:1546–s1558.