Genome-wide signatures of differential DNA methylation in pediatric acute lymphoblastic leukemia

The DNA methylation landscape in ALL

To obtain an initial view of the variation in CpG site methylation in our dataset, we subjected the complete set of methylation data to principal component analysis (PCA). T-ALL, BCP ALL, and non-leukemic samples formed separated clusters using the principal components 1 and 2 (Figure 1A). Although only two components were needed to capture >60% of the variation in the dataset (Figure 1B), higher order components separated the subtypes of BCP ALL from each other (not shown). Although the non-leukemic reference samples originated from different blood cell populations, they clustered together, clearly separated from the ALL samples. Unsupervised cluster analysis across all of the CpG sites revealed distinct methylation patterns that separated ALL cells according to their cytogenetic and immunophenotypic subtype. The evident difference between ALL cells and the non-leukemic blood cells, and the similarity between the non-leukemic cells in the heatmap (Figure 1C) provide the rationale to use these cells as a non-leukemic reference cell panel to detect differential methylation.

Differential DNA methylation

We searched for differentially methylated CpG sites (DMCs) in the ALL cells by comparing the β-values (methylation values ranging from 0.0 to 1.0) in non-leukemic reference samples to the ALL samples of each individual subtype. CD19+, CD34+, and BM samples were used as the reference panel for BCP ALL and CD3+, CD34+, and BM were used as the reference panel for T-ALL. For calling a CpG site as differentially methylated, we required a minimum absolute ∆β-value of 0.2 and a false discovery rate (FDR)-adjusted Wilcoxon rank-sum P-value of <0.01 for the difference. This analysis revealed between 21,799 and 58,157 DMCs in the ALL subtypes, distributed across 5,956 to 8,245 gene regions (Table 2; in Additional file 1: Table S1). In total, 9,406 of the DMCs annotated to 2,023 gene regions and 2,979 CpG islands were observed across all the ALL subtypes and were thereby considered 'constitutive' (Additional file 2: Table S2). The vast majority of the constitutive DMCs (98.6%) were hypermethylated in the ALL cells compared with the non-leukemic reference cells (Figure 2A). The number of DMCs that were unique for each ALL subtype according to the applied criteria varied independently of the number of samples in a subtype, from 16,841 CpG sites in 895 unique gene regions in T-ALL to 271 CpG sites in 36 unique gene regions in the t(9;22) subtype (Table 2). As expected, the heterogeneous BCP ALL samples with unknown cytogenetic aberrations labeled as 'undefined' and those with 'non-recurrent' abnormalities did not display unique differential methylation patterns. The methylation patterns between BCP ALL subtypes differed substantially, with high methylation levels in samples harboring MLL rearrangements, which is opposite to a recent finding of predominant hypomethylation in adult ALL with MLL rearrangements [9], while the high hyperdiploid (HeH) samples were predominantly hypomethylated in our study (Table 2; Figure 2B-I), as has been previously described in pediatric BCP ALL for HeH [11]. The distribution between hyper- and hypomethylation between the subtypes of pediatric BCP ALL in our study is in agreement with the findings in a recent study of 50,000 CpG sites that used an alternative method for DNA methylation analysis [16]. For the DMCs, the absolute average β-value difference between ALL cells and reference cells for the subtype-specific DMCs was approximately 0.50, which is in agreement with allele-specific gains or losses of DNA methylation in ALL compared to normal cells (Figures 2A-I; Additional file 3: Figures S1A-F).

DMC signature	DMCsa	Genesb	Genes unique	DMCs unique	Unique DMCs	Unique DMCs
(number of patients)			to subtypec	to subtype	+ DNAm (%)d	-DNAm (%)e
Constitutive (774)	9,406	2,023	NA	NA	NA	NA
T-ALL (101)	58,157	8,245	895	16,841	15,487 (92.0)	1,365 (8.0)
MLL/11q23 (28)	31,403	7,142	300	1,763	1,285 (72.9)	478 (27.1)
dic(9;20) (20)	53,680	9,009	202	2,370	1,561 (65.9)	809 (34.1)
HeH (187)	42,779	7,773	271	3,014	268 (8.9)	2,746 (91.1)
t(1;19)TCF3/PBX1(23)	21,799	5,956	107	1,110	272 (24.5)	838 (75.5)
t(12;21)ETV6/RUNX1(163)	45,589	7,973	156	2,114	1,126 (53.3)	988 (46.7)
t(9;22)BCR/ABL1(19)	23,871	6,047	36	271	140 (51.7)	131 (48.3)
iAMP21 (10)	44,726	8,614	272	2,656	997 (37.5)	1,659 (62.5)
Undefined (105)	39,262	7,059	3	56	8 (14.3)	48 (85.7)
Non-recurrent (100)	42,109	7,434	2	27	14 (51.9)	13 (48.1)

