DNA methylation and hydroxymethylation patterns in acute myeloid leukemia patients with mutations in DNMT3A and IDH1/2 and their combinations

Abstract

BACKGROUND:

Aberrant epigenetic patterns are a hallmark of acute myeloid leukemia (AML). Mutations in profound epigenetic regulators DNMT3A and IDH1/2 often occur concurrently in AML.

OBJECTIVES:

The aim was to analyze DNA methylation, hydroxymethylation and mRNA expression profiles in AML with mutations in DNMT3A and IDH1/2 (individually and in combinations).

METHODS:

Infinium MethylationEPIC BeadChip (Illumina) covering 850,000 CpGs was utilized. The validation of hydroxy-/methylation data was done by pyrosequencing. HumanHT-12 v4 Expression BeadChip (Illumina) was used for expression examination.

RESULTS:

Hierarchical clustering analysis of DNA hydroxy-/methylation data revealed clusters corresponding to DNMT3A and IDH1/2 mutations and CD34+ healthy controls. Samples with concurrent presence of DNMT3A and IDH1/2 mutations displayed mixed DNA hydroxy-/methylation profile with preferential clustering to healthy controls. Numbers and levels of DNA hydroxymethylation were low. Uniformly hypermethylated loci in AML patients with IDH1/2 mutations were enriched for immune response and apoptosis related genes, among which hypermethylation of granzyme B (GZMB) was found to be associated with inferior overall survival of AML patients ( $P=$ 0.035).

CONCLUSIONS:

Distinct molecular background results in specific DNA hydroxy-/methylation profiles in AML. Site-specific DNA hydroxymethylation changes are much less frequent in AML pathogenesis compared to DNA methylation. Methylation levels of enhancer located upstream GZMB gene might contribute to AML prognostication models.

1. Introduction

Acute myeloid leukemia (AML) is a predominantly fatal hematological malignancy and displays great molecular heterogeneity. Numerous genomic next-generation sequencing based studies have greatly contributed to a more comprehensive understanding of AML pathogenesis.  However,  AML  exhibits  less

Table 1

Characteristics of AML samples examined by DNA hydroxy-/methylation and gene expression profiling

Sample	Sex	%	NPM1	Activated signalling					DNA methylation				Myeloid	Spliceosome
		blast	activated										transcription
			signalling										factor
			NPM1	FLT3-ITD	FLT3-TKD	NRAS	JAK2	CBL	DNMT3A	IDH1	IDH2	TET2	CEBPA	SRSF2
			(%VAF)	length	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)	(%VAF)
DNMT3Amut_1	F	58	A 37.2%	$+$ 30 bp	negative	negative	negative	negative	R882H 46.4%	negative	negative	negative	negative	negative
DNMT3Amut_2	F	53	negative	$+$ 51 bp	negative	negative	negative	negative	R882H 43.7%	negative	negative	negative	negative	negative
DNMT3Amut_3	F	58	A 38.8%	$+$ 21 bp	negative	negative	negative	negative	R882H 46.8%	negative	negative	negative	negative	negative
DNMT3Amut_4	F	71	A 40.5%	$+$ 57 bp	negative	negative	negative	negative	R882C 46.6%	negative	negative	negative	negative	negative
DNMT3Amut_5	F	52	negative	$+$ 156 bp	negative	negative	negative	negative	R882H 46.4%	negative	negative	M1701I 46.8%	negative	negative
DNMT3Amut_6	M	54	A 39.0%	$+$ 81 bp	negative	negative	negative	negative	R882H 43.0%	negative	negative	negative	negative	negative
DNMT3Amut_7	M	65	negative	negative	negative	negative	negative	negative	R882H 44.7%	negative	negative	negative	negative	negative
DNMT3Amut_8	F	64	negative	$+$ 42 bp	negative	negative	negative	negative	R882H 49.8%	negative	negative	K536fs 49.6%	negative	negative
IDH1mut_1	F	84	A 39.0%	negative	negative	G12V 21.2%	negative	Y368D 2.3%	negative	R132C 42.0%	negative	negative	negative	negative
IDH1mut_2	M	89	D 43.4%	negative	D835V 29.5%	G12D 2.3%	negative	negative	negative	R132H 39.4%	negative	negative	negative	negative
					D835Y 7.2%
