Automated CUT&Tag profiling of chromatin heterogeneity in mixed-lineage leukemia

Acute myeloid and lymphoid leukemias often harbor chromosomal translocations involving the KMT2A gene, encoding the KMT2A lysine methyltransferase (also known as mixed-lineage leukemia-1), and produce in-frame fusions of KMT2A to other chromatin-regulatory proteins. Here we map fusion-specific targets across the genome for diverse KMT2A oncofusion proteins in cell lines and patient samples. By modifying CUT&Tag chromatin profiling for full automation, we identify common and tumor-subtype-specific sites of aberrant chromatin regulation induced by KMT2A oncofusion proteins. A subset of KMT2A oncofusion-binding sites are marked by bivalent (H3K4me3 and H3K27me3) chromatin signatures, and single-cell CUT&Tag profiling reveals that these sites display cell-to-cell heterogeneity suggestive of lineage plasticity. In addition, we find that aberrant enrichment of H3K4me3 in gene bodies is sensitive to Menin inhibitors, demonstrating the utility of automated chromatin profiling for identifying therapeutic vulnerabilities. Thus, integration of automated and single-cell CUT&Tag can uncover epigenomic heterogeneity within patient samples and predict sensitivity to therapeutic agents.

Ten percent of acute leukemias harbor chromosomal translocations involving the KMT2Agene encoding lysine methyltransferase 2A (also referred to as mixed-lineage leukemia-1)1. In its normal role, KMT2A catalyzes methylation of the Lys4 residue of the histone H3 nucleosome tail (H3K4) and is required for fetal and adult hematopoiesis2. The N-terminal portion of KMT2A contains a low-complexity domain that mediates protein–protein interactions, an AT-hook/CXXC domain that binds DNA and multiple chromatin-interacting domains (PHD domains and a bromo domain), whereas the C-terminal portion contains a transactivation domain that interacts with histone acetyltransferases and a SET domain that catalyzes histone H3K4 methylation3,4,5,6. The KMT2A precursor protein is cleaved to a 320-kDa N-terminal fragment (KMT2A-N) and a 180-kDa C-terminal fragment (KMT2A-C) that form a stable dimer7,8.

KMT2A contributes to leukemogenesis through oncogenic chromosomal rearrangements involving the DNA-binding domain in the N-terminal portion of KMT2A with a diverse array of other chromatin-regulatory proteins9,10. Although more than 80 translocation partners have been identified in KMT2A-rearranged (KMT2Ar) leukemias, fusions involving the AF9, ENL, ELL, AF4 and AF10 transcriptional elongation factors account for the majority of cases1,10. These fusion partners regulate RNA polymerase II elongation (ELL and AF4), recruit the DOT1L H3K79 histone methyltransferase (AF10), or both (AF9 and ENL)11,12,13,14. Additionally, ENL and AF9 interact with the CBX8 chromobox protein to neutralize the Polycomb repressive complex 1 (PRC1) gene-silencing complex15,16.

Previous work has suggested that KMT2A fusion proteins bind different genomic loci depending on the fusion partner to drive different leukemia subtypes17,18. For example, AF4 fusions are more common in acute lymphoid leukemia (ALL), and AF9 fusions are associated with acute myeloid leukemia (AML)1. In addition, KMT2A rearrangements are also prevalent in mixed-phenotype acute leukemia (MPAL), and numerous examples of KMT2Ar leukemias that interconvert between lineage types have been documented17,19,20,21. However, because methods for efficiently and reliably profiling KMT2A fusion-binding sites in scarce input patient samples are lacking, the relationship among KMT2A fusions, chromatin structure and lineage plasticity has been challenging to fully characterize. Here we establish a chromatin profiling platform that efficiently profiles oncogenic fusion proteins, transcription-associated complexes and histone modifications in cell lines and patient samples. By integrating these results with findings from related single-cell methods, we characterize the regulatory dynamics of KMT2Ar leukemias. We identify groups of target genes for fusion oncoproteins that show divergent patterns of active and repressive chromatin within the same sample. These patterns suggest that KMT2A fusion proteins activate distinct oncogenic networks within different cells of the same tumor and may explain the lineage plasticity associated with KMT2Ar leukemia. In addition, we find that distinct fusion partners display different affinity for various transcriptional cofactors, which predicts cancer sensitivity to therapeutic compounds.

Mapping the binding sites of diverse KMT2A fusion proteins