^aDifferentially methylated CpG sites (DMCs) with mean Δβ-values >0.20 and false discovery rate-corrected Wilcoxon rank-sum P-values <0.01.

^bThe number of gene regions to which the DMCs are annotated.

^cGenes with DMCs are considered unique to a subtype if only that subtype had significant DMCs in the gene region.

^dHypermethylated DMCs, + DNA methylation (DNAm).

^eHypomethylated DMCs, - DNA methylation (DNAm).

Functional genomic distribution of differentially methylated CpG sites

The hypermethylated DMCs were enriched in CpG islands, while hypomethylated DMCs were primarily annotated to 'open sea' regions, independent of whether they were constitutive or subtype-specific (Figure 3A). The subtype-specific differences were more frequently observed in CpG island 'shores' and 'shelves', which display a large variation in β-value between ALL samples (Additional file 3: Figure S2). Both constitutive and subtype-specific DMCs in proximal promoter regions (transcription start sites and 5’ untranslated regions) of genes were commonly hypermethylated, but a greater enrichment of subtype-specific hypomethylation was observed in gene bodies and in intergenic regions (Figure 3B). To explore putative functional roles for the DMCs, we intersected the genomic coordinates of the constitutive and subtype-specific DMCs with regions defined by chromatin-immunoprecipitation of six histone marks and DNase1 hypersensitivity (DHS) assays in relevant primary cell types such as CD19+, CD3+, and CD34+ cells [17, 18]. Although the histone code in normal blood cells may not reflect that in ALL cells, the genomic distribution of histone marks is useful for annotating functional regions of the genome. This analysis revealed differences in enrichment between constitutive and subtype-specific DMCs to functional genomic regions with marks of repressed or active chromatin (Figure 3C). The 9,406 constitutive DMCs were enriched more than two-fold in regions marked by repressive H3K9me3 and H3K27me3, or bivalently by H3K27me3 and H3K4me3, which marks active chromatin (P < 0.001). On the contrary, the subtype-specific DMCs were enriched more than two-fold in regions of active chromatin marked by DHS, H3K4me3, and H3K4me1 (P < 0.001; Figure 3C). These observations suggest that subtype-specific methylation of CpG sites has specific functional roles.

The constitutive DMCs were enriched in genes in the transcriptional regulatory network in embryonic stem cells (P = 3.53 × 10^-3) and in genes that regulate or are regulated by transcription factors involved in embryonic development: NANOG (P = 9.7 × 10^-6), OCT4 (P = 4.9 × 10^-5), SOX2 (P = 2.3 × 10^-6), and REST (P = 4.75 × 10^-13) (Additional file 2: Table S3). While no enrichment to known pathways was observed for the subtype-specific DMC signatures, all of the DMC signatures were enriched for genes with biological functions in cancer, cellular development, cellular growth and proliferation, and cell-to-cell signaling (P < 0.05).

DMCs as regulators of gene expression