IDH1mut_3	F	65	A 10.8%	$+$ 57 bp	negative	Q61K 36.8%	negative	negative	negative	R132H 37.8%	negative	negative	negative	negative
IDH1mut_4	M	65	negative	negative	negative	negative	negative	negative	negative	R132C 33.4%	negative	negative	negative	negative
IDH2mut_1	F	90	A 42.3%	$+$ 18 bp	negative	negative	negative	negative	negative	negative	R140Q 44.5%	negative	negative	negative
IDH2mut_2	M	64	A 43.5%	negative	negative	negative	negative	negative	negative	negative	R140Q 42.8%	negative	negative	P95delinsRP
														35.3%
IDH2mut_3	M	84	A 41.5%	negative	negative	negative	negative	negative	negative	negative	R140Q 46.1%	negative	negative	negative
IDH2mut_4	M	75	A 37.3%	negative	negative	negative	negative	negative	negative	negative	R140Q 48.1%	negative	negative	negative
DNMT3A&IDH1mut_1	M	50	A 42.0%	negative	D835Y 20.2%	negative	negative	negative	R882H 45.2%	R82K 49.5%	negative	negative	P197delinsHPP	negative
													21.0%
DNMT3A&IDH1mut_2	F	78	A 44.5%	negative	negative	G13D 2.0%	negative	negative	R882C 45.3%	R132S 47.5%	negative	negative	negative	negative
						Q61K 2.9%
DNMT3A&IDH1mut_3	F	86	negative	negative	negative	negative	negative	negative	R882H 32.4%	R132C 29.9%	negative	negative	negative	negative
DNMT3A&IDH1mut_4	F	62	A 37.3%	negative	negative	G13D 38.0%	negative	negative	R882H 44.7%	R132H 35.0%	negative	negative	negative	negative
DNMT3A&IDH2mut_1	M	86	negative	$+$ 36 bp	negative	negative	negative	negative	R882H 42.7%	negative	R140Q 46.4%	negative	negative	negative
DNMT3A&IDH2mut_2	M	50	B 18.9%	negative	negative	negative	negative	C416S 4.9%	R882C 40.5%	negative	R140Q 36.6%	negative	negative	negative
								R420Q 67.6%
DNMT3A&IDH2mut_3	M	57	negative	negative	negative	negative	negative	negative	R882C 45.3%	negative	R140Q 43.4%	negative	negative	negative
DNMT3A&IDH2mut_4	M	64	negative	negative	negative	negative	negative	negative	R882P 29.0%	negative	R140Q 34.0%	negative	negative	negative

Variant allele frequency (VAF) and position of mutations present in selected acute myeloid leukemia (AML) diagnostic samples is shown. AML patients chosen for this study were selected from 258 consecutive non-APL intensively treated AML patients that were examined for their mutational status by ClearSeq AML kit (Agilent Technologies, USA) in the concurrent study [Folta et al., Br. J. Haematol. 2019, doi: 10.1111/bjh.15916, in press]. F – female, M – male, VAF – variant allele frequency.