Characterizing the chromatin localization of oncogenic fusion proteins has often been limited by the inability of chromatin immunoprecipitation and sequencing (ChIP–seq) to be used with small amounts of patient samples. To efficiently compare the binding sites for wild-type KMT2A and fusion proteins, we applied AutoCUT&RUN22 across a panel of four KMT2Ar leukemia cell lines and eight primary KMT2Ar patient samples sorted for CD45+blasts. This collection spans the spectrum of KMT2Ar leukemia subtypes with diverse KMT2Atranslocations that create oncogenic fusion proteins with the transcriptional elongation factors AF4 (SEM, RS4;11, 1° ALL-1, 1° MPAL-2), AF9 (1° AML-3, 1° MPAL-1), ENL (KOPN-8, 1° AML-2), AF6 (ML-2) and AF10 (1° AML-4, 1° AML-5) as well as a relatively rare fusion to the cytoplasmic GTPase SEPT6 (1° AML-1) (Supplementary Table 1). With the exception of ML-2, an AML-derived cell line, these samples also contained a wild-type copy of the KMT2A locus. For comparison, we also profiled KMT2A localization in untransformed human CD34+hematopoietic stem and progenitor cells (HSPCs), in H1 human embryonic stem cells and in the K562 leukemia cell line, each of which lacks KMT2A translocations. Antibodies to the C-terminal portion recognized only wild-type KMT2A-C, while antibodies to the N-terminal portion recognized both wild-type KMT2A-N and the fusion proteins (Fig. 1a). Therefore, binding sites unique to oncogenic fusion proteins could be identified by comparing the chromatin profiles obtained with antibodies specific to C-terminal and N-terminal KMT2A. We used AutoCUT&RUN to profile replicate samples with two different antibodies to the N terminus and two different antibodies to the C terminus of KMT2A, and correlation analysis of the sequencing results showed high reproducibility (r = 0.78 ± 0.19; Extended Data Fig. 1a).

a, A general strategy for mapping KMT2A fusion proteins. The wild-type KMT2A protein (black) is cleaved (white lines) into KMT2A-N and KMT2A-C proteins. Common oncogenic lesions (black arrowhead) produce in-frame translation of oncogenic KMT2A with numerous fusion partners (gray). Antibodies to C-terminal KMT2A (blue) recognize wild-type KMT2A-C. Antibodies to N-terminal KMT2A (red) recognize wild-type KMT2A-N and oncogenic KMT2A fusion proteins. b, Example of a wild-type KMT2A-binding site (EIF4E) and an oncoprotein-binding site (HOXA locus). Black scale bars, 10 kb. c, Scatterplot comparing KMT2A peak width and relative enrichment of the KMT2A N versus C terminus in control (CD34+ and K562) samples and KMT2Ar (SEM and ML-2) samples. d, Heatmap comparison of KMT2A signal over 10-kb windows centered on wide KMT2A-binding sites (top) and narrow KMT2A-binding sites (bottom) in CD34+ HSPCs and KMT2A fusion protein-binding sites (top) and wild-type KMT2A-binding sites (bottom) in SEM cells. e, Pie charts showing the fraction of wide KMT2A peaks (CD34+ and H1 cells) and KMT2A fusion-bound sites (KMT2Ar samples) overlapping transcriptional start sites (TSSs), gene bodies and intergenic regions. P values were computed using Fisher’s exact test. f, Box plots of KMT2A-N and KMT2A-C signal showing that the antibody to N-terminal KMT2A is enriched relative to the antibody to the C-terminal portion of KMT2A at fusion-binding sites. The center line indicates the median, box limits represent the first and third quartiles, and whiskers show all data within 1.5 times the interquartile range (IQR) of the lower and upper quartiles; outliers are not shown. P values were computed using a two-sample t test (two sided). CD34+, n = 131; K562, n = 65; ML-2, n = 144; SEM, n = 91; RS4;11, n = 92; KOPN-8, n = 192; 1° ALL-1, n = 156; 1° AML-1, n = 349; 1° AML-2, n = 423; 1° AML-3, n = 103; 1° AML-4, n = 270; 1° AML-5, n = 186; 1° MPAL-1, n = 248; 1° MPAL-2, n = 189. g, Principal-component analysis (PCA) of fusion oncoprotein-binding sites in KMT2Ar samples. The first two components are shown.

Full size image