To investigate whether the DMCs influence gene expression and to determine which of the annotation classes of DMCs are involved in the regulation of gene expression, we compared the DNA methylation levels of each constitutive and subtype-specific DMC with gene expression data. First, we determined the correlation between the methylation levels of constitutive DMCs and mRNA expression levels obtained using digital gene expression sequencing of 28 ALL samples, including T-ALL and five BCP ALL subtypes, and five reference samples [19] (Additional file 2: Table S4). The β-values of only a small proportion (<1%) of the constitutive DMCs (n = 85) correlated with up- or down-regulation of the mRNA expression levels of 41 genes (permuted P ≤ 0.05 and fold change ≥2) (Additional file 2: Table S5). This observation was expected since 79% of the constitutive DMCs were annotated to regions containing the repressive H3K27me3 or H3K9me3 marks in healthy blood cells and thus genes in these regions were presumably not widely expressed (Figure 3C). Secondly, we determined which of the subtype-specific DMCs correlated with microarray-based gene expression data for 93 of the ALL samples of the t(12;21), HeH, t(1;19), t(9;22), dic(9;20), MLL/11q23 and T-ALL subtypes (Additional file 2: Table S6). We found that, on average, 15% (range 10 to 21%) of the β-values for the subtype-specific DMCs annotated to genes correlated with gene expression levels (permuted P ≤ 0.05 and fold change ≥2) (Additional file 2: Tables S7 to S13). The proportion of DMCs and gene annotations in t(12;21) that were correlated with gene expression in our study were highly similar to those in a recent, small methylation study on the t(12;21) BCP ALL subtype [12]. Ten of the 17 genes suggested in the earlier study based on their correlation to be drivers of leukemogenesis were also highlighted in our study (Additional file 2: Table S14).

We used the functional annotation of the DMCs correlated with gene expression to explore their putative functional roles, and found hypermethylated DMCs that correlated with down-regulation of gene expression to be enriched in DHS regions, active promoters (H3K4me3), and enhancers (H3K27ac/H3K4me1) (Figure 4A; Additional file 3: Figure S3). On the contrary, hypomethylation of gene bodies was highly correlated with either up- or down-regulation of gene expression. DMCs that were highly correlated with gene expression included genes with functions in epigenetic regulation and previously known subtype-specific gene expression in ALL (Figure 4B). For example, we observed an inverse correlation between the β-value and gene expression for the UHRF1 gene, which encodes a methyl CpG binding protein that has high affinity for hemi-methylated DNA and was highly expressed in the ALL samples, independent of their subtype, while it was not expressed in reference samples [20]. DNA methylation of NCOR2, which is a transcriptional co-repressor that acts through covalent modification of histones [21], was positively correlated with gene expression in T-ALL. We also show up-regulation of known subtype-specific genes such as BIRC7 in t(12;21) [12, 22] and DDIT4L in HeH [19], and previously unobserved subtype-specific expression of PHACTR3 in t(1;19) and UAP1 in the dic(9;20) subtype.

Differential DNA methylation in relapsed ALL

Next we compared the genome-wide DNA methylation levels between paired samples at diagnosis and relapse from 27 patients, and in five of the patients after a second relapse (Additional file 2: Table S15). We used PCA to visualize the genome-wide methylation patterns of the sample pairs. Plots of the first two principal components showed similar changes in DNA methylation levels between diagnosis, first, and second relapse in all patients (Figure 5A; Additional file 3: Figure S4). We observed a similar pattern in 10 paired BCP ALL samples at diagnosis and relapse from the Quebec childhood ALL (QcALL) cohort that were included for verification of our results (Figure 5B).

In total, we identified 6,612 DMCs in 1,854 gene regions in the 27 paired diagnosis-relapse ALL samples (Additional file 2: Table S16). Although only 773 (12%) DMCs at relapse overlapped with the constitutive DMCs, the gene region annotations of both signatures were remarkably similar, and included 1,186 (64%) of overlapping gene regions. Hence, like the genes in the constitutive signature, the genes in the relapse signature were enriched for the transcriptional regulatory network in embryonic stem cells and in the Wnt/β-catenin signaling pathways (P = 2.8 × 10^-7, 1.8 × 10^-4; Additional file 3: Figures S5 and S6), to genes regulated by REST, SOX2, NANOG and OCT4 (P < 6.6 × 10^-10), and to regions with the repressive H3K27me3 mark or bivalent H3K4me3/H3K27me3 marks (P < 0.001; Figure 5C).

The methylation levels of each of the relapse DMCs increased in each of the ALL pairs, with the highest levels after the second relapse (Figure 5D). The β-values of the CpG sites in the relapse signature were highly similar in the Nordic Society of Pediatric Hematology and Oncology (NOPHO) and QcALL sample sets (Figure 5E), suggesting that this signature of DMCs is common to relapsed ALL samples, regardless of subtype and treatment protocol. To visualize individual β-value changes in the paired samples, the top 25 ranking DMCs from the relapse signature are plotted in the paired samples (Figure 5F). Regional analysis surrounding CpG sites in each of the top 25 genes showed that nearby CpG sites displayed concordant (increased) methylation changes at relapse (Additional file 3: Figure S7).