genetic lesions when compared to most solid tumors [1]. Presumably other factors, such as epigenetic changes, may significantly contribute to initiation and progression of AML. Indeed, mutations of genes involved in DNA methylation and hydroxymethylation pathways are found in approx. 40% of AML patients at diagnosis [1]. Here, we focused on mutations in DNA methyltransferase 3A gene (DNMT3A ${}^{\text{mut}}$ ) and isocitrate dehydrogenases types 1 and 2 genes (IDH1/2 ${}^{\text{mut}}$ ) and their effect on DNA hydroxy-/methylation. DNMT3A ${}^{\text{mut}}$ is connected with DNA hypomethylation, whereas IDH1/2 ${}^{\text{mut}}$ are involved in hypermethylation signature [4, 9]. Furthermore, DNMT3A ${}^{\text{mut}}$ is often present concomitantly with either IDH1 ${}^{\text{mut}}$ or IDH2 ${}^{\text{mut}}$ and little is known about effects of these opposing mutations on DNA hydroxy-/methylation status. Just recently, the antagonism of DNMT3A ${}^{\text{mut}}$ and IDH1/2 ${}^{\text{mut}}$ has been described for the first time [8]. Hence, we aimed to examine influence of clearly defined molecular background on DNA hydroxy-/ methylation and expression patterns in AML patients and identify genes with biologically relevant changes in terms of potential prognostic impact.

2. Materials and methods

2.1 Patients

This study conforms with The Code of Ethics of the World Medical Association, it was approved by the Institutional Ethics Committee and all patients provided their full consent. Twenty-four de novo cytogenetically normal (CN) AML patients selected from the concurrent collaborative mutational study [Folta et al., Br. J. Haematol. 2019, doi: 10.1111/bjh.15916, in press; for more details see a legend of Table 1] were divided according to their mutational status into 5 groups: DNMT3A ${}^{\text{mut}}$ ( $n=$ 8), IDH1 ${}^{\text{mut}}$ ( $n=$ 4), IDH2 ${}^{\text{mut}}$ ( $n=$ 4), DNMT3A&IDH1 ${}^{\text{mut}}$ ( $n=$ 4) and DNMT3A&IDH2 ${}^{\text{mut}}$ ( $n=$ 4). As shown in Table 1, samples were picked up based on the variant allele frequency (VAF) of particular mutations with preference for high VAF in the same codon for each mutation.

2.2 Sample preparation

DNA and RNA from peripheral blood leukocytes were extracted using AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany). Mononuclear cells of healthy blood donors were harvested from buffy coats by Ficoll gradient centrifugation (Histopaque, Sigma-Aldrich). CD34+ cells were then isolated using MicroBeads kits (Miltenyi Biotec, Bergish Gladbach, Germany) and DNA was extracted using MagCore system (RBCBioscience, New Taipei City, Taiwan).

2.3 Methylation and hydroxymethylation assessment

DNA (1.2 $\mu$ g) of all samples was processed through both bisulfite (BS) and oxidative bisulfite (oxBS) conversion using TrueMethyl-Seq kit (CEGX, Cambridge, UK). For each sample, BS and oxBS-converted DNA (140–250 ng) was analysed on Infinium MethylationEPIC BeadChip (Illumina, CA, USA). Raw IDAT files were processed with R package ChAMP. $\beta$ -values (range 0–1, correspond to methylation percentage) were normalised using Functional normalisation [6] and corrected for batch effect with ComBat function [10]. OxBS $\beta$ -values were used for assessing DNA methylation. Difference between BS and oxBS $\beta$ -values ( $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}})$ represented DNA hydroxymethylation. Hierarchical clustering analysis was performed with R package pvclust [14].

2.4 Pyrosequencing

BS and oxBS converted DNA (10–20 ng) was amplified using PyroMark PCR kit (Qiagen). For primer sequences see Table S1. Pyrosequencing was performed on PyroMark Q24 instrument (Qiagen).

2.5 Gene expression analysis

Expression profiles of all AML samples and CD34+ cells from healthy controls ( $n=$ 6) were established using HumanHT-12 v4 Expression BeadChip (Illumina). R packages limma, Biobase and sva were used for expression data analysis. Raw IDAT files were processed by NormExp background correction and quantile normalization was utilized. Batch effect was corrected using ComBat function [10]. Expression levels of GZMB (Hs188051_m1) and CHFR (Hs00943495_m1) were assessed with TaqMan Gene Expression Assays (Life Technologies, MA, USA) and GAPDH (Hs02758991_g1) was used as a reference gene. Amplification was carried out with TaqMan Universal MasterMix II (Life Technologies) using StepOne instrument (Life Technologies).