To evaluate our KMT2A dual-antibody approach, we first compared the KMT2A N-terminal and C-terminal profiles between our KMT2A wild-type and KMT2Ar samples. As expected, in H1 and K562 cells and CD34+ HSPCs, KMT2A-N and KMT2A-C showed nearly identical patterns of enrichment across the genome (r = 0.82 ± 0.08; Fig. 1b and Extended Data Fig. 1a,b). Strikingly, in H1 cells, KMT2A binding was generally focused in narrow peaks directly over TSSs, whereas in K562 cells and CD34+ progenitors additional regions showed wide peaks of both KMT2A-N and KMT2A-C extending from TSSs across gene bodies. Many of the genes with a wide KMT2A distribution in CD34+ progenitors (for example, HOXA9, RUNX2, MEIS1, MEF2C) are master regulators of hematopoietic cell fate (Fig. 1c) and have previously been defined as KMT2A fusion oncoprotein targets in leukemias18,23. In the KMT2Ar leukemia samples, the correlation between the KMT2A N-terminal and C-terminal profiles was significantly lower than in the control samples (r = 0.53 ± 0.21; Extended Data Fig. 1a,c), and many of the wide KMT2A-bound regions showed an enriched KMT2A N-terminal signal relative to the KMT2A C terminus (Fig. 1b,c). To systematically define fusion protein-binding sites across our collection of samples, we used Gaussian mixture modeling to partition KMT2A peaks into two different distributions based on both the width of the KMT2A peaks and the enrichment-normalized ratio of KMT2A-N to KMT2A-C signal (KMT2A N/C score; Fig. 1c and Extended Data Fig. 2). In CD34+ HSPCs, 131 of 6,336 KMT2A-bound sites were called as wide peaks (mean = 8.2 ± 4 kb) and a two-component Gaussian mixture model failed to partition the KMT2A-bound sites based on N/C score (Fig. 1c,d and Extended Data Fig. 2a–e), suggesting that there is similar enrichment of KMT2A-N and KMT2A-C proteins, consistent with wild-type KMT2A binding. In comparison, in the SEM cell line encoding a KMT2A–AF4 fusion, 195 of 8,259 KMT2A-bound regions were called as wide (mean = 13.1 ± 10 kb), and about half of these wide peaks (91/195) were enriched for KMT2A-N relative to KMT2A-C, which we interpret as fusion oncoprotein-binding sites (Fig. 1c,d and Extended Data Fig. 2f–j). In line with this interpretation, 61 of 91 of the oncoprotein target genes we identified in SEM cells overlapped with KMT2A–AF4 target genes that were previously identified in SEM cells using ChIP–seq23 (Extended Data Fig. 2k). The ML-2 cell line has a deletion of the wild-type KMT2A allele and harbors only a KMT2A–AF6 fusion oncoprotein. As expected, the majority of KMT2A-bound sites had a high KMT2A N/C score (144/211; Fig. 1c and Extended Data Fig. 2l). This persistent localization of the KMT2A–AF6 fusion oncoprotein to chromatin in ML-2 cells demonstrates that binding of the oncoprotein is not dependent on wild-type KMT2A.

Next, we examined how the distribution of KMT2A differed between wild-type and oncofusion proteins. In CD34+ HSPCs, 81% of wide peaks overlapped a gene TSS, whereas in KMT2Ar samples significantly fewer fusion oncoprotein peaks overlapped a TSS (30%), and significantly more (50%) overlapped a gene body (Fig. 1e). In comparison, in the control H1 human embryonic stem cell line, only 15 of 17,000 KMT2A peaks were called as wide and these peaks were significantly less enriched on gene TSSs than the wide KMT2A peaks we identified in CD34+ HSPCs and significantly less enriched in gene bodies than the oncoprotein-binding sites in KMT2Ar leukemia samples (Fig. 1e and Extended Data Fig. 2m). This pattern of KMT2A fusion oncoproteins spreading across target gene bodies is consistent with previous reports23,24. By comparing enrichment of the KMT2A N terminus and C terminus across fusion oncoprotein-binding sites in all KMT2Ar samples, we found that in all cases the N terminus was significantly more enriched than the C terminus (Fig. 1f), and in all cases the fold difference between the N-terminal and C-terminal signal at oncoprotein-bound regions in KMT2Ar samples was greater than in the wide KMT2A-bound regions in CD34+ HSPCs (Extended Data Fig. 3a). Taking these findings together, we conclude that our KMT2A N-terminal versus C-terminal antibody multiplexing approach identifies regions bound by diverse KMT2A fusion oncoproteins.

