Integrative methylome-transcriptome analysis unravels cancer cell vulnerabilities in infant MLL-rearranged B cell acute lymphoblastic leukemia

Abstract

B cell acute lymphoblastic leukemia (B-ALL) is the most common childhood cancer. As predicted by its prenatal origin, infant B-ALL (iB-ALL) shows an exceptionally silent DNA mutational landscape, suggesting that alternative epigenetic mechanisms may substantially contribute to its leukemogenesis. Here, we have integrated genome-wide DNA methylome and transcriptome data from 69 patients with de novo MLL-rearranged leukemia (MLLr) and non-MLLr iB-ALL leukemia uniformly treated according to the Interfant-99/06 protocol. iB-ALL methylome signatures display a plethora of common and specific alterations associated with chromatin states related to enhancer and transcriptional control in normal hematopoietic cells. DNA methylation, gene expression, and gene coexpression network analyses segregated MLLr away from non-MLLr iB-ALL and identified a coordinated and enriched expression of the AP-1 complex members FOS and JUN and RUNX factors in MLLr iB-ALL, consistent with the significant enrichment of hypomethylated CpGs in these genes. Integrative methylome-transcriptome analysis identified consistent cancer cell vulnerabilities, revealed a robust iB-ALL-specific gene expression-correlating dmCpG signature, and confirmed an epigenetic control of AP-1 and RUNX members in reshaping the molecular network of MLLr iB-ALL. Finally, pharmacological inhibition or functional ablation of AP-1 dramatically impaired MLLr-leukemic growth in vitro and in vivo using MLLr-iB-ALL patient-derived xenografts, providing rationale for new therapeutic avenues in MLLr-iB-ALL.

Keywords: Epigenetics; Genetics; Leukemias; Oncology; Transcription.

Figure 1. Global DNA methylation status of the different iB-ALL subtypes.

(A) Violin plots reflecting the global DNA methylation levels of CpG sites identified by WGB-Seq. Graph represents the percentages of CpG methylation distribution of the genome segmented in 10 Kbp genomic windows. The DNA methylation status significantly differs among all iB-ALL subtypes and the average methylation of healthy BCPs. ***P < 0.001, 2-sided Wilcoxon’s rank sum test. (B) Circos-plot representation of DNA methylation levels along the genome. CpG methylation was averaged in 10 Mbp genomic windows and the average DNA methylation value for each iB-ALL subtype is represented as a histogram track. Inner lines identify MLL-AF4⁺ (blue) and MLL-AF9⁺ (green) translocation events. (C) Line plots depicting the DNA methylation profile for the FLT3 gene. The CpG context and the CpG site location are mapped to the corresponding genomic coordinates, as indicated in the lower panels. Areas with significant (q < 0.05) differential methylation between MLL-AF4⁺/MLL-AF9⁺ and non-MLLr iB-ALL or healthy BCPs methylomes are shaded. ***P < 0.001. (D) Venn diagrams representing the total number of DMRs with consistent hyper- or hypomethylation changes for each of the indicated comparisons. (E) Heatmap indicating the log₂ odds ratio enrichment of different DNA repetitive regions for significant hyper- or hypomethylated DMRs, respectively (q < 0.05).

Figure 2. Identification of differentially methylated sites in iB-ALL subgroups.

(A) Schematic indicating the number of samples analyzed with the Human MethylationEPIC microarray platform. (B) PCA for 758,932 CpG sites across all samples included in the DNA methylation study. (C) Heatmap representation depicting the methylation status of the 10,000 most variable CpG sites (y axis) for the different iB-ALL samples, BCPs, and naive B cells (x axis). (D) Barplot displaying the number of common (black) and specific (colored) significantly hyper- or hypomethylated CpG sites observed in the indicated comparisons (FDR < 0.05, mean β difference > 0.25). Venn diagram represents the number of overlapping dmCpGs between naive B cells and iB-ALL (β > 0.25). dmCpGs overlapping naive B cells and iB-ALL samples were discarded for downstream analyses because they represent methylation changes naturally occurring during B cell differentiation. The inset shows the total number of hyper- and hypomethylated CpG sites observed in each condition. (E) Pairwise Pearson’s correlation analysis indicating the degree of similarity between the different sample groups. A total of 77,596 dmCpGs observed at any iB-ALL vs. BCP condition were used for proper comparisons, and all comparisons were statistically significant (P < 0.001).