Table 2

Summary of methylation and hydroxymethylation results

Figure 1.

Hierarchical clustering of samples based on DNA hydroxy-/methylation and gene expression data. Hierarchical clustering of DNA hydroxy-/methylation and gene expression was performed by pvclust package [14] using multiscale bootstrap resampling ( $n=$ 1,000 replications), the correlation method for distance calculation and Ward’s clustering method. A. Hierarchical clustering of DNA methylation data represented by oxBS $\beta$ -values passing quality filters (716,847 CpGs) revealed 3 main clusters according to presence of either DNMT3A ${}^{\text{mut}}$ or IDH1/2 ${}^{\text{mut}}$ and healthy controls. The concurrent presence of DNMT3A ${}^{\text{mut}}$ with IDH1/2 ${}^{\text{mut}}$ resulted in preferential clustering to healthy controls’ group. B. Hierarchical clustering of $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}$ values of 536,617 probes ( $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}<$ 0.03335304 were taken as unhydroxymethylated and assigned value of 0.0001; probes with undetectable hydroxymethylation for all samples were excluded) separated DNMT3A ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ patients from healthy donors, IDH1 ${}^{\text{mut}}$ and most of DNMT3A&IDH1/2 ${}^{\text{mut}}$ patients. C. Hierarchical clustering of gene expression data represented by 10,067 probes (detection $p$ -value $<$ 0.05) showed a cluster specific for healthy controls samples. Otherwise, no clearly defined clusters of AML with particular mutations existed. Supportively, technical replicates of 2 samples from two different arrays (DNMT3A ${}^{\text{mut}}$ 1-1&2; DNMT3A ${}^{\text{mut}}$ 4-1&2) clustered next to each other.

2.6 Statistics

Kaplan-Meier curves and two-sided log-rank test were used to estimate the overall survival (OS) and event-free survival (EFS) using GraphPad Prism 7 (GraphPad Software, La Jolla, CA, USA) at a level of significance of 0.05. Cox regression analysis was performed in the SPSS software (IBM, Armonk, NY, USA). Validation of hydroxy-/methylation data and correlation of DNA methylation and expression datasets was assessed based on Pearson correlation coefficient (R) and F-test of overall significance.

3. Results and discussion

3.1 Analysis of EPIC array data

3.1.1 DNA methylation

After quality filtering steps, 716,847 probes for each sample remained for a subsequent analysis. Hierarchical clustering revealed 3 main DNA methylation clusters (Fig. 1A) – DNMT3A ${}^{\text{mut}}$ specific, IDH1/2 ${}^{\text{mut}}$ specific and healthy controls cluster. Majority of DNMT3A&IDH1/2 ${}^{\text{mut}}$ samples (63%) clustered preferentially to the healthy controls’ group.