We then compared oncoprotein target sites among leukemias with different KMT2A fusions. We found that 81 of 440 (~18%) of all fusion oncoprotein target genes were shared by five or more of the KMT2Ar leukemia samples we profiled, representing 12% of the total sequence space occupied by the fusion proteins (Extended Data Fig. 3b,c). As expected, the group of genes we identified as the most frequent KMT2A fusion targets across our collection of samples included several genes that are known to be required for KMT2Ar leukemia, including HOXA9, MEIS1, MEF2C, MBNL1 and JMJD1C25,26,27,28,29 (Extended Data Fig. 3c and Supplementary Table 2). PCA of KMT2A N/C scores across all oncoprotein-binding sites indicated that both the specific fusion partner and the myeloid versus lymphoid lineage bias of the tumor may influence tumor-specific localization of the oncofusion protein (Fig. 1g). For example, all KMT2A–AF4 samples clustered together in the PCA plot and grouped with a sample from a patient with ALL and a sample from a patient with MPAL. By contrast, the ALL cell line KOPN-8 carries a KMT2A–ENL fusion protein and partitioned away from KMT2A–AF4-bearing leukemias. Primary AML samples bearing KMT2A–AF9, KMT2A–AF10 and KMT2A–ENL fusions formed a second cluster, apart from the KMT2A–SEPT6-containing primary AML and the primary KMT2A–AF9-bearing MPAL sample. Thus, tumors bearing KMT2A–AF4 fusions share a distinct binding profile, but other oncofusion proteins such as KMT2A–ENL and KMT2A–AF9 also appear to be influenced by the lineage bias of the tumor.

Chromatin landscape of KMT2Ar leukemia samples

To economically characterize the global chromatin landscape of tumors at a scale that could be generally applied to patient samples, we developed AutoCUT&Tag, a modification of our previous AutoCUT&RUN robotic platform22. CUT&Tag takes advantage of the high efficiency and low background of antibody-tethered Tn5 tagmentation-based chromatin profiling relative to previous methods, such as ChIP–seq and CUT&RUN30. The standard CUT&Tag protocol requires DNA extraction before library enrichment by PCR. However, we recently developed conditions for DNA release and PCR enrichment without extraction (CUT&Tag-direct)31. In this improved protocol, a low concentration of SDS is used to displace bound Tn5 from tagmented DNA, and subsequent addition of the non-ionic detergent Triton X-100 quenches the SDS to allow for efficient PCR. This streamlined protocol makes CUT&Tag compatible with robotic handling of samples in a 96-well plate format and generates profiles with data quality comparable to that produced by benchtop CUT&Tag (r = 0.79 ± 0.093; Extended Data Fig. 4a–c).

To define the chromatin features around KMT2A fusion-binding sites, we used AutoCUT&Tag to profile the active chromatin modifications H3K4me1, H3K4me3, H3K36me3, H3K27ac and H4K16ac as well as initiating RNA polymerase II marked by Ser5 phosphorylation of the C-terminal domain (RNAP2S5p). In addition, we profiled the silencing histone modifications H3K27me3 and H3K9me3. Together, these eight modifications distinguish active promoters, enhancers, transcribed regions, developmentally silenced chromatin and constitutively silenced chromatin32 and provide a straightforward picture of the regulatory status of a genome (Fig. 2a). Replicate profiles for each mark in control CD34+ samples and KMT2Ar leukemia samples were very similar and were merged for further analysis (H3K27me3, r = 0.93 ± 0.051; H3K4me3, r = 0.96 ± 0.015; H3K4me1, r = 0.90 ± 0.037; H3K9me3, r = 0.83 ± 0.060; H3K27ac, r = 0.80 ± 0.077; H3K36me3, r = 0.95 ± 0.021; H4K16ac, r = 0.97 ± 0.012; RNAP2S5p, r = 0.77 ± 0.107) (Extended Data Fig. 4d).

We first compared the chromatin features associated with sites bound by wild-type KMT2A to sites bound by KMT2A oncofusion proteins across all samples. In line with localization of the KMT2A fusion proteins to actively transcribed genes, we found that the active promoter marks H3K4me3, RNAP2S5p and H3K27ac were all present at oncofusion protein-binding sites (Extended Data Fig. 5a–c). H3K4me3 was also enriched at some promoters in the ML-2 cell line (for example, LPO and LYZ in Fig. 2b), which lacks the KMT2A methyltransferase domain, indicating that another H3K4me3 methyltransferase is responsible. In comparison to sample-matched wild-type KMT2A-bound sites, H3K27ac was enriched at oncofusion protein-binding sites in all samples, but this difference was statistically significant in only 6 of 11 samples (SEM, KOPN-8, 1° AML-1, 1° AML-2, 1° AML-3 and 1° AML-5; Extended Data Fig. 5b). The H3K4me3 mark was significantly enriched at oncofusion protein-binding sites in five of the samples (SEM, RS4;11, 1° AML-1, 1° AML-2 and 1° MPAL-2) and significantly depleted in five of the other samples (1° ALL-1, 1° AML-3, 1° AML-4, 1° AML-5 and 1° MPAL-1; Extended Data Fig. 5c). Oncofusion protein-binding sites lacked H3K27me3 and H3K9me3 (Extended Data Fig. 5d,e) but were enriched in H3K4me1 and H3K36me3, both of which mark transcribed gene bodies, and this enrichment was significant in 9 of 11 and 6 of 11 of the KMT2Ar leukemia samples, respectively (Extended Data Fig. 5f,g). Enrichment of these marks is expected with mistargeting of KMT2A fusions to gene bodies33.