Figure 3. Different DNA methylation patterns in iB-ALL subgroups according to the genomic location and CpG context.

(A and B) Stacked barplots depicting the relative frequency of significant hyper- or hypomethylated CpGs in relation to their CpG context (A) or CpG location (B). (C) Bubble plots showing enrichment of TFBS under the indicated conditions as determined by the information obtained from the GTRD database. Top (yellow bubbles) and bottom (blue bubbles) graphs represent TFBS enriched in hyper- and hypomethylated CpGs, respectively. Bubble color denotes statistical significance (–log₂ adjusted P value). Dot size indicates the log₂ or enrichment of a particular TFBS data set as compared with the background distribution of the EPIC platform.

Figure 4. Aberrant DNA methylome in iB-ALL is associated with alterations at particular chromatin states.

(A) Heatmaps displaying histone mark enrichment analyses of hyper- and hypomethylated CpGs specific for naive B cells, each iB-ALL subtype. or common to a least 2 iB-ALL subtypes, as compared with healthy BCPs. Color scales represent odds ratio of significant dmCpGs obtained in previous analyses across 6 common histone modifications from the NIH Roadmap Epigenome consortium as compared with the background distribution of the Human MethylationEPIC platform. Legend indicates the type of B-lineage hematopoietic data sets used for comparisons. (B) Heatmaps representing chromatin state enrichment analyses of hyper- and hypomethylated CpGs in the aforementioned conditions. Color ranges indicate the odds ratio of the significant dmCpGs observed across 18 chromatin states obtained from the NIH Roadmap Epigenome consortium.

Figure 5. Gene coexpression network analysis reveals the AP-1 complex members JUN and FOS as potential pathogenic contributors in MLLr iB-ALL.

(A) Schematic depicting the number of samples analyzed by RNA-Seq. (B) PCA generated from the gene expression matrix (log₂ fragments per kilobase of transcript per million mapped reads [FPKM] values) across all samples analyzed by RNA-Seq. (C) Tanglegram representation of the cluster relationship between paired DNA methylation and gene expression data (entanglement = 0.31, cophenetic correlation score = 0.4). (D) Barplot indicating the number of common (black) and specific (colored) significantly up- or downregulated genes observed for each comparison versus healthy BCPs. Inner graph indicates the total number of DEGs under each condition. Venn diagram represents the number of overlapping DEGs between naive B cells and iB-ALL compared with healthy BCPs. P < 0.001, 1-tailed hypergeometric test. (E) Violin plots showing the coefficient of variation of the gene expression (log₂ fold change) of the DEGs for each condition compared with healthy BCPs. P < 0.001 for all comparisons, 2-sided Wilcoxon’s rank sum test. (F) Scatter plot depicting the normalized enrichment score (NES) and the contribution (represented by dot size) of particular gene modules in each group.

Figure 6. Integration of DNA methylation and gene expression data.

(A) Schematic depicting the workflow for integrating DNA methylation and RNA-Seq data using the ELMER algorithm. (B) Barplot depicting the number of gene expression–correlating hyper- or hypomethylated CpGs (absolute Pearson’s correlation > 0.5). (C) Violin plots showing the distribution of gene expression changes (log₂ fold change of the indicated groups versus healthy BCPs) for those genes with consistent correlation with dmCpG probes. (D and E) Scatter plots indicating the correlation between DNA methylation and gene expression for the genes LMO2 (D) and BRIP1 (E), respectively. Pearson’s correlation (Cor) score is indicated for each comparison. (F) Scatterplot of Reactome pathway enrichment analyses. Genes with consistent correlation with DNA methylation were used for enrichment calculation versus the background data set, which included all the genes with detectable expression in our RNA-Seq data set (18,668). Color range denotes the significance of the represented ontology (adjusted P value), while dot size indicates the ratio between the number of hits identified and the total number of hits in a given ontology. (G) Heatmap representation of enrichment of TFBSs in gene expression–correlating hyper- or hypomethylated CpGs. Color range indicates enrichment or underrepresentation of a given motif (log₂ odds ratio) as calculated by the ELMER algorithm.