Numbers of differentially methylated positions (DMPs) are depicted in Table 2A. We observed a profound hypomethylation in DNMT3A ${}^{\text{mut}}$ and a hypermethylation tendency in IDH1 ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ , which is in line with previous observations [4, 9]. DNMT3A&IDH1/2 ${}^{\text{mut}}$ samples displayed fewer DMPs (less than a half compared to DNMT3A ${}^{\text{mut}}$ and IDH1/2 ${}^{\text{mut}}$ alone) with a more balanced hypo- and hypermethylation ratio in DNMT3A&IDH2 ${}^{\text{mut}}$ and a hypomethylation trend in DNMT3A&IDH1 ${}^{\text{mut}}$ . This has confirmed the antagonistic effect of these two opposing mutations described recently [8]. IDH2 ${}^{\text{mut}}$ patients demonstrated considerably higher numbers of DMPs in comparison with IDH1 ${}^{\text{mut}}$ patients (37,349 vs. 13,830 DMPs), which is in accordance with a thesis that these two mutations may represent, to some extent, different biological entities [16]. On the other hand, using the on-line annotation tool GREAT [12] we discovered that apoptosis and immune response pathways were affected within all samples with IDH1 ${}^{\text{mut}}$ or IDH2 ${}^{\text{mut}}$ even when present together with DNMT3A ${}^{\text{mut}}$ (for a list of 215 common hyper-DMPs in bed format – see Table S2). GZMB gene (encoding for granzyme B that is crucial for apoptosis induction [2]) was among the apoptosis-related genes universally hypermethylated in all IDH ${}^{\text{mut}}$ patients. Another gene, CHFR (checkpoint with forkhead- associated and RING finger domains), was found hypermethylated in IDH1 ${}^{\text{mut}}$ samples. CHFR is a mitotic checkpoint and tumors suppressor gene, which was recently described as having a negative prognostic value when downregulated and/or hypermethylated [7, 17].

3.1.2 DNA hydroxymethylation

$\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}$ (corresponding to hydroxymethylation levels) was calculated for each probe on the 850K array that passed quality filters. Threshold for hydroxymethylation detection was set to 0.03335304 based on the 95th percentile of negative $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}$ values across all examined samples [11]. Lower $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}$ were considered as unhydroxymethylated. Hierarchical clustering of all samples based on probes for which at least one sample had $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}>$ 0.03335304 ( $n=$ 536,617 CpGs) separated DNMT3A ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ patients from healthy donors, IDH1 ${}^{\text{mut}}$ and most of DNMT3A&IDH1/2 ${}^{\text{mut}}$ patients (Fig. 1B). Surprisingly, hydroxymethylation patterns of IDH1 ${}^{\text{mut}}$ patients were more similar to healthy controls and patients bearing DNMT3A&IDH1/2 ${}^{\text{mut}}$ .

A rather conservative threshold of $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}⩾$ 0.2 was set for detection of biologically relevant DNA hydroxymethylation sites. Numbers of differentially hydroxymethylated CpGs (hmCpGs) were on average two orders of magnitude below the numbers of aberrantly methylated CpGs as seen in Table 2B. DNMT3A ${}^{\text{mut}}$ patients possessed the lowest number of hmCpGs. This might be explained by the relative absence of potential sites for hydroxymethylation (i.e. the observed extreme hypomethylation in DNMT3A ${}^{\text{mut}}$ ). On the contrary, IDH2 ${}^{\text{mut}}$ patients harboured the highest number of hmCpGs. This is somewhat surprising considering that mutant IDH isoforms cause inhibition of TET2-mediated hydroxymethylation [4]. However, driving a clear assumption is impeded by generally low numbers of hmCpGs detected in blood (compared to brain samples) [13]. And again, there was a prominent difference between IDH1 ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ in terms of hmCpG numbers (715 vs. 1,931 hmCpGs). Interestingly, MYB, a well-known proto-oncogene, was among the genes both hydroxymethylated and upregulated in IDH2 ${}^{\text{mut}}$ patients (Fig. S1). We also detected RNF216 hydroxymethylation in all AML samples with exception of DNMT3A ${}^{\text{mut}}$ . RNF216 (ring finger protein 216) plays a fundamental role in cellular homeostasis and its upregulation was associated with progression of colorectal cancer [15]. For the list of hmCpG sites ( $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}⩾$ 0.2) in particular AML subgroups see Table S3.

3.2 Validation of DNA methylation and hydroxymethylation data

BS-pyrosequencing was utilized for DNA methylation validation of GZMB (cg00771752) and CHFR gene (cg00338702). For both examined loci, the Pearson correlation coefficient ( $R=$ 0.91, $p$ -value $<$ 1.87E-5 and $R=$ 0.72, $p$ -value $<$ 3.52E-7, resp.) confirmed consistency between microarray and BS-pyrosequencing results.

