Epigenetic changes in human model KMT2A leukemias highlight early events during leukemogenesis.

Chromosomal translocations involving KMT2A gene are one of the most common genetic alterations found in pediatric acute myeloid leukemias (AML) although the molecular mechanisms that initiate the disease remain incompletely defined. To elucidate these initiating events we have used a human model system of AML driven by the KMT2A-MLLT3 (KM3) fusion. More specifically, we investigated changes in DNA methylation, histone modifications, and chromatin accessibility at each stage of our model system and correlated these with expression changes. We observe the development of a profound hypomethylation phenotype in the early stages of leukemic transformation after KM3 addition along with loss of expression of stem cell associated genes along with skewed expression in other genes such as S100A8/9 implicated in leukemogenesis. In addition, early increases in the expression of the lysine demethylase KDM4B was functionally linked to these expression changes as well as other key transcription factors. Remarkably, our ATAC-seq data showed that there were relatively few leukemiaspecific changes and the vast majority corresponded to open chromatin regions and transcription factor clusters previously observed in other cell types. Integration of the gene expression and epigenetic changes revealed the adenylate cyclase gene ADCY9 as an essential gene in KM3-AML, and suggest the potential for autocrine signalling through the chemokine receptor CCR1 and CCL23 ligand. Together, our results suggest that KM3 induces subtle changes in the epigenome while co-opting the normal transcriptional machinery to drive leukemogenesis.

Introduction

Leukemias are a class of blood disorders characterized by the uncontrolled proliferation of hematopoietic stem and progenitor cells (HSPC) that have acquired a block in the normal process of differentiation.[1-3] Recent large-scale studies have highlighted the age-related genetic differences between pediatric and adult acute myeloid leukemia (AML) and have also revealed the high level of genetic heterogeneity in this disease.[4] Not only are chromosomal translocations more frequent in younger patients, but recurrent mutations in specific genes implicated in the disease show a bias in either pediatric patients (e.g., NRAS, KIT, KRAS) or adults (DNMT3A, NPM1, IDH1/2).[5] As a result, defining the combined role of these mutations in each patient, especially in pediatric patients in whom predisposition variants may play a greater role, is extremely challenging and increases the complexity of identifying novel potential targeted therapeutics. Studies underscore the fact that even though the mutational burden in AML is one of the lowest of any cancer[4] there is still tremendous heterogeneity between leukemic patients which creates challenges when trying to define the genetic determinates of the disease.

As noted, chromosomal fusions are common in pediatric AML, particularly those involving the lysine specific methyl transferase 2A gene (KMT2A; also known as the mixed lineage leukemia [MLL] gene). KMT2A is rearranged in approximately 10% of all leukemias[6] but the frequency in infant acute lymphocytic leukemia (ALL, >70%[7]) and AML (>35%[7]) is much higher. The KMT2A gene encodes a large (500 kDa) and complex 38-exon protein[8,9] that can be cleaved by taspase-1 into two separate fragments (KMT2A-C and KMT2A-N).[10] The conserved SET domain located at the C-terminus of the protein is responsible for the methylation of the histone H3 at lysine 4 (H3K4). Previous studies have shown that KMT2A is part of a large macromolecular complex composed of different proteins[11] that function to improve the stability of KMT2A allowing the complex to regulate the transcriptional activation of HOX genes. The KMT2A gene can be fused in-frame with more than 120 different partners, creating fusion proteins that are typically associated with poor prognosis leukemias.[6] One of the most common fusion partners is the MLLT3 gene (also known as AF9), which is found in 30% and 13% of KMT2Arearranged AML and ALL,[12] respectively; this fusion is associated with an intermediate risk. The KMT2AMLLT3 (KM3) fusion protein is part of the DOT1L complex (DOTCOM) and leads to aberrant expression of specific target genes marked by H3K79 methylation.[13] Despite the identification of several key gene targets required for transformation, the complete molecular mechanisms utilized by this oncogenic fusion are still unclear.

The strong oncogenic potential of KMT2A fusions makes them ideal drivers for in vivo experimental model systems to explore molecular mechanisms involved in leukemogenesis.[14] To specifically overcome the problems of patients’ genetic heterogeneity and scarcity of samples, we have developed a single-donor, human model leukemia system using healthy CD34+ cord blood HSPC.[15] This adaptation of a prior human model[16] uses a retrovirus to deliver the human KM3 fusion gene which we have shown is sufficient to generate a human leukemia.[15] Transduced cells are cultured in vitro for 30-40 days prior to injection into immune-deficient mice which develop a leukemia 24-30 weeks later. Because we can sequence the transcriptome/genome/exome of the initial CD34+ cells used to generate the model AML, the genetic background (e.g., single nucleotide variants) of the initial donor is well defined, allowing the potential role of acquired mutations to be assessed precisely. In addition to these advantages, this model also allows the study of genetic mechanisms involved in the initiation of the disease, when the fusion is first introduced but before the complete transformation of the cells. Our previous genomic analyses demonstrated that the oncogenic fusion alone is sufficient for the development of leukemias, without the requirement for any recurrent secondary mutations,[15] despite the presence of these in pediatric patients with KMT2A-mutated AML. We have also used these data to identify a number of genes that represent novel biomarkers in patients with KMT2A fusions.[17] Interestingly, while approximately onethird of these biomarker genes were expressed in our model leukemia system shortly after cells were transduced with the KM3 fusion, the majority were only expressed after xenotransplantation. This observation, coupled with the limited potential for in vitro growth of these cells, suggests that cells partially transformed (or “primed”) by KM3 may require additional in vivo signals (e.g., from the bone marrow niche) to complete their leukemic transformation. The single-donor model leukemias therefore not only recapitulate the behavior and phenotype of the disease, but also provide an experimental system to study genetic and epigenetic mechanisms involved in the disease, and a unique insight into early stages of transformation.