DNA methylation for predicting relapse-free survival in ALL

Finally, we utilized CpG sites that constitute the four signatures of differential methylation defined in this study to search for DMCs that are predictive of relapse-free survival of ALL patients. For this purpose, relapse-free survival in each ALL subtype further stratified into standard risk (SR), intermediate risk (IR), high risk (HR), and infant (I) treatment groups was analyzed against the β-values of the DMCs comprising the constitutive, subtype-specific, subtype-specific correlated with gene expression, and relapse signatures using nearest shrunken centroids classification (Additional file 2: Table S17; Additional file 3: Figure S8) [23]. Four of the methylation signatures allowed for prediction of relapse-free survival with an area under the receiver operating characteristic (ROC) curve (AUC) >0.60 (Figure 6A). After permutation testing, subtype-specific DMCs in the group of ALL patients with the t(12;21) translocation that had been treated according to the standard risk (SR) protocol (n = 71) were found to be associated with relapse (P = 0.033). In addition, the subtype-specific sites in patients with the t(9;22) translocation treated on the high risk (HR) protocol (n = 18) and 11q23/MLL rearrangements treated on the infant protocol (n = 14) had indicative P-values of 0.062 and 0.098, respectively, despite the small number of samples in these groups (Additional file 2: Table S17). The relapse signature in all patients treated according to the infant protocol was not statistically significant (P = 0.22).

DNA and RNA samples