OxBS-pyrosequencing of 2 selected loci of MYB (cg23073132) and RNF216 (cg05617317) did not completely reproduce hydroxymethylation results obtained from EPIC array. $\mathrm{\Delta }{\beta }_{\text{BS-oxBS}}$ calculated from BS and oxBS values measured by pyrosequencing were significantly smaller than those detected by EPIC array, although the Pearson correlation coefficient between these two methods was good ( $R=$ 0.70, $p$ -value $=$ 0.024 and $R=$ 0.82, $p$ -value $=$ 0.002, resp.). This might be again caused by the generally low hydroxymethylation levels detected in blood samples compared to those previously identified in brain samples by other authors [13].

3.3 Influence of comutations (other than DNMT3A and IDH1/2) on DNA hydroxy-/methylation profile

Despite our effort to select molecularly homogenous diagnostic AML cohort, certain degree of heterogeneity could not be avoided due to the nature of AML often presenting more mutations simultaneously. In our AML mutational subgroups, NPM1 mutations were the most often occurring comutation (15/24 AML). Interestingly, only one of IDH1/2 ${}^{\text{mut}}$ samples (IDH1 ${}^{\text{mut}}$ _4) did not bear NPM1 mutation and it then clustered outside the IDH ${}^{\text{mut}}$ DNA methylation cluster. Similarly, there was a tendency for separation of DNMT3A ${}^{\text{mut}}$ samples according to the presence of NPM1 comutations. This observation is in line with previously described role of NMP1 mutations in hypermethylation phenotype [5]. Unlike NPM1 mutation, FLT3-ITD (10/24 AML) had no obvious effect on DNA methylation signature. Patients with the concurrent presence of DNMT3A and TET2 mutations ( $n=$ 2) did not display mixed DNA methylation profile as observed in DNMT3A&IDH1/2 ${}^{\text{mut}}$ samples and did not significantly differ from the rest of DNMT3A ${}^{\text{mut}}$ patients ( $n=$ 6). This might be due to the stronger effect of oncometabolite 2-hydroxyglutarate (2HG) on the activity of Tet2 protein in comparison with the effect of heterozygous TET2 mutations, where some activity of enzyme might be still preserved. Other comutations were represented less frequently and there was no apparent influence on DNA methylation profiles.

Regarding DNA hydroxymethylation, presence of TET2 mutation in two DNMT3A ${}^{\text{mut}}$ samples (DNMT3A ${}^{\text{mut}}$ _5 and DNMT3A ${}^{\text{mut}}$ _8) was reflected by their clustering next to each other. Other comutations did not apparently impact DNA hydroxymethylation patterns.

Figure 2.

GZMB methylation levels in connection with AML prognosis. A. Overall survival of AML patients based on GZMB methylation levels. B. Event-free survival of AML patients based on GZMB methylation levels.

3.4 Gene expression profiling

Hierarchical clustering of 10,067 probes (Table S4) with detection $p$ -value $<$ 0.05 did not reproduce the clusters obtained from DNA methylation or hydroxymethylation data (Fig. 1C). Only samples from healthy donors formed a consistent group. We assume that gene expression profiles (GEP) do not copy the layout of DNA hydroxy-/methylation clusters because not all DNA methylation alterations are directly translated into GEP changes. Besides, other (epi)genetic mechanisms might be involved in gene expression regulation.

3.5 Evaluation of GZMB and CHFR methylation levels as a potential prognostic biomarkers