Histone modification profiling holds the potential to reveal similarities and distinctions between leukemias by reporting their transcriptional status. For example, H3K4me3 reports gene promoter activity and was enriched at marker genes that matched the immunophenotypic characterization of each leukemia (Fig. 2b). To determine how the global distribution of these marks varied between KMT2Ar leukemia samples, we first identified regions enriched for each modification in our collection of KMT2Ar leukemia samples as well as CD34+ HSPCs using the SEACR peak-calling method34 and performed PCA to cluster samples according to their modification-specific similarities. Overall, active chromatin features marked by H3K4me1, H3K4me3, H3K36me3, H4K16ac or RNAP2S5p clustered samples according to their ALL, AML or MPAL lineage designation (Fig. 2c–e and Extended Data Fig. 6a,b), suggesting that a similar repertoire of active genes is used in each leukemia subtype. By contrast, PCA based on H3K27ac or H3K27me3 CUT&Tag profiles partitioned samples into groups largely unrelated to leukemia subtype (Fig. 2f and Extended Data Fig. 6c), and only the 1° AML-1 sample was distinguished by H3K9me3 (Extended Data Fig. 6d). H3K27me3 is an epigenetically inherited histone modification that is linked to developmental progression as cells determine their identity. Thus, these distinct H3K27me3 leukemia landscapes may be related to hematopoietic transitions that are defective in each tumor.

We next examined the lineage-specific variation in gene and regulatory element usage as indicated by the global chromatin landscape of each of the marks we profiled by performing t-SNE of these elements followed by density peak clustering35. This analysis revealed that H3K4me1-marked regions were highly variable between lineage subtypes, with a substantial fraction of elements marked specifically in AML samples falling to one side of the t-SNE plot (8,221/56,267), ALL-specific elements partitioned to the other side of the plot (8,466/56,267) and CD34+ HSPC elements grouped in the middle (7,141/56,267) (Fig. 2g and Extended Data Fig. 6e,f). A fraction of both the AML- and ALL-specific elements were also marked by H3K4me1 in CD34+ cells and the primary MPAL samples we profiled (Fig. 2g and Extended Data Fig. 6e,f). This regulatory overlap implies that MPAL leukemias share features with both ALL and AML and that KMT2Ar leukemia samples maintain H3K4me1 at regulatory elements used during normal hematopoiesis.

While only about half of the H3K4me1-marked elements were similarly labeled across all samples (~50%, 25,973/56,267), a much larger fraction of H3K4me3 (~75%, 10,958/13,998) and H3K36me3 (~85%, 22,858/26,759) peaks were common across leukemia subtypes, indicating that these subtypes largely share gene expression repertoires (Fig. 2h,i and Extended Data Fig. 6g–j). Grouping H3K4me3-marked promoter regions by t-SNE also partitioned 64 AML- and 508 ALL-specific elements to opposite sides of the t-SNE graph and identified 1,918 elements that were shared with the MPAL samples and CD34+ HSPCs (Fig. 2h and Extended Data Fig. 6g,h); however, as compared to H3K4me1, where we identified 23,828 lineage-specific peaks among the 56,267 total peaks (~40%), a smaller proportion of H3K4me3-marked features showed any lineage specificity (2,490/13,998 peaks, or ~18%). This is consistent with previous reports that regulatory elements marked by H3K4me1 generally show more cell type specificity than promoter elements marked by H3K4me3 (refs. 36,37).

Similarly to the t-SNE analysis of H3K36me3-marked regions, t-SNE analysis of H3K27ac-, H4K16ac- and H3K9me3-marked regions did not partition the genome by lineage identity (Extended Data Fig. 7a–c). By contrast, both RNAP2S5p and H3K27me3 peaks showed diversity similar to that observed with H3K4me1 (Fig. 2j and Extended Data Fig. 7d). Analysis by t-SNE with H3K27me3 did not partition elements according to lineage subtype (Fig. 2j). Rather, AML and ALL samples had a significantly greater proportion of the genome that was marked by H3K27me3 than CD34+ cells (Fig. 2j and Extended Data Fig. 7e), suggesting that these tumor types are more differentiated. In line with this interpretation, MPAL samples had significantly fewer regions marked by H3K27me3 than the ALL or AML samples and were considered to have a higher degree of lineage plasticity (Fig. 2j and Extended Data Fig. 7e). We conclude that high-throughput CUT&Tag profiling provides a powerful tool to characterize KMT2Ar leukemias and that profiling the developmentally repressed genome reveals tumor-specific differences that are not apparent by profiling the active genome.