Figure 7. E2F and AP-1 interacting factors control the methylation status of downstream target motifs.

(A) Ideogram representing the genomic location of E2F5 expression–correlating dmCpG sites. n denotes the number of correlating dmCpGs identified with ELMER algorithm. (B) Boxplot depicting the average methylation of the significant E2F5 expression–correlating CpG probes across different groups. (C) Boxplot reflecting the expression of E2F5 in the indicated groups. (D) Scatter plot showing the correlation between average DNA methylation of E2F5 motif targets and the expression of E2F5. Colored dots: blue, BCP; red, MLL-AF4⁺; green, MLL-AF9⁺; yellow, non-MLLr. (E) Violin plots indicating the distribution of gene expression changes (log₂ fold change of the indicated groups versus healthy BCPs) of target genes with E2F5 motifs obtained with ELMER algorithm.*P < 0.05; ***P < 0.001, 2-sided Wilcoxon’s rank sum test. All correlated genes with E2F5 motifs included in any of the iB-ALL subgroups were used for the representation of the B cell gene expression distribution. The “random” group includes a random sampling of the same number of genes included in the B cell group, but using the original gene expression matrix including all genes with detectable expression in the RNA-Seq data set. (F–J) Same as A–E, but for the TF FOSL2.

Figure 8. FOSL2 contributes to the modulation of the MLL-AF4⁺ methylome.

(A) Western blot confirming FOSL2 knockout in 2 clones of CRISPR/Cas9-edited SEM cells. Total protein refers to the loading control measured by Ponceau staining. (B) Violin plot displaying the average methylation landscape of SEM-WT and SEM-FOSL2^KO cells ***P < 0.001, 2-sided Wilcoxon’s rank sum test. (C) Histogram indicating the number of hyper- and hypomethylated CpG sites between SEM-FOSL2^KO and SEM-WT cells (mean β difference > 0.40). (D) HOMER motif enrichment analysis highlighting the top enriched TFBS identified in the context of hypermethylated CpGs in SEM-FOSL^KO cells. (E) Plot illustrating the enrichment analysis obtained from gene expression data from SEM-FOSL2^KO and SEM-WT and the gene set corresponding to the AP-1 pathway. The NES and the significance scores were calculated using the preranked mode of GSEA. (F) Violin plots displaying the distribution of gene expression changes (log₂ fold change of SEM-FOSL2^KO versus SEM-WT cells) of FOSL2 target genes (number of annotated dmCpGs > 3, n = 107) or a random subset on genes (n = 500), including all genes with detectable expression in the RNA-Seq data set. ***P < 0.001, 2-sided Wilcoxon’s rank sum test. (G) Barplots representing RNA-Seq expression results (in FPKM) of NUDT21 and FSCN1 genes in a pool of SEM-WT and SEM-FOSL2^KO cells. (H) Plots depicting the genomic localization and β value score of the interrogated MethylationEPIC probes in the vicinity of NUDT21 and FSCN1 genes. Colors denote the methylation landscape of the indicated groups, and MLL-AF4⁺ iB-ALL patient samples were included for proper comparisons. Significant dmCpGs (mean β difference > 0.40) are highlighted in red. Bottom panel displays the genomic context of the aforementioned probes and genes, including the location of FOSL2 motifs (GTRD database) and CpG islands. TSSs are highlighted in orange.

Figure 9. MLL-AF4 regulates the expression of AP-1 members and RUNX1 and reshapes the methylome landscape of CD34⁺ cells.