The human granzyme B protein (GrB) encoded by GZMB gene plays a crucial role in induction of apoptosis [2]. We measured DNA methylation of GZMB enhancer (located approx. 40 kb upstream) covering two adjacent CpG sites in 104 diagnostic AML patients (for patients characteristics see Table S5). Healthy donors samples ( $n=$ 10) were examined in parallel as a control group. Based on DNA methylation levels measured in healthy controls, we set an upper normal methylation threshold of 45%. We considered a sample as unmethylated (UM) – both CpGs displayed methylation ratio $<$ 45%; intermediately methylated (IM) – one CpG $<$ 45% and one CpG $>$ 45%; and hypermethylated (HM) – both CpGs $>$ 45%. Altogether, there were 40 AML within the UM group, 15 AML in IM group and 49 AML in HM group. As GZMB encodes for a potent apoptotic inducer, we speculated whether its differential DNA methylation could be associated with altered expression on mRNA level and/or distinct prognosis. As for GZMB expression, we did not identify a clear relationship with DNA methylation levels. A weak positive or negative correlation was only observed for UM and HM samples, respectively ( $R=$ 0.3, $p$ -value $=$ 0.12 and $R=-$ 0.4, $p$ -value $=$ 0.02, resp.). However, with regard to overall (OS) and event-free survival (EFS), we found that GZMB hypermethylation was linked with inferior prognosis (Fig. 2). Using multivariate Cox regression analysis, GZMB methylation, together with age, FLT3-ITD status and transplantation in the 1 ${}^{\text{st}}$ remission, was found to be a factor influencing the prognostication model for OS. This, however, was not confirmed for EFS, where only prognostic group according to ELN [3], FLT3-ITD status and transplantation in the 1 ${}^{\text{st}}$ remission remained statistically significant (Table 3). Because GZMB hypermethylation was discovered in IDH1/2 ${}^{\text{mut}}$ samples, we tested whether the negative prognostic value of GZMB hypermethylation could be explained by presence of IDH1/2 mutations. However, only 35% of AML patients within the HM group possessed these mutations. Therefore, we assume that there might be other molecular mechanisms responsible for obtaining GZMB hypermethylation pattern, which is then related with inferior OS. Next, we tested whether the differences in DNA methylation levels are features of leukemic blasts themselves or reflections of non-leukemic cell infiltration. For 23 AML patients (from original 104 AML cohort), we obtained data on GZMB methylation also from sorted AML blast population. Reassuringly, for all samples we got the same DNA methylation category as measured in peripheral blood samples. Regarding the absence of a clear GZMB gene expression/methylation correlation and yet the influence of the particular enhancer hypermethylation on OS, we hypothesize it might be explained by either of two mechanisms. Firstly, the enhancer assigned to GZMB gene might not be in fact regulating GZMB gene, but some other biologically relevant gene for AML pathogenesis. However, GZMB seems to be the most probable gene regulated by this enhancer due to its location and orientation (it was predicted by GREAT as well as GeneHancer – database of human enhancers and their inferred target genes). Still, it cannot be excluded that other gene(s) might be under control of this locus. Secondly, it is possible that DNA hypermethylation of the discussed enhancer arises early in hematopoiesis. If this DNA hypermethylation switch occurs in some AML patients at the level of hematopoietic stem cell (HSC), it might be then transferred to both lineages – myeloid and lymphoid. In myeloid lineage, GZMB gene is not actively transcribed. Therefore, the existence or non-existence of hypermethylated enhancer does not further decrease GZMB expression. In lymphoid lineage, on the other hand, the presence of hyper- or hypomethylated enhancer might directly influence the level of GZMB expression by decreasing or enhancing the basal transcription from GZMB promoter. Hence in the lymphoid lineage, the hypermethylated GZMB enhancer can lead to decreased GZMB transcription when compared to cells without hypermethylation. This can theoretically explain the worse prognosis for patients with hypermethylated GZMB enhancer that may be coupled with decreased GZMB activity in lymphoid lineage (T-cells/NK-cells). We plan to test these hypotheses to fully unveil the mechanism behind the association of prognosis and hypermethylation of this enhancer region. We also examined gene expression and DNA methylation levels of CHFR gene. Only four percent of AML exhibited downregulation of CHFR expression by at least 2 SD (standard deviation) from average CHFR expression in healthy controls – Fig. S2. CHFR-downregulated AML had an average DNA methylation level of 93% (cg00338702), whereas the average methylation level in the rest of the samples was only 19% ( $p$ -value $<$ 0.0001). Two recent papers were devoted to CHFR expression in terms of DNA methylation. In one of them, authors described CHFR hypermethylation as a frequent event [7]. However, the location of the region assigned by the authors to CHFR was actually located in the promoter region of another gene (ZNF605). The second study claimed that DNA methylation is rare and not responsible of CHFR downregulation, which was described as a negative prognostic factor [17]. We confirmed that CHFR methylation is a rare event in AML, though we only observed CHFR downregulation in connection with hypermethylation. Hence, we could not reproduce the results regarding the role of CHFR expression downregulation as a biomarker of unfavourable prognosis in AML, because only a small proportion of AML featured such downregulation.