Bivalent chromatin at KMT2A fusion protein target sites

In addition to marking promoters that are engaged in active transcription, H3K4me3 is present at a limited subset of transcriptionally repressed ‘bivalent’ (that is, ‘poised’) promoters that are also marked by H3K27me3 (refs. 38,39). In our collection of leukemia samples, we observed both H3K4me3 and H3K27me3 at some promoters that were called as KMT2A fusion protein targets (Fig. 2a, left). Additionally, we observed genes that were bound by the oncofusion protein in the majority of KMT2Ar leukemia samples but were not called as targets in specific samples; we termed this group ‘missing targets’ (Fig. 2a, right). To systematically define bivalent promoters within our collection of samples, we quantified the abundance of H3K4me3 and H3K27me3 in 2-kb windows centered on gene TSSs of marked and unmarked promoters for each modification. By intersecting these groups, we identified approximately 2,000–5,000 bivalent promoters in each of the KMT2Ar leukemia samples (Fig. 2k and Supplementary Table 3). Interestingly, we found that approximately 33% (129/396) of promoters for missing targets were called as bivalent, whereas approximately 24% (267/1,097) of KMT2A fusion target promoters were bivalent and only 14% of wild-type KMT2A target promoters were bivalent (Fig. 2l). Thus, oncofusion protein target promoters are enriched for a bivalent chromatin signature, suggesting that expression of these genes may fluctuate among cells within a sample.

Cell-to-cell chromatin heterogeneity at KMT2A fusion targets

To test whether the bivalent chromatin signature at KMT2A fusion target promoters is due to heterogeneity among cells, we performed single-cell CUT&Tag on four KMT2Ar cell lines and four primary KMT2Ar leukemia samples. Antibody binding and pA-Tn5 tethering were performed on bulk samples, and individual cells were then arrayed in microwells on the ICELL8 platform for barcoded PCR library enrichment30. We optimized the median number of unique reads per cell while maintaining a high fraction of reads in peaks (FRiP) on the ICELL8 by varying the amount of SDS detergent used to release Tn5 after tagmentation and the amount of Triton X-100 used to quench SDS before PCR (Extended Data Fig. 8a,b). Using this approach, we profiled 1,137–3,611 cells for the H3K4me3, H3K27me3 and H3K36me3 histone modifications. After excluding cells with fewer than 300 fragments, single-cell CUT&Tag for H3K4me3, H3K27me3 and H3K36me3 yielded a median of 4,972, 13,025 and 3,962 unique reads per cell, respectively (Extended Data Fig. 8c). As a second quality-control step, we called peaks on the aggregate data of all cells profiled for each mark and removed cells that had a FRiP value below the normal distribution (Extended Data Fig. 8d,e). Profiles for each single cell were then split into 5-kb bins tiled across the genome, and cells were projected in UMAP space on the basis of this binning (Fig. 3a–c). Encouragingly, cells taken from the same leukemia sample and profiled in different experiments were clustered together in UMAP space, indicating that the data quality was consistent between batches (Extended Data Fig. 8f–h). This approach resolved clusters for samples based on H3K4me3 or H3K27me3 profiling but not H3K36me3 profiling (Fig. 3a–c). This implies that the leukemia samples differ in both sets of active promoters and silenced regions.

To examine intratumoral heterogeneity in the H3K4me3 and H3K27me3 signals, we first used the archR single-cell software package40,41 to calculate imputed gene scores for all genes according to the UMAP projection of all cells. We then determined the normalized dispersion of the imputed scores in cells from the same sample (Fig. 3d,e). Strikingly, bivalent missing targets showed significantly higher H3K4me3 dispersion in the SEM, KOPN-8, 1° ALL-1, 1° AML-2, 1° MPAL-1 and 1° MPAL-2 samples than tumor-matched control genes (Fig. 3f). This implies that levels of the H3K4me3 active promoter mark in these genes vary among cells within KMT2Ar leukemias.

Next, we examined variation in the repressive H3K27me3 mark at bivalent oncoprotein target genes. In 1° MPAL-1 cells, the normalized dispersion of H3K27me3 was significantly higher in bivalent missing target genes, and in the 1° MPAL-2 sample the normalized H3K27me3 dispersion was higher in bivalent target genes than in tumor-matched control genes (Fig. 3f). Some bivalent genes varied among cells for both the H3K4me3 and H3K27me3 modifications. For example, the HOXA9 gene was a missing target in 1° MPAL-1 cells (Fig. 2a) but showed high dispersion in both the H3K4me3 and H3K27me3 signals (Fig. 3d,e). Thus, bivalency of chromatin marks is associated with heterogeneity among cells within a sample.