(A) Boxplot depicting the average expression of E2F5, FOSL2, andRUNX1 genes in healthy untransduced CD34⁺ cells or in CD34⁺ cells transduced with either human:murine chimeric MLL-AF4 or human MLL-AF9. ***P < 0.001, 2-sided Welch’s t test. n = 3; n = 6; n = 3, respectively. (B) UCSC Genome Browser tracks representing the binding pattern of MLL-AF4 (in CD34⁺ cells) or MLLN/AF4C (in SEM cells) in the vicinity of E2F5, FOSL2, or RUNX1 genes. Data represents ChIP-Seq signals obtained from NCBI database Gene Expression Omnibus [GEO] GSE84116 and GSE74812, respectively. (C) Barplots indicating RT-PCR relative fold-change of E2F5, FOSL2, or RUNX1 expression between nonedited CD34⁺ cells (CD34^control) and CRISPR-edited CD34⁺ cells carrying locus-specific t(4;11)/MLL-AF4⁺ (CD34^{CRISPR t(4;11)}). Barplots represent mean± SD. *P < 0.05, 2-sided Welch’s t test. (D) Scatterplot indicating the number of hyper- or hypomethylated CpG sites observed upon CRISPR t(4:11)/MLL-AF4⁺ edition in CD34⁺ cells (mean β difference > 0.20). (E) Venn diagrams representing the number of overlapping genes decorated with dmCpGs (top) or the number of overlapping dmCpGs (bottom, positional overlap) between MLL-AF4⁺ patients and CD34⁺ t(4;11)/MLL-AF4⁺ (CD34^{CRISPR t(4;11)}) cells, as compared with healthy BCPs or nonedited CD34⁺ cells, respectively. P < 0.001. One-tailed hypergeometric test was used for all the comparisons. (F) Plots displaying enrichment of TFBS in the context of hypo- and hypermethylated CpGs in CD34⁺ t(4;11)/MLL-AF4⁺ (CD34^{CRISPR t(4;11)}) cells as determined by the information obtained from the GTRD database. The y axis represents the –log₁₀ adjusted P value enrichment of particular TFBS data set as compared with the background distribution of the EPIC platform.

Figure 10. Targeting AP1 complex impairs the growth of MLL-AF4⁺ B-ALL cells in vitro and in vivo.

(A) Impaired in vitro proliferation of dnFOS-SEM cells. A representative experiment is shown (n = 3). Inset represents RT-PCR confirming dnFOS expression in transduced SEM cells. RPL19 was used as a loading control. (B) Impaired clonogenic capacity of dnFOS-SEM cells. (C) BLI imaging showing significant decreased tumor burden in NSG mice transplanted with dnFOS-SEM cells in vivo (n = 5 mice/group). (D) Endpoint macroscopic images of spleens from mice transplanted with WT- and dnFOS-SEM cells. (E) Endpoint leukemic burden in PB of mice transplanted with WT- and dnFOS-SEM cells. (F–J) Same as A–E, but for the CRISPR-edited FOSL2^KO-F103 SEM cells. *P < 0.05; **P < 0.01; ***P < 0.001, 2-sided Welch’s t test. Barplots represent mean±SD.

Figure 11. Pharmacological inhibition of the AP-1 complex specifically impairs the proliferation of MLLr leukemic cells.

(A) Dose-dependent impaired in vitro proliferation of CD34^{CRISPR t(4;11)} cells upon treatment with the AP-1 chemical inhibitors SR11302 (left) or T5224 (right). The DMSO condition was common for both treatments and is represented in both graphs for proper visualization purposes. (B and C) In vitro proliferation of several MLLr (B) and non-MLLr (C) leukemic cell lines treated with increasing doses of the T5224 inhibitor. Resulting P values were adjusted for multiple comparisons using Holm’s method. *P < 0.05; **P < 0.01; ***P < 0.001, 1-sided Welch’s t test. Line plots represent mean ± SD (n = 3), and color denotes the indicated concentration for the different drug treatments.

Figure 12. The AP-1 inhibitor T5224 synergistically cooperates with VXL to partially reduce leukemia aggressiveness in vivo using MLLr iB-ALL xenografts.

(A) In vivo experimental design of the preclinical MLLr iB-ALL PDX models treated with VXL alone or in combination with the T5224 inhibitor. MRD was evaluated in BM 1 week after completion of the 2 cycles of VXL-based chemotherapy (day 20). Subsequent relapses were monitored in BM, spleen, and PB. (B) BM leukemic engraftment determined before treatment (day 0) in order to equally randomize mice based on BM leukemic burden to receive VXL or VXL+T5224 treatment. (C) Dot plots monitoring MRD in BM 1 week after completion of the 2 cycles of VXL±T5224 (day 20). (D) Proportion of leukemic mice and levels of leukemic engraftment in BM, spleen, and PB at the end of the follow-up period (mice were sacrificed when blasts were detected in PB or at day 80). For all the graphs, horizontal black line represents the mean and dots represent individual animals. **P < 0.01, 1-sided Welch’s t test.