BM aspirates or peripheral blood samples were collected from pediatric ALL patients enrolled in the NOPHO ALL92 or ALL2000 protocols [5]. Clinical follow-up data were obtained from the NOPHO registry. The median follow-up time for patients in continuous complete remission was 9.1 years (range 4.6 to 18 years). Lymphocytes were isolated from ALL samples at diagnosis (n = 764), remission (n = 86), first relapse (n = 27), and second relapse (n = 5) by Ficoll-isopaque centrifugation (Pharmacia, Uppsala, Sweden; Table 1). All samples included in the study contained >80% leukemic blasts at diagnosis (average 91%) and relapse (average 90%), and <5% at remission. For validation, a sample set of DNA samples that were isolated at diagnosis, remission, and relapse from 10 children with pediatric BCP ALL from the QcALL cohort was used. Clinical information for QcALL and relapse samples is available in Additional file 2: Table S15. CD19+ B cells and CD3+ T cells were isolated from peripheral blood mononuclear cells of healthy Swedish blood donors using positive selection (CD19 Microbeads #120-050-301 and CD3 Microbeads #130-050-101) and MACS cell separation reagents (Miltenyi Biotec, Bergisch Gladbach, Germany). Pooled CD34+ cells isolated from five healthy blood donors were purchased from 3H Biomedical (Uppsala, Sweden) [48]. DNA and RNA were extracted as previously described [19, 28]. The study was approved by the Regional Ethical Review Board in Uppsala, Sweden and was conducted according to the guidelines of the Declaration of Helsinki. The patients and/or their guardians provided informed consent.

DNA methylation assay

DNA was treated with sodium bisulfite (EZ DNA methylation Gold, Zymo Research, Irvine, CA, USA) and DNA methylation levels were measured using the Infinium HumanMethylation 450k BeadChip assay (Illumina, San Diego, CA, USA). The ALL samples and controls were randomly distributed across the arrays, all arrays were measured using the same HiScan instrument, and no evidence for batch effects was observed in the β-values (data not shown). The methylation β-value distribution between Infinium type I and II probes was normalized using peak-based correction (Additional file 3: Figure S12) [49]. The data were filtered by removing the data from probes on the X and Y chromosomes and with genetic variation affecting probe hybridization (Additional file 3: Figure S13). After filtering, methylation data for 435,941 CpG sites remained for further analysis (Additional file 1). A subset of diagnostic ALL samples (n = 364) were previously analyzed on a custom GoldenGate DNA methylation array (Illumina) [10]. DNA methylation values of 207 CpG sites interrogated by both arrays evaluate reproducibility of the β-value measurements (Additional file 3: Figure S14). Additional details about the methylation assay, probe filtering, and technical validation can be found in Additional file 4. The DNA methylation data are available at the Gene Expression Omnibus (GEO) with accession number GSE49031.

Annotation of CpG sites

CpG sites were annotated to RefSeq genes and CpG islands according to the Human Methylation 450k manifest file version 1.1. The distribution of probes that passed our stringent filtering is shown in relation to CpG islands, gene regions, and corresponding β-value distributions are shown in Additional file 3: Figures S15 and S16. When a CpG site had more than one gene-level annotation, that is, was present in both the transcription start site and the first exon, both annotations were used.

The following publicly available chromatin datasets from primary CD19+, CD3+, or CD34+ cells were obtained from the NIH Roadmaps Epigenomics Project: DHS regions, H3K27me3, H3K36me3, H3K4me3, H3K9me3, and H3K4me1 (Additional file 2: Table S18) [17]. Peaks were called using the MACS software using default settings [50]. H3K27ac peaks were downloaded from the UCSC table browser [51] derived from H1-hESC and GM12878 cell lines [18]. CpG sites were annotated for the chromatin marks by overlapping genomic location with a peak in at least two of the replicates analyzed (Additional file 2: Table S1).

Analysis of differential DNA methylation

DMCs were determined using the non-parametric Wilcoxon rank-sum test. They were determined in T-ALL using remission BM, CD3+, and CD34+ cells as reference and in BCP ALL using remission BM, CD19+, and CD34+ cells. The Wilcoxon signed-rank test was used to identify methylation differences between paired samples at diagnosis and relapse. Minimal cut-off values for the mean absolute differences in DNA methylation (∆β) of 0.2 were applied to highlight CpG sites with large differences between groups. CpG sites with standard deviations >0.10 in the reference control group (n = 33,533 sites) were removed from DMC lists to minimize DMCs occurring due to cell type-specific variability (Additional file 3: Figure S2).

Correlation between DNA methylation and gene expression

Genome-wide digital mRNA gene expression (DGE) sequencing data from 28 ALL patient samples and five non-leukemic reference samples were generated as previously described (Additional file 2: Table S4) [19]. RNA expression levels for 93 ALL patient samples were measured with Affymetrix U1333 Plus 2.0 arrays (Additional file 2: Table S6). Raw data were processed and normalized using the robust multichip average (RMA) algorithm [19, 52]. The expression datasets are publicly available at GEO under series GSE47051. Details on the gene expression assays can be found in Additional file 4. For each DMC signature, the correlation between β-value and log2 transformed gene expression was evaluated using the Pearson’s correlation coefficient. Statistical significance of each DMC was calculated by permuting the data 10,000 times and comparing the correlation coefficient in the unpermuted data to the permuted coefficients. In each dataset, the permuted P-values were adjusted for multiple testing using the Benjamini and Hochberg approach for controlling FDR [53].

Data analysis and visualization

Data analysis was carried out in the R environment [54]. The R code for the analyses performed in this study is available at GitHub [55]. One-sided Fisher’s exact tests were used to assess the significance of the enrichment of DMCs to functionally annotated regions, using the annotation of the 450k array as background. Pathway analysis and enrichment for upstream regulators was performed using software from Ingenuity Pathway Analysis (Ingenuity® Systems, Redwood City, CA, USA) and significance was evaluated with the Fisher’s exact test. All P-values were adjusted for multiple testing by FDR unless otherwise stated. Analysis of relapse-free survival for constitutive and relapse DMC signatures was performed on all patients. Relapse-free survival for the subtype-specific signatures was evaluated individually for T-ALL and BCP ALL separated into the cytogenetic subtypes 11q23/MLL, HeH, t(1;19), t(12;21), and t(9;22). Each subtype was further stratified according to standard, intermediate, high risk, or infant treatment protocols [5]. The patients with dic(9;20) and iAMP21 were not analyzed for relapse-free survival due to the small number of patients in each treatment group. Nearest shrunken centroids classifiers were designed to discriminate between the classes and evaluated with repeated cross-validation [23]. AUC was used to measure predictive performance and statistical significance was evaluated by permuting the data 1,000 times. Each CpG site was scored by its coefficient after shrinkage and the significance was evaluated by permutation testing, as described above. Further details on the relapse-free classification procedure can be found in Additional file 3: Figure S8 and Additional file 4.