Table 3

Multivariate Cox regression analysis of prognostic factors in AML patients for overall and event-free survival

Variable	Multivariate analysis
	HR	95% CI	$P$ -value
A. Overall survival
GZMB methylation	1.446	1.027–2.036	0.035	*
Prognostic group according	1.612	0.952–2.729	0.076	ns
to ELN [3]
Age at diagnosis	1.036	1.008–1.065	0.012	*
Leukocytes at diagnosis	1.000	0.995–1.005	0.933	ns
Percentage of blasts	0.993	0.982–1.003	0.169	ns
FLT3-ITD	3.257	1.381–7.679	0.007	**
NPM1 mutation	0.652	0.317–1.339	0.244	ns
Transplantation in the 1 ${}^{\text{st}}$ CR	0.378	0.191–0.748	0.005	**
B. Event-free survival
GZMB methylation	1.202	0.891–1.622	0.229	ns
Prognostic group according	2.465	1.501–4.048	0.0001	***
to ELN [3]
Age at diagnosis	1.013	0.991–1.036	0.246	ns
Leukocytes at diagnosis	1.001	0.996–1.006	0.604	ns
Percentage of blasts	1.000	0.991–1.010	0.960	ns
FLT3-ITD	2.696	1.273–5.711	0.01	*
NPM1 mutation	0.737	0.385–1.411	0.356	ns
Transplantation in the 1 ${}^{\text{st}}$ CR	0.178	0.088–0.360	0.0001	***

Prognostic variables and corresponding hazard ratios (HR), confidence intervals (CI) and $P$ -values for multivariate analysis in terms of overall survival (A) and event-free survival (B). ns – not significant, CR – complete remission.

4. Conclusions

In this study, we found out that specific mutational background in AML patients is connected with distinct DNA hydroxy-/methylation patterns based on presence of DNMT3A ${}^{\text{mut}}$ or IDH1/2 ${}^{\text{mut}}$ . Furthermore, mixed DNA hydroxy-/methylation profile for the concurrent presence of these mutations was observed. DNA hydroxymethylation was low in numbers and levels of hmCpGs proposing that it is less important in AML pathogenesis when compared to DNA methylation involvement. To the best of our knowledge, this is the first report providing a list of genomic regions affected by aberrant hydroxymethylation with single-nucleotide resolution in AML. We also identified profound difference in numbers of sites affected by DNA methylation and hydroxymethylation in IDH1 ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ patients suggesting different biological background of these mutations. In terms of biologically relevant DNA methylation changes, hypermethylation of enhancer located upstream GZMB gene, encoding an apoptotic inducer, was found to be connected with inferior overall survival. Further detailed studies clarifying the role of DNA hydroxymethylation in AML and reasons for DNA hydroxy-/methylation variances in IDH1 ${}^{\text{mut}}$ and IDH2 ${}^{\text{mut}}$ are needed.