Grouping bivalent target genes according to the Pearson correlation of their imputed gene scores across cells of a given leukemia sample separated two groups by either H3K4me3 or H3K27me3 profiling (Fig. 3g,h and Extended Data Fig. 9a–g). For example, the missing target gene HOXA9 had elevated H3K4me3 scores in a small fraction of 1° MPAL-1 leukemia cells (~15%) (Fig. 3i and Extended Data Fig. 9h). The TAPT1 gene clustered together with HOXA9(Fig. 3g) and, as expected, had the highest H3K4me3 scores in the same cells as HOXA9 (r = 0.87) (Fig. 3i). By contrast, genes that were anticorrelated with HOXA9, such as CPEB2 (r = −0.74) and MEIS1 (r = −0.42), had the weakest H3K4me3 signal in cells where HOXA9 was active (Fig. 3i). This suggests that there are two exclusive gene expression programs activated by KMT2A fusion oncoproteins. Furthermore, we found that the imputed H3K27me3 scores also formed inverse patterns of gene association from H3K4me3, where genes such as HOXA9 that were rarely marked by H3K4me3 in 1° MPAL-1 leukemia cells showed elevated H3K27me3 scores in the majority of tumor cells (~55%) (Fig. 3j and Extended Data Fig. 9i). These groups of divergent KMT2A fusion oncoprotein targets may contribute to the phenotypic plasticity of KMT2Ar leukemias.

AutoCUT&Tag profiling predicts drug sensitivity of leukemias

We reasoned that the distinct binding sites of KMT2A fusion proteins might be driven in part by the cofactors with which the fusion oncoproteins associate. Therefore, we used AutoCUT&Tag to map the distributions of ENL and DOT1L, two chromatin proteins that interact with KMT2A fusion proteins33. Regions bound by KMT2A fusion proteins were enriched for DOT1L and ENL in all samples as compared to sample-matched wild-type KMT2A-bound sites (Fig. 4a–c). DOT1L has been proposed to be a central component of oncogenic transformation by KMT2A fusion proteins in certain leukemias14,42, and we found that the DOT1L histone methyltransferase was significantly more enriched at oncofusion protein targets in KOPN-8, 1° AML-3 and 1° MPAL-1 samples than in the other leukemias we profiled (Fig. 4b). Both 1° AML-3 and 1° MPAL-1 carry a KMT2A–AF9 fusion protein, whereas KOPN-8 carries a KMT2A–ENL fusion, suggesting that the AF9 fusion partner recruits particularly high levels of DOT1L to oncoprotein target loci while other leukemias can be variable for DOT1L recruitment. This was also illustrated by the KMT2A–SEPT6 fusion (1° AML-1), where only modest enrichment of DOT1L and ENL at fusion-binding sites was observed.

Several studies have suggested that KMT2A–AF9-bearing leukemias are particularly sensitive to pharmacological inhibition of DOT1L methyltransferase activity43,44. We hypothesized that elevated DOT1L signal at oncoprotein target sites might be indicative of sensitivity to DOT1L inhibitors. Indeed, we found that KOPN-8 cells were more sensitive to the DOT1L inhibitor EPZ-5676 (half-maximal inhibitory concentration (IC50) = 22.45 µM) than either SEM (IC50 = 110.45 µM) or RS4;11 (IC50 = 29.71 µM) cells (Fig. 4d). Previous reports have also shown that KOPN-8 cells are sensitive to the DOT1L inhibitor EPZ-0007477 (ref. 44). We suspect that the EPZ-5676 IC50 values obtained here are higher than those previously reported for KMT2Ar cell lines because we only exposed cells to the inhibitor for 72 h rather than 4–10 d45,46. In line with the increased sensitivity of KOPN-8 cells to DOT1L inhibitors, after 72 h of EPZ-5676 treatment, we found a significant depletion of the active histone mark H3K27ac at oncofusion protein-binding sites in KOPN-8 cells but not at KMT2A–AF4-bound sites in the SEM or RS4;11 leukemia cell lines (Fig. 4e). Thus, this pharmacological agent specifically alters the chromatin of oncofusion protein targets in KMT2Ar leukemia samples.