In the present study, we leverage these advantages to perform a detailed epigenetic analysis of the various stages of the model and, through correlation with model and expression data from patients with leukemia, define the epigenetic changes driven by the oncogenic fusion that contribute to leukemia development. Analysis of changes in the patterns of DNA methylation confirmed previous observations of profound hypomethylation in KMT2A-translocated AML. Interestingly, however, B-cell ALL driven by the same fusion show a smaller number of differentially methylated cytosines relative to HSPC methylation, which are predominately hypermethylated. Despite the generally poor global correlation between DNA methylation and gene expression, we identified ADCY9, a member of the adenylyl cyclase family, as a gene essential for KM3-AML growth that also exhibits coordinated changes in DNA methylation and gene expression. We further characterized changes in chromatin accessibility through an assay for transposaseaccessible chromatin with high-throughput sequencing (ATAC-sequencing), along with alterations in histone modifications including H3K4me3 and H3K79me2. Interestingly, our ATAC-sequencing analysis revealed that the vast majority of regions of open chromatin are shared with normal HSPC and monocytes, with very few being specific to the leukemic cells. Even those regions that are unique to the leukemia show a high degree of overlap with regions of open chromatin seen in other cell types and are also characterized by the presence of known clusters of transcription factor binding sites. Integration of the epigenetic data with expression data highlighted a role for the histone demethylase KDM4B during the initial stages of leukemic transformation through its impact on S100A8/S100A9 expression levels, whose relative ratios have been demonstrated to be critical for blocking differentiation of leukemic cells.[18] Collectively, these results suggest that the leukemic transformation of normal HSPC by KM3 involves a very subtle epigenetic shift that also implicates co-option of normal transcriptional networks.

Methods

Patients’ samples and model leukemia generation

All pediatric and adult AML patients’ samples and additional clinical information used in this study were collected by the Banque de Cellules Leucémiques du Québec (BCLQ) with informed consent and project approval from the Research Ethics Boards of CHU Ste-Justine, Maisonneuve-Rosemont Hospital and Université de Montréal. RNA was extracted from samples preserved in TRIzol (Invitrogen) according to standard protocols, while DNA for methylation-sequencing (methyl-sequencing) analysis was obtained from viably conserved samples. The characteristics of the patients’ AML samples are shown in Online Supplementary Table S1. The model leukemias used in this study were generated as previously described.[15]

Human CD34+ hematopoietic stem and progenitor cell isolation and cell culture

CD34+ human HSPC were isolated from fresh umbilical cord blood, collected under informed consent at Sainte-Justine pediatric hospital (distributed by HémaQuébec’s Public Cord Blood Bank [Montreal, Quebec, Canada]) or CHU de Quebec and Hotel Dieu de Levis. CD34+ cord blood cells were positively selected as previously described[15] and processed immediately without in vitro culture. Briefly, cord blood was diluted with phosphate-buffered saline (PBS) + citrate 0.6% (1:2) and mixed with Ficoll. After spinning the mixture at 2,100 rpm, for 30 min at room temperature with the lowest deceleration speed, white blood cells found at the interphase were collected and diluted with PBS (1:2). Cells were then spun at 2,100 rpm for 10 min and resuspended in PBS + citrate 0.6% + bovine serum albumin 0.5%. An isolation kit (Miltenyi Biotec; #130-100-453) composed of beads recognizing CD34+ cells was used according to the manufacturer’s protocol. The Posseld method on an AutoMACS® machine was then used to isolate CD34+ cells and the final yield was quantified using a hemocytometer. Leukemia cells lines (KG1a, THP-1, NOMO1 and MOLM-13) were grown in RPMI medium supplemented with 10% heat-inactivated fetal bovine serum, penicillin (100 U/mL) and streptomycin (100 mg/mL) at 37°C in 5% CO2. Cell lines were transfected using a lentiviral protocol as previously described.[15] The shRNA sequences used are shown in Table 1.

Table 1.

List of primers used for shRNA-mediated knockdown of ADCY9 expression.

Chromatin immunoprecipitation, RNA-sequencing and ATAC-sequencing

Cells were cross-linked using formaldehyde 1% (for 7 min, room temperature) before glycine quenching (0.125 nM, 5 min, room temperature), cell lysis and sonication (Covaris S2 sonicator). Protein G-coupled beads (Protein G Dynabeads™) were incubated overnight on a rotator at 4°C with 3 mg of antibodies against different histone marks (H3: #ab1791; H3K4me3: #ab8580-100; H3K79me2: #ab3594). Sonicated chromatin was incubated with 250 mL of beads for 4 h on a rotator at 4°C before sequential washing with: ‘Low Salt’ buffer (0.5% NP40, 15 mM KCl, 10 mM Tris pH 8.0, 1 mM EDTA), three different ‘High Salt’ buffers (0.5% Triton, 10 mM Tris pH 8.0, 100 mM [2] or 400 mM [3] or 500 mM [4] NaCl), twice in LiCl buffer (0.5% NP40, 250 mM LiCl, 10 mM Tris pH 8.0, 1 mM EDTA). ChIPed material was then eluted (buffer: 10 mM Tris pH 8.0, 1 mM EDTA) and incubated overnight at 65°C with shaking (1,200 rpm) to reversecrosslink. An RNaseA was added (for 1 h at 37°C) and Proteinase K digestion was performed (1 h at 37°C) before DNA was purified by phenol-chloroform. A TruSeq ChIP Library Preparation Kit (Illumina, #IP-202-9001) was used to construct DNA libraries as described in the manufacturer’s instructions, with polymerase chain reaction (PCR) for a total of 12 cycles using the Illumina indexed library primers. An Illumina NextSeq 500 was used to sequence samples (75 bp paired-end) according to standard Illumina protocols. For all RNA-sequencing experiments, RNA was purified from cells (1x106) in TRIzol reagent, purified on RNeasy columns (Qiagen) and then used as input for Illumina stranded Tru-seq protocols. For CD34+ cells transduced with the KM3 fusion and cultured in vitro, RNA was collected after 30 days from GFP+ cells while for all shRNA knockdown experiments RNA was collected 4 days after lentiviral transduction. Lastly, RNA from the xenografted model AML was collected from the bone marrow of mice sacrificed 24-30 weeks after injection. All sequencing in this study was performed by the genomics core facility at the Institute for Research in Immunology and Cancer. The ATAC-sequencing protocol used was based on a method published by Buenrostro et al.[19] and was performed using 50,000 cells. Samples were sequenced (75 bp paired-end) according to Illumina protocols. The primers used are shown in Table 2.

Table 2.

List of primers used to prepare assays for transposase-accessible chromatin (ATAC)-sequencing libraries.

DNA methylation capture sequencing

Genomic DNA from cell pellets from either KM3 patients’ samples, pooled cord blood donors (freshly isolated), cultured cells (collected after 30 days) or model AML cells (collected 24 weeks after injection) was isolated using the PureLink Genomic DNA Mini Kit according to the manufacturer’s protocol (Life Technologies, cat. #K1820-01). Three micrograms of DNA from each sample were fragmented to 140-180 bp DNA fragments using Covaris S2 (parameters: 10% duty cycle at intensity 5 for 6 cycles of 60 s with 200 cycles set at sweeping mode). Methylsequencing libraries for DNA methylation analysis were prepared using the SureSelectXT Human Methyl-Seq Target enrichment system (Agilent Technologies, cat. #5190-4661) according to the manufacturer’s instructions. Briefly, DNA was fragmented, purified and used for end-repair and adenylated prior to ligation of methylation adaptors (synthesized with 5’-methylcytosine instead of cytosine). The adaptor-ligated DNA samples were purified using AMPure XP beads (Beckman Coulter Genomics, cat. #A63880). Nine hundred nanograms of each DNA sample were hybridized with Agilent SureSelect Methyl-Seq biotinyled RNA baits for 24 h at 65°C. The libraries were captured using Dynal MyOne Streptavidin T1 magnetic beads (Invitrogen, cat. #65601). The captured DNA was subjected to bisulfite conversion with the EZ DNA Methylation Gold Kit (Zymo Research, cat. #D5005) as described in the manufacturer’s instructions, followed by PCR for a total of 16 cycles using the Illumina indexed library primers. The final libraries were quantified by 2100 Bioanalyzer chips (Agilent Technologies) with quality control on a MiSeq (Illumina) sequencer. Libraries were sequenced on multiple lanes of a HiSeq2000 flow cell at the IRIC genomics platform and the FastQC R package was used to assess the quality of the reads.

Bioinformatics analysis of ChIP-sequencing and ATACsequencing

The quality of the raw sequence data was checked by FastQC (version 0.11.4). Trimmomatic (version 0.32) was used to remove adapters and keep reads with quality scores >30. Paired-end sequencing reads were then mapped to the human reference genome (GRCH37/hg19) using BWA (version 0.7.10) with the default parameters. PCR duplicates were removed by Samtools (version 1.9). For ATAC-sequencing data, reads from mitochondrial DNA were removed (samtools idxstats input.bam | cut -f 1 | grep -v MT | xargs samtools view -b input.bam > output.bam). The EaSeq software package (version 1.04)[20] was used to analyze data and integrate gene expression and data from epigenetic marks from chromatin immunoprecipitation (ChIP)-sequencing experiments. To calculate the enrichment of each histone mark, the analysis focused on specific regions in the genome with a distance from the transcription start site of -1,000 to +3,000 bp for H3K79me2, and -2,000 to +2,000 bp for H3K4me3. MACS2 (version 2.2.4)[21] was the peak calling algorithm used to determine the list of significant enrichment peaks, following the default para - meters (callpeak –f BAMPE –g hs –q 0.05). Intervene[22] was used to find peaks in common in our datasets (bedtools intersect –a – b –f 0,50 –r) and generate Venn diagrams (Pybedtools). A ngs.plot[23] was used to visualize global repartition of mapped reads across the genome. Previously published data for lineagespecific ATAC-sequencing in a range of normal blood cells[24] were accessed through the UCSC track hub (https://s3-us-west-1.amazonaws. com/chang-public-data/2016_NatGen_ATAC-AML/hub.txt) DNAse hypersensitive sites and transcription factor clusters generated from the ENCODE consortium V3 were analyzed by custom Perl Scripts. Public ATAC-sequencing peaks for all cell types were identified using a threshold of a normalized score >60 for all tracks for starts/ends of peaks. Fastq files for BLUEPRINT ATAC-sequencing data (accessions EGAD00001002710 and EGAD00001002709) were downloaded, remapped and analyzed with MACS2 as described above. Gene set enrichment analysis was performed using software (version 3.0) from the Broad Institute[25] and a pre-ranked list based on significantly differentially expressed genes between CD34+ and CD34+KM3 cells using the R package DESeq2.[26]

Mapping methyl-sequencing reads and methylation analysis

Paired end sequencing reads (2x100 bp) were mapped to the human reference genome (GRCH37/hg19) using Bismark with a directional alignment mode and Bowtie2 as the alignment process. All putative duplicate reads were removed by previously described Perl scripts.[27] The methylation calls were extracted with Bismark-methylation-extractor and the analysis was performed using custom python scripts, R package MethylKit 0.9.2 and SeqMonk 1.45.2 as the visualization tool (). Hierarchical clustering of the samples based on the differentially methylated cytosines was performed using the computed Pearson correlation distance between samples. For each cell type (CD34+, CD34+KM3, KM3- AML model) the DNA methylation data corresponded to the combination of methyl-sequencing libraries from two biological replicates. Analysis of methylated cytosines was restricted to positions covered by a minimum of eight reads coverage for each biological replicate. To identify differentially methylated cytosines (DMC) a logistic regression model was applied to compare the fraction of methylated cytosines across the three types of cells and a sliding linear model (SLIM) method was used to correct P-values for multiple testing and derive q-values.[28] To define the DMC, we specified two parameters: q-values below 0.01 and a threshold of percentage methylation difference ≤40%. All functional annotations for genomic regions were defined based on GRCh37/hg19 assembly downloaded from the UCSC website. Gene bodies were defined as the transcribed regions from the start to the end of transcription sites for each RefSeq entry and promoter regions were defined as 1 kb of sequence upstream of the transcription start site for each RefSeq transcript.

Results

Transcriptional consequences of KM3 expression

To characterize the short-term transcriptional changes induced by the KM3 fusion, we analyzed the RNAsequencing data from the initial unmanipulated cord blood cells and these cells after their transduction with the KM3 fusion. As previously described,[16] untransduced cells are largely lost in culture after about 1 month of in vitro culture and the remaining cells show induction of genes associated with a myeloid phenotype (e.g., CCR1, CD14) (Figure 1A). Interestingly, transduced cells also display the loss of normal or leukemia stem cell markers (e.g., CD34 and GPR56, respectively) and the loss of expression of other genes previously associated with AML induction in other subtypes (e.g., MN1[30]). This observation is supported by expression data from patients with KMT2A-translocated AML in the TARGET study,[4] which showed an equivalent loss of MN1 and GPR56[4] (Figure 1B) and a similar loss of expression in studies of GPR56 seen in adult AML.[31] A global view of the expression changes through gene set enrichment analysis confirmed the loss of stem cell phenotype and the gain in expression of HOXA9-MEIS1 targets consistent with leukemic transformation (Online Supplementary Figure S1).

Figure 1.

Gene expression changes in the (A) Volcano plot of differentially expressed genes (red dots are genes with a P value <10-5, and a log2 fold change >3) from RNA-sequencing data of CD34+ cells and CD34+KM3 cells. Highlighted genes include known markers for stem cells, myeloid cells or previously reported biomarkers for KM3-AML. (B) Scatterplot of MN1 and GPR56 expression in approximately 1,000 cases of pediatric AML in the TARGET cohort. Red dots represent patients with KMT2A-rearranged AML. (C) Summary of transcription factor binding motifs found to be enriched in promoters of upregulated genes in (A) with either ShinyGO v0.6134 or MAGICTRICKS.[33] (D) Gene expression levels of transcription factor families highlighted in (C) in either CD34+ cells (blue) or CD34+KM3 cells (red). (E) High levels of expression of TLR4, S100A8 and S100A9 are seen after introduction of the KMT2A-MLLT3 fusion and are maintained in the model AML. (F) Loss of KDM4B expression through shRNA knockdown in the KMT2A-MLLT3+ cell line (THP-1) induces growth inhibition (*P<0.005, **P<0.001) and (G) reduced expression of S100A8/A9, CEBPb and other cell-cycle related genes. (H) shRNA knockdown of KDM4B in CD34+ cells transduced with the KM3 fusion show gene expression decreases very similar to those in THP-1 cells. FPKM: fragments per kilobase of exon per million mapped reads

To examine the regulation underlying the transcriptional changes seen, we studied the transcription factor binding motifs upstream of the differentially expressed genes. Analysis of enriched motifs using several different software (including HOMER,[32] MAGICTRICKS,[33] and ShinyGO[34]) produced differing lists of potential transcription factors (Figure 1C), but factors in the CEBP and FOS families, along with NFE2, were highlighted consistently, at differing ranks. Analysis of RNA-sequencing data confirmed that the expression of all members of the CEBP and FOS families increased during in vitro culture (Figure 1D), along with transient changes in ETS, MAF, and JUN families (Online Supplementary Figure S2). We also observed an increase in the expression of S100A9 and S100A8 (Figure 1E), which have been described as direct targets of CEBPb.[35] It has been recently demonstrated that the ratio of these calcium-binding proteins is critical for blocking AML differentiation,[18] with S100A8 opposing the differentiation promoted by S100A9, which acts through TLR4. Interestingly, and in agreement with these findings, while the expression of both S100A8 and S100A9 increased when KM3 was added to CD34+ cells, there was a larger increase in S100A8, changing the ratio of S100A8:S100A9 from 1.1:1 to 3.2:1, favoring inhibition of differentiation. Other CEBP members, including CEBPb, have also previously been implicated in the process of normal myeloid differentiation,[36] which suggests that transcription factors normally present are coopted to activate genes that are part of a leukemic transcriptional program. At the same time, these factors and others upregulated are also involved in immune and inflammatory pathways, suggesting that there may be additional transcriptional programs that are activated and superimposed on the myeloid differentiation that we observe. To determine whether changes in the expression of chromatin-modifying enzymes might be correlated with our RNA-sequencing data, we looked for enzymes whose expression increased at least 2-fold during in vitro culture and whose expression was maintained in both the model and patients’ AML (Online Supplementary Table S2). This analysis highlighted a single gene with this profile, KDM4B, which has been implicated in murine AML, albeit with functional redundancy.[37] The shRNA-mediated knockdown of KDM4B in THP-1 cells showed a significant impact on growth (Figure 1F) and RNA-sequencing from these cells revealed a profound downregulation of S100A8 and S100A9, CEBPb, P2RY2 (a KM3-AML biomarker[17]), and other cell cycle regulators (Figure 1G). We next performed the same experiment but with CD34+ cells transduced with the KM3 fusion to validate the changes in our model system. The data showed very similar gene expression changes in these genes with the exception of CEBPb, although other CEBP family members were downregulated (Figure 1H). Together, these data support the hypothesis that the early epigenetic changes, at least partially mediated through KDM4B, lead to activation of S100A8/9 as well as other genes required for sustained proliferation. Given the similarity of their DNA binding consensus sequences, defining the precise contribution of different CEBP family members in this process will require further focused studies.

Sequential changes in DNA methylation are not globally correlated with expression changes

Given the known role of KM3 as an epigenetic regulator and the demonstrated importance of DNA methylation changes in other subtypes of AML,[38,39] we next sought to investigate changes in DNA methylation in our model leukemias and genetically matched patients’ samples. Using methyl-sequencing, capture-based, bisulfite sequencing,[40] we assessed DNA methylation levels genome-wide at 84 Mbp of target sequence regions including CpG islands, promoters and published differentially methylated regions. Biological duplicates at each step of our model leukemia system were compared to bisulfite-treated DNA from three adult KM3-AML patients. In line with expression data from our model leukemias,[15] the methylation data between replicates were highly consistent (Online Supplementary Figure S3) and were also consistent with published data regarding global methylation levels (Online Supplementary Figure S4).[41]

Hierarchal clustering of DMC (CpG coverage ≥8, qvalue ≤0.1 and methylation differences ≥40%) demonstrated that relatively small numbers of changes occurred after introduction of the KM3 fusion gene and in vitro culture (Figure 2A). The number of genes with DMC present varied from 7,645 (CD34+ vs. CD34+KM3) to 13,947 (CD34+ vs. model AML); however a large percentage of these genes (40% and 20%, respectively) contained only a single DMC. Interestingly, relative to the initial CD34+ cells with the KM3 fusion, there was a striking increase (>6 fold) in the number of DMC detected in the final model leukemias and patients’ samples. This suggests that the initial epigenetic changes induced by the fusion gene are relatively limited, with the majority of changes occurring only during in vivo growth. Along with the previous stepwise changes in gene expression patterns we have reported,[15] these results again support the importance of microenvironmental cues received in vivo. With respect to the direction of change, in agreement with recent observations in patients with KMT2A-AML,[41] DMC observed in both model and patient KM3-AML were profoundly biased (~89%) towards demethylation (Figure 2B). In opposition to this trend, model leukemias generated from the same donor cells but which developed a lymphoid phenotype (B-cell ALL) had primarily hyper-methylated cytosines, and fewer DMC compared to the AML samples. The majority (>75%) of DMC were located in gene bodies and particularly in promoter regions and putative functional elements at various stages of the model AML and in patients (Online Supplementary Figures S5 and S6).

Figure 2.

DNA methylation changes in model leukemias. (A) Hierarchical clustering of differentially methylated cytosines (DMC) identified by methyl-sequencing analysis of model leukemias at different stages and KM3 patients’ samples (pAML). Numbers and arrows represent the number of DMC relative to the initial CD34+ cells. (B) The KM3 fusion induces a strong hypomethylation phenotype in both patient and model acute myeloid leukemias (pAML and mAML, respectively), whereas model B-cell acute lymphoblastic leukemias (mB-ALL) predominately exhibit hypermethylation. (C) As noted in other global studies, the overall correlation between DNA methylation and expression is relatively limited, likely due in part to the complexity of the localization and size of the changes. (D) Differentially methylated regions around genes known to be suppressed (GATA2) or activated (HOXA9) in AML can be seen comparing CD34+ cells versus AML patients’ cells. (E, F) ADCY9 exhibits hypomethylation in its promoter region specifically in the model AML, compared to CD34+ cells, which is correlated with gene expression levels, specifically in KMT2A-MLLT3 AML models as well as patients (E) and external data from TARGET (F). (G) Knockdown of ADCY9 using multiple shRNA shows a significant impact specifically on KM3-AML cell proliferation but not on AML cells lacking the fusion (KG1a; *P<0.005, **P<0.001).

To clarify the importance of the localization of these changes and their potential impact, we next examined the correlation between DNA methylation changes and gene expression. Despite the general view regarding the impact of DNA methylation changes in the promoter regions and transcription levels, we observed only a weak global correlation between the two (Figure 2C). While we could find examples of expected changes (e.g., in previously described KMT2A target genes such as GATA2, HOXA7 and HOXA9[42,43]) (Figure 2D), overall, there was a very poor global correlation. This observation mirrors those made in recent comprehensive single-cell studies of DNA methylation changes during normal hematopoiesis in humans[44] as well as previous bulk studies.[45] We were able to identify new examples of genes with evidence for epigenetic regulation such as the atypical adenylate cyclase gene ADCY9. The promoter region of ADCY9 exhibits hypomethylation specifically in the model AML compared to normal CD34+ cells and model B-cell ALL (Figure 2E) and the gene is only expressed in the KMT2Arearranged model and patients’ leukemias (Figure 2F). ADCY9 is a poorly characterized and distantly related member of the adenylate cyclase gene family and no specific inhibitors have yet been described, although it was recently identified as a fusion partner of the KMT2D gene in a case of ALL.[46] We therefore performed shRNA knockdowns in several KM3-AML cell lines to determine whether ADCY9 expression was relevant to these leukemias. We found that all the KM3-AML lines showed impaired proliferation with multiple shRNA in contrast to no effect in a non-KM3-AML control line (Figure 2G), demonstrating that ADCY9 is required in these cells.

Beyond regulating specific genes, we next wanted to examine whether global methylation changes might have an impact on gene expression through their localization within the DNA binding sites of transcription factors. To do this, we selected all the DMC with coverage that exhibited consistently higher (≥40%) methylation levels in CD34+ cells than in either CD34+KM3 or model leukemias. Transcription factors in the ENCODE ChIPsequencing dataset (covering 338 factors in 130 cell lines) that overlapped the positions of these cytosines were then used to generate a network (using GOnet[47]) based on shared gene ontology terms (Online Supplementary Figure S7). This network contains a large collection of factors that: (i) are highly expressed in normal bone marrow, (ii) exhibit a large loss of DNA methylation directly in their binding sites after KM3 addition, and (iii) are enriched for activity in regulating hematopoietic and myeloid differentiation. The results indicate that despite the poor overall correlation between DNA methylation and the expression of individual genes, the changes that occur are globally relevant for leukemic transformation.

Epigenetic analysis reveals potential CCR1-CCL23 autocrine signaling in KM3-acute myeloid leukemia

Given that the KM3 fusion is a well-known epigenetic regulator, we next asked whether the presence of the gene fusion had a specific impact on the chromatin organization in our model leukemias. To investigate this, we performed ChIP-sequencing on relevant histone marks (H3K4me3, H3K79me211) and ATAC-sequencing on samples of each step of our leukemia model (CD34+, CD34+KM3, KM3-AML). Enrichment of the different histone marks correlated with gene expression levels (Online Supplementary Figure S7) and expected spatial distribution, and the majority of peaks seen in the model leukemias were common to both the initial and transduced cells (Figure 3A). In addition, we analyzed the combined expression and epigenetic data to compare the initial CD34+ cells and CD34+KM3 cells to highlight critical early events potentially implicated in the disease. Using the EaSeq package[20] we identified a set of genes (n=331) (Online Supplementary Table S2) with coordinated epigenetic and expression changes (Figure 3B). Among the genes with the most striking changes was the C-C motif chemokine receptor 1 (CCR1) gene, a G-protein coupled receptor family member. CCR1 is expressed in different types of blood cells, can bind a number of different ligands and is induced during normal monocyte development.[48] This gene was of specific interest because our previous work had shown that CCL23, a ligand only recognized by CCR1, is also a biomarker for KMT2A-translocated AML.[15] RNA-sequencing data from our models and genetically matched patients’ samples confirmed that the expression of both CCR1 and CCL23 is increased dramatically upon the integration of the KM3 fusion, while neither is expressed in normal CD34+ cells (Figure 3C). External RNA-sequencing data from the BloodSpot database (Online Supplementary Figure S9A) or TARGET cohort[4] also confirm the KMT2A-specific expression of CCL23 (Online Supplementary Figure S9B) and there is a small but significant difference in survival of pediatric AML patients with high levels of CCL23 (Online Supplementary Figure S9C). Although previous studies have highlighted the importance of CCR1 in different cancers (e.g., breast[49] and prostate[50] cancer) no studies have looked at its role in AML. Although we saw an inhibition of growth specifically in KM3-AML cell lines using a potent inhibitor of CCR1 (data not shown), this inhibition was only observed at levels (10-20 mM) much higher than the reported half maximal inhibitory concentration (IC50) (9 nM) of the inhibitor.[51] The very high expression levels of both the receptor and ligand may explain this observation; however, specific genetic models will be required to define the biological importance of this potential autocrine-signaling axis.

Figure 3.

Epigenetic analysis changes in (A) Venn diagram of peaks determined by MACS2 in different stages of the model system (CD34+, CD34+KM3, KM3 AML model) for H3K4Me3 and H3K79Me2 (n=3 for chromatin immumoprecipitation [ChIP]-sequencing of each mark). (B) A scatterplot of ChIP-sequencing enrichments for each histone mark in an area ± 2,000 bp (H3K4me3) or -1,000 to +3000 bp (H3K79me2) around the transcription start site for each gene is shown, with individual genes colored by the log2 ratio (CD34+/CD34+KM3) of expression for each gene. (C) Boxplots show the expression of CCR1 and its ligands in different stages of the model leukemias and in pediatric KMT2A-AML patients. The dotted green line represents an arbitrary threshold for an expressed gene (FPKM >1) (D) Annotation of ATAC-sequencing peaks unique to model leukemias (mAML) show that the majority are either in intronic or intergenic regions. (E) Hierarchical clustering of ATAC-sequencing peaks in published data from different normal blood cell types compared to leukemic cells. NK: natural killer cells; CD4T: CD4+ T lymphocytes; CD8T: CD8+ T lymphocytes; ERY: erythrocytes; AML: acute myeloid leukemia (model); Mono: monocytes; CD34 HSPC: CD34+ hematopoietic stem and progenitor cells; GMP: granulocyte-monocyte progenitor cells; B: B lymphocytes; CLP: common lymphoid progenitor; MEP: megakaryocyteerythroid progenitor cells; HSC: hematopoietic stem cells; MPP: multipotential progenitor cells; LMPP: lymphoid-primed multipotent progenitors cells; CMP: common myeloid progenitor cells. (F) A list of transcription factors is shown ranked by their abundance in the mAML and with their corresponding rank in the entire dataset and their expression level in the model system stages in the columns to the right. (G) UCSC snapshot showing an example of a gene with ATAC-sequencing peaks in CD34+KM3, mAML, CD34+ HSPC, hematopoietic stem cells and monocytes. The blue vertical line shows a peak that is unique to the mAML samples, whereas red vertical lines highlight common peaks that also overlap ENCODE transcription factor clusters.

KM3-mediated epigenetic changes are minimal relative to normal blood cells

Location analysis of the annotated ATAC-sequencing peaks unique to the model AML cells revealed the expected peaks at promoters of active genes, although the majority of open chromatin regions were either intergenic or intronic (Figure 3D). As a point of comparison for our ATAC-sequencing data, we downloaded previously published and normalized ATAC-sequencing data generated from normal blood cells[24] along with raw ATACsequencing data from normal CD34+ cells and monocytes from the BLUEPRINT project.[52] Using peaks identified via MACS2 in our ATAC-sequencing data or remapped published data, or through supervised thresholding of normalized UCSC tracks, we compared regions of open chromatin in normal and leukemic blood cells. As with the histone modifications, the majority of peaks identified in the model leukemia samples were shared with one or more normal types of blood cell. Considering the 10,415 peaks identified by MACS2 that were specific to the model AML (Online Supplementary Figure S10), analysis of potential transcription factor binding motifs once again identified factors, such as PU.1, that are critical for normal myeloid cell differentiation (Online Supplementary Figure S11). Hierarchical clustering of all ATAC-sequencing peaks from the CD34+KM3 and model AML revealed that the chromatin landscape in the model AML is close to that of normal monocytes and CD34+ HSPC (Figure 3E). The clear separation of CD34+KM3 cells, despite the small absolute number of unique peaks (and more peaks common to CD34+ HPSC than to monocytes), suggests that the chromatin of these cells is subtly affected by their culture in vitro.

Using data from the ENCODE consortium, we further examined the transcription factors that have been experimentally observed to bind to leukemia-specific open chromatin regions. For this, we ranked the abundance of individual factors in the entire dataset, and in the leukemia-specific peaks, to account for the difference in frequencies of binding sites for individual factors (Figure 3F). This analysis revealed an enrichment of factors in the leukemia-specific peaks that did not necessarily correlate with an increase in expression of the factor. Given the ChIP-sequencing data, this suggests that the novel, leukemia-specific peaks may create competition for the binding of a range of expressed transcription factors, potentially co-opting their activity. In other cases, while there was no difference in the ranking, the factor itself was more highly expressed; in cases in which such factors have been demonstrated to be essential for leukemia development (e.g., FOS[53]), the higher expression level may be important to ensure a high/continuous occupancy at the open regions in order to drive the leukemic transcriptional program.

Across the genome, we observed many examples of genomic regions in which ATAC-peaks common to both leukemic and normal cells were seen in the vicinity of other leukemia-specific peaks (Figure 3G). Given our observations, we questioned whether all of the leukemia-specific peaks we observed were potentially functional or whether they might represent stochastic changes in accessibility. To address this question, we compared data for chromatin accessibility (DNaseI hypersensitive sites) and ChIP-sequencing data from the ENCODE consortium to the ATAC-sequencing peaks that were unique to the model leukemias. Interestingly, the peaks determined by either MACS2 or supervised thresholding showed high levels of overlap (65% and 75%, respectively) with both previously defined transcription factor clusters, DNaseI hypersensitive sites and overlapped chromatin marks associated with enhancers seen in other cell types. These data argue that despite their unique presence in leukemic cells, the large majority of leukemia-specific open chromatin regions likely represent enhancers.

Discussion

The KMT2A-MLLT3 gene fusion is one of the most common genetic anomalies found in pediatric AML and is an epigenetic regulator that is essential, directly or indirectly, for the alterations of chromatin structure[13,54,55] and DNA methylation seen in the disease.[56] Despite its extensive characterization, less is known regarding its role in the initial transformation of HSPC into leukemic cells. Our current study leverages our previously described human model leukemia system to integrate epigenetic and expression changes to better characterize these changes. With respect to DNA methylation changes, our data confirm previously published results[57] showing that an aberrant hypomethylation profile is found in KM3 leukemias. Interestingly, relative to the initial CD34+ cells, CD34+KM3 and B-cell ALL have a hypermethylation phenotype, whereas the model leukemia and patients’ AML show the reverse. While the significance of all of these global differences is not yet clear, previous singlecell studies showed that differentially methylated promoter regions in normal blood cells were enriched in binding sites for transcription factors promoting myeloid differentiation.[44] Our data (e.g., Online Supplementary Figure S7) agree with this finding, which may reflect a default lymphoid differentiation pathway in the absence of demethylation of these regions. Recent studies have highlighted the importance of DNA methyltransferase in KM3-mediated leukemogenesis in which DNMT3A and DNMT3B loss-of-function mutations accelerate leukemia progression by enhancing cell cycle progression.[58] Forced expression of DNMT3B delayed KM3 leukemia development[59] whereas the deletion of DNMT1 prevented KM3 AML progression.[60] These findings and our own data demonstrate that the loss of de novo methylation, but not necessarily maintenance methylation, is required for the induction of AML.

Analysis of the correlated epigenetic and gene expression changes highlighted a number of genes of interest, including ADCY9, an atypical member of the adenylate cyclase family. ADCY9 shows low but relatively broad expression in most tissues except blood (of all lineages[61]) (Online Supplementary Figure S12), and is capable of inducing cAMP accumulation,[62] although its activation with G protein subunits appears to be cell-type specific. Our functional studies demonstrate that ADCY9 function is critical for the growth of KM3-AML cells, and levels of cAMP in leukemias[63,64] have been shown to be important. Nevertheless, whether the growth defects seen in response to the loss of ADCY9 expression are mediated through its ability to produce cAMP remains unclear and will require further study. The integration of our data also highlighted the potential importance of the chemokine receptor gene CCR1 and one of its ligands, CCL23. This receptor-ligand pair were noteworthy in part because of their strong level of upregulation, but also because the expression of both increased immediately after the addition of the gene fusion, prior to the complete transformation of the cells. The CCR1 receptor can recognize at least ten different ligands including CCL3 and CCL5, which are both ubiquitously expressed in normal blood cells. In contrast, CCL23 is only expressed at very low levels in promyelocytes[65] (Lin- FSChiSSCintCD34-CD15intCD49dhiCD33hiCD11b-CD16-) and granulocyte-monocyte progenitors[65] (Lin- CD34+CD38+CD45RA+CD123+) but is highly expressed in KMT2A-AML according to a number of independent datasets. The impact of CCR1 signaling in AML remains poorly defined, although one consequence appears to be activation of the Gai pathway resulting in the inhibition of cAMP production.[66] In this context, it is interesting to note that ADCY9 is specifically upregulated in KMT2AAML (Figure 2F), potentially as counterbalance to the activity of CCR1 and Gai and to maintain cAMP levels.

Perhaps most surprisingly with respect to chromatin structure and accessibility, while we expected to see widespread changes, in fact there were relatively few modifications specific to the model leukemia cells compared to normal blood cells. This observation agrees with those of a recent study showing that KMT2A-MLLT3 helps to preserve gene expression of the cellular states in which it is expressed.[67] As noted, most of the leukemiaspecific regions overlap open chromatin regions and transcript factor clusters seen in large-scale ENCODE datasets. Although speculative, these observations suggest the interesting hypothesis that there may be a finite number of regions in the genome that are competent to be open in vivo across all cell types. The addition of the KM3 fusion, while activating the expression of a subset of genes, may also cause re-localization of chromatin, leading to the opening of regions not normally seen in blood cells. Thus, while chromatin remodeling is required to give access to endogenously expressed transcription factors in some of the novel locations that we characterized in our study, others may simply be a response to this remodeling. The distribution of ATAC peaks that we observed, largely in introns and intergenic regions (Figure 3D), also raises the question of the potential biological role as novel enhancers. Similar distributions in ATACsequencing data have been noted[68] and intronic enhancers have been shown to regulate the expression of genes in a tissue-specific manner,[69] which may explain why many of the ATAC-sequencing peaks not present in blood cells nevertheless corresponded to open chromatin regions in other cell types. These results open a new avenue to look at some of these regions of potential interest in the context of AML, to better understand their role in enabling access to genes which may be essential for leukemic growth.

In conclusion, here we have leveraged the advantages of our model system to characterize the immediate impact of the KM3 fusion without the complicating factor of additional mutations and other genetic abnormalities typically found in patients’ genomes. While our results highlight the subtlety of the changes that occur once the oncogene is introduced, they also support the role of previous genes implicated in leukemogenesis and highlight a novel receptor/ligand pair with potential functional relevance for KMT2A-AML. The novel insight into the nature and extent of epigenetic reprogramming that occurs in the disease will help to guide future studies to identify novel potential therapeutic targets.