To extend this analysis, we profiled the transcriptional scaffold protein Menin, which interacts with the N-terminal portion of KMT2A and with oncofusion proteins, by using AutoCUT&RUN. SEM, RS4;11 and KOPN-8 cells had similar levels of Menin at KMT2A fusion-bound sites, but the SEM cell line was more sensitive to the Menin inhibitor VTP50469 (IC50= 2.4 nM) than RS4;11 (IC50 = 18.7 nM) or KOPN-8 (IC50 = 14.7 nM) cells (Fig. 5a–c). We then used AutoCUT&Tag to profile H3K4me3 in VTP50469-treated cells (Fig. 5d). We called H3K4me3 peaks that showed a significantly depleted signal after drug treatment (Extended Data Fig. 10), and these depleted sites were highly enriched in the gene bodies of oncofusion protein targets (Fig. 5e,f). Finally, we examined chromatin accessibility and the presence of initiating RNA polymerase II (RNAP2S5p) in drug-treated cells using Pol-CUTAC47,48 and found that many oncofusion- and Menin-bound sites were normally highly accessible and bound by initiating RNA polymerase II (Fig. 5g). This supports the idea that oncofusion protein-induced transcription in KMT2Ar leukemias is highly sensitive to Menin inhibition.

Here we have applied high-throughput chromatin profiling to KMT2Ar leukemias to delineate fusion protein-specific targets and to identify chromatin features that are characteristic of myeloid, lymphoid and mixed-lineage leukemias. To economically profile these features, we took advantage of the high signal-to-noise ratio and low sequencing depth requirements inherent to CUT&RUN and CUT&Tag and fully automated both methods on a standard liquid handling robot. As CUT&Tag requires only thousands of cells for informative histone modifications47, AutoCUT&Tag is suitable for profiling of samples for a wide range of studies, including developmental and disease studies, and screening of patient samples. The enhanced throughput and consistency of the AutoCUT&RUN and AutoCUT&Tag platforms for chromatin profiling make these technologies suitable for profiling patient specimens.

By also performing AutoCUT&RUN on KMT2A fusions and components of the Super Elongation and DotCom complexes, we have elucidated the details of mechanisms that likely contribute to the heterogeneity of these tumors. We found that the most common KMT2A fusion proteins, including KMT2A–AF4, KMT2A–AF9, KMT2A–ENL and KMT2A–AF10, all colocalize with the DOT1L and ENL proteins in gene bodies. This suggests that interaction of the C-terminal domains of AF4, AF9, ENL and AF10 with transcriptional elongation complexes likely recruits fusion proteins from the promoter into the gene body. In line with the possibility that these interactions have a pivotal role in oncogenic transformation, wild-type ENL protein is required for tumor growth in numerous KMT2Ar cell lines49.

Using AutoCUT&Tag to profile histone modifications in leukemia samples, we identified frequent KMT2A fusion oncoprotein sites with bivalent chromatin features. At some sites, bivalent chromatin features correlated with heterogeneity among cells of the same tumor, which suggests that the heterogeneity in gene expression seen in populations of mixed-lineage leukemia cells is rooted in chromatin dynamics. We identify a group of KMT2A oncoprotein target genes that are shared in the majority of KMT2Ar leukemias but are missed in a subset of samples. In several of the KMT2Ar leukemia samples we profiled, these missing targets were among the genes that showed the highest variation in active and repressive chromatin marks within the tumor, suggesting that these missing targets may be bound and activated by the oncoprotein in a limited subset of cells within the tumor, causing them to fall below the levels necessary for detection in our bulk KMT2A profiling assays. This heterogeneity we observed at KMT2A oncoprotein target genes has implications for how resistance to therapies may develop, if only a subset of cells are susceptible to specific anticancer agents.

Heterogeneity in leukemias may arise if an early cancerous cell divides and differentiates into two related cell types. Alternatively, certain leukemias may sporadically switch between cell types21,50. Our single-cell profiling reveals that some leukemias display both active and repressive chromatin states at KMT2A fusion target loci that differ among individual cells. Kinetic analysis of chromatin dynamics within cell populations will be needed to determine whether bivalency reflects differentiation or sporadic switching, with implications for therapeutic strategies to limit relapse.

Multiple compounds targeting chromatin proteins have shown promise as therapeutics for certain leukemias44,51. Profiling the targets of these compounds distinguishes certain KMT2Ar leukemias in which DOT1L is enriched at KMT2A fusion oncoprotein target sites, thus providing a strategy for selecting patients for whom treatment with DOT1L-targeting compounds is suitable. We also identified samples where KMT2A fusion oncoprotein target genes are broadly enriched for H3K4me3 in gene bodies and also bound by initiating RNA polymerase II. These leukemias are particularly sensitive to treatment with the Menin inhibitor VTP50469 and again demonstrate the utility of chromatin profiling for selecting therapeutic treatments. Incorporating AutoCUT&RUN and AutoCUT&Tag into longitudinal clinical trials could thus provide a route to assess the efficacy of epigenetic medicines. In addition, these technologies are highly consistent between replicates and increase the number of samples that can be processed and sequenced in parallel by an order of magnitude relative to conventional chromatin profiling, suggesting that it could be feasible to apply AutoCUT&RUN and AutoCUT&Tag for patient diagnosis.