DNA methylation profile in chronic myelomonocytic leukemia associates with distinct clinical, biological and genetic features

ABSTRACT

Chromosomal abnormalities are detected in 20–30% of patients with chronic myelomonocytic leukemia (CMML) and correlate with prognosis. On the mutation level, disruptive alterations are particularly frequent in chromatin regulatory genes. However, little is known about the consequential alterations in the epigenetic marking of the genome. Here, we report the analysis of genomic DNA methylation patterns of 64 CMML patients and 10 healthy controls, using a DNA methylation microarray focused on promoter regions. Differential methylation analysis between patients and controls allowed us to identify abnormalities in DNA methylation, including hypermethylation of specific genes and large genome regions with aberrant DNA methylation. Unsupervised hierarchical cluster analysis identified two main clusters that associated with the clinical, biological, and genetic features of patients. Group 1 was enriched in patients with adverse clinical and biological characteristics and poorer overall and progression-free survival. In addition, significant differences in DNA methylation were observed between patients with low risk and intermediate/high risk karyotypes and between TET2 mutant and wild type patients. Taken together, our results demonstrate that altered DNA methylation patterns reflect the CMML disease state and allow to identify patient groups with distinct clinical features.

KEYWORDS:

• Chronic myelomonocytic leukemia
• DNA methylation
• hypermethylation
• prognosis
• TET2

Introduction

Chronic myelomonocytic leukemia (CMML) is a rare clonal hematopoietic malignancy with features from both myelodysplastic syndromes (MDS) and myeloproliferative neoplasms (MPN) [1]. The 2008 World Health Organization (WHO) classification distinguishes two categories (CMML-1 and CMML-2) according to blast percentage in bone marrow (BM) or peripheral blood (PB), while the original French-American-British (FAB) Group distinguishes two disease variants, myelodysplastic and myeloproliferative CMML, based on leukocyte count [MD-CMML (<13 × 10⁹/L) and MP-CMML (>13 × 10⁹/L)) [1,2]. Clinical course in CMML is variable, with a median overall survival (OS) of 20 months and progression to acute myeloid leukemia (AML) in 15–30% of cases [1].

Large scale chromosomal abnormalities are detected in 20–30% of patients at diagnosis and correlate with clinical outcomes [3,4]. Most frequent alterations recur in 3–8% of cases and include trisomy 8 (+8), abnormalities of chromosome 7 (del(7q)/-7), loss of Y chromosome (-Y), or an overall complex karyotype [3,4]. According to the CMML-specific cytogenetic risk classification, proposed by the Spanish Group of MDS, up to 80% of CMML patients fall into the low risk cytogenetic category, characterized by normal karyotype and isolated -Y. Patients with complex karyotype, +8 or del(7q)/-7 are classified into the high risk group, while the intermediate risk category includes all other alterations [3].

Mutations in epigenetic regulatory genes are present in more than 60% of CMML patients and include disruptive mutations in genes regulating DNA methylation, such as in TET2 (50-60%) and DNMT3A (5-10%), and post-translational histone modification, such as in ASXL1 (40-50%) and EZH2 (5-10%) [5–8]. In cancer, there are complex changes in DNA methylation [9]. Frequently, global DNA demethylation is accompanied by the de novo methylation of CpG islands. In particular, hypermethylation of CpG islands overlapping gene promoters has been associated with the transcriptional silencing of tumor suppressor genes (TSG), leading to clonal advantage of the affected cell [10]. DNA methylation has been widely investigated in hematological neoplasms [10–12]. MDS is characterized by local DNA hypermethylation that increases with disease severity [13–15]. In AML, DNA methylation correlates with specific cytogenetic abnormalities [14] and the hypermethylation of specific genes such as CDH1, CDKN2B and DBC1 has been shown to have an impact on the outcome of MDS and AML patients [13–15]. In CMML, hypermethylation of genes such as CDKN2B and TIF1γ has been reported [16,17]. In addition, a signature of DNA methylation was found to predict the response of CMML patients to hypomethylating agents [18]. It is further clear that TET2 mutations correlate with changes in DNA methylation in CMML and other diseases [19–22]. These studies show that DNA methylation plays a role in the pathogenesis of CMML. However, little is known about genome-wide DNA methylation changes and their relevance for disease progression.

In this study, we performed DNA methylation microarrays in 64 CMML patients with the aim to investigate the differences in the methylation profile between CMML patients and healthy controls, as well as to evaluate the possible prognostic impact of DNA methylation abnormalities. In addition, we have explored the differences in DNA methylation between patients with different cytogenetic and molecular features.

Results

Characteristics of CMML patients

A total of 64 patients with CMML were included in the study. Median follow-up of living patients was 36 months (range: 5.7 to 83.3 months). Main clinical and biological characteristics of patients are summarized in Table 1. Median age at diagnosis was 71 years and there was a 1.8:1 male:female ratio. According to the FAB criteria [2], 47 (73%) patients were classified as MD-CMML and 17 (27%) as MP-CMML. Considering the 2008 WHO classification [1], 46 (72%) cases corresponded to CMML-1 and 18 (28%) to CMML-2. Progression to AML was observed in 25 (39%) patients. Conventional G-banding cytogenetics (CC) was performed at diagnosis and an informative result was obtained in all but one patient. This patient had no chromosomal abnormalities when analyzed by SNP array (SNP-A). According to the CMML-specific cytogenetic risk stratification [3], most patients (50, 78%) had low risk cytogenetic features: normal karyotype (47) and isolated -Y (3). The rest of cases displayed alterations that corresponded to intermediate (7, 11%) or high risk (7, 11%) groups (Table 1). Targeted deep sequencing of 83 myeloid-related genes had been previously performed on a subset of these patients [23]. The following genes were mutated in >10% of analyzed cases: TET2 (34/46, 73%), ASXL1 (26/57, 46%), SRSF2 (21/57, 37%), RUNX1 (13/46, 28%), CBL (8/57, 14%), ZRSR2 (6/46, 13%) and EZH2 (5/46, 11%) (Table 1).

Table 1. Main clinical and hematological characteristics of CMML patients at diagnosis (n = 64).

Variable	Median (range)	n (%)
Age, years	71 (41–89)
<70		25/64 (39)
≥70		39/64 (61)
Gender
Male		41/64 (64)
Female		23/64 (36)
FAB classification
Myelodysplastic (CMML-MD)		47/64 (73)
Myeloproliferative (CMML-MP)		17/64 (27)
WHO classification
CMML-1		46/67 (72)
CMML-2		18/64 (28)
Hemoglobin level, g/dL	11.7 (6.5–15.2)
<10 g/dL		16/64 (25)
≥10 g/dL		48/64 (75)
Leukocyte count, x10^9/L	8.7 (3.2–88.0)
<13 × 10^9/L		47/64 (73)
≥13 × 10^9/L		17/64 (27)
Platelet count, x10^9/L	135.0 (5.6–481.0)
<100 × 10^9/L		19/64 (30)
≥100 × 10^9/L		45/64 (70)
Neutrophil count, x10^9/L	4.3 (0.5–61.6)
<1.8 × 10^9/L		9/60 (15)
≥1.8 × 10^9/L		51/60 (85)
Blasts in BM, %	3.0 (0.0–19.0)
<10%		57/64 (89)
≥10%		7/64 (11)
Progression to AML
Yes		25/64 (39)
No		39/64 (61)
Cytogenetics
Low risk karyotype		50/64 (78)
Normal karyotype		47/64
Isolated -Y		3/64
Intermediate risk karyotype		7/64 (11)
High risk karyotype		7/64 (11)
Gene mutations
TET2		34/46 (73)
ASXL1		26/57 (46)
SRSF2		21/57 (37)
RUNX1		13/46 (28)
CBL		8/57 (14)
ZRSR2		6/46 (13)
EZH2		5/46 (11)
CPSS risk group [53]
Low		32/64 (50)
Intermediate-1		13/64 (20)
Intermediate-2		14/64 (22)
High		5/64 (8)

AML: acute myeloid leukemia; BM: bone marrow; CPSS: CMML-specific prognostic scoring system; FAB: French-American-British; Uninformative CC: cases with no metaphases; WHO: World Health Organization.

Methylation profiling of CMML

DNA methylation was assessed with the Infinium HumanMethylation450 BeadArray). All samples passed quality control criteria. After preprocessing, a total of 48,186 probes, which overlapped with SNPs (47,007 sites) or corresponded to unreliable measurements (1,179 sites), were removed from further analyses. Unsupervised analysis of the resulting methylation profiles was performed by principal component analysis (PCA), which segregated the malignant samples from the control samples (Supplementary Figure 1A). Based on methylation data, significant (P<0.01) association between the following clinical traits was observed: CMML-specific prognostic scoring system (CPSS) was associated with FAB subtype, WHO category, cytogenetic features, and AML progression; FAB and WHO classifications were also correlated with AML progression (Supplementary Figure 1B).

Figure 1. Hierarchical clustering according to DNA methylation in CMML distinguishes two main clusters. Heatmap showing the most variable (>0.4, P < 0.05) differentially methylated probes across the series of CMML patients. Each row represents a single probe and each column corresponds to a separate patient. Each probe has been normalized by the mean β-value in all samples and the color in the heatmap represents the scaled methylation value. Below is represented in color-code the distribution of the clinical features from the patients (WHO classification, FAB classification, AML progression, CPSS risk groups, and cytogenetic features and gene mutations). White gaps represent missing data. AML: acute myeloid leukemia; CPSS: CMML specific prognostic scoring system; FAB: French-American-British; Int: intermediate; WHO: World Health Organization.

Differential methylation between CMML and healthy controls

To identify changes in DNA methylation associated with CMML, we compared the DNA methylation patterns between CMML patients (n = 64) and healthy controls (n = 10) by differential methylation analysis. Differentially methylated probes included 1,844 probes hypermethylated and 255 probes hypomethylated in CMML patients compared to controls (>0.3 methylation difference, P<0.05). Hypermethylated probes affected a total of 1,100 genes and 636 promoter or enhancer regions, while hypomethylated probes affected 165 genes. To investigate whether hypermethylation in CMML patients was associated with specific biological processes, we performed Gene Ontology (GO) analysis of all the hypermethylated genes, but the analysis did not reveal a significant enrichment in any GO term. Differentially methylated blocks (DMBs) analysis across the genome (see Materials and Methods) revealed that aberrant methylation was not distributed randomly across chromosomes. Bin screening using the 10 Mb filter revealed 14 DMBs in which DNA methylation was significantly altered across the whole region (Supplementary Table 1). The region in which DNA methylation showed the highest significant difference, using both 10 Mb and 1 Mb filters, was located in chromosome 5q31.3 and included the protocadherin (PCDH) gene cluster. This cluster, comprising 53 genes and 3 pseudogenes, expands across approximately 1 Mb and contains the three alpha, beta, and gamma PCDH clusters tandemly linked on this region (Supplementary Figure 2). Other highly significant DMBs included a region on chromosome 17q24.3-q25.3 that contained genes with several altered sites, such as P4HB, RPTOR, and TBC1D16, as well as a region on chromosome 6p22.1-p21.2, including DAXX, DDAH2, PPP1R18, PRRT1, TNF, TRIM39, and ZBTB12. We also identified several genes (72) aberrantly methylated in CMML patients that play a role in the WNT signaling pathway, including APC2, CDH1, CDH13, CDH23, CREBBP, CSNK1A1, CSNK1D, FBXW11, LRG4, PLCB3, WNT5, and several PCDH genes (Supplementary Table 2). Aberrant DNA methylation was also detected in many genes involved in inflammation mediated by chemokine and cytokine signaling pathways (e.g., CCR1, CCR2, CCR3, CCRL2, CX3XR1, CXCR1, CXCR2, IFNGR2, PIK3D, PRKCZ), PDGF signaling (e.g., ETS1, HRAS, MAPK3, PDGFD, PDGFRB, PIK3CD, PIK3R1, PRKCA, SHC1, SRC), and apoptosis (e.g., BCL2A1, BIRC3, CRADD, DAXX, FADD, TNF, TNFRSF1A, TNFRSF1B, TNFSF10, TRAF2).

Figure 2. Kaplan-Meier survival curves in CMML patients. (A) Overall survival (OS) and progression free survival (PFS) curves according to the two groups identified by hierarchical clustering. (B) OS and PFS curves according to BCL2 DNA methylation levels.

Table 2. Comparison of the main clinical and biological parameters of CMML patients between the two groups identified by hierarchical cluster analysis.

Variable		Group 1 n (%)	Group 2 n (%)	P value
Sex	Male	17 (68)	24 (62)	0.790
	Female	8 (32)	15 (39)
Age (years)	<70 years	15 (60)	10 (26)	0.009
	≥70 years	10 (40)	29 (74)
WHO classification	CMML-1	11 (44)	35 (90)	<0.001
	CMML-2	14 (56)	4 (10)
FAB classification	MD-CMML	16 (64)	31 (80)	0.246
	MP-CMML	9 (36)	8 (20)
Hemoglobin level	<10 g/dL	10 (40)	6 (15)	0.039
	≥10 g/dL	15 (60)	33 (85)
Leukocyte count	<13 × 10^9/L	16 (64)	31 (80)	0.246
	≥13 × 10^9/L	9 (36)	8 (20)
Platelet count	<100 × 10^9/L	3 (12)	16 (41)	0.023
	≥100 × 10^9/L	28 (88)	23 (59)
Neutrophil count	<1,8 × 10^9/L	4 (16)	9 (23)	0.544
	≥1,8 × 10^9/L	21 (84)	30 (77)
Blasts in BM	<10%	20 (80)	37 (95)	0.100
	≥10%	5 (20)	2 (5)
Cytogenetic risk	Low risk	15 (60)	35 (90)	0.011
	Int/high risk	10 (40)	4 (10)
CPSS	Low	7 (28)	25 (64)	<0.001
	Intermediate-1	5 (20)	8 (21)
	Intermediate-2	8 (32)	6 (15)
	High	5 (20)	0 (0)
AML progression	Yes	14 (56)	11 (28)	0.036
	No	11 (44)	28 (72)
TET2	Wild type	9 (69)	4 (12)	<0.001
	Mutated	4 (31)	30 (88)
ASXL1	Wild type	7 (35)	24 (65)	0.100
	Mutated	13 (65)	13 (35)
SRSF2	Wild type	17 (85)	19 (51)	0.020
	Mutated	3 (15)	18 (49)
RUNX1	Wild type	11 (85)	23 (68)	0.301
	Mutated	2 (15)	11 (32)
CBL	Wild type	48 (90)	31 (84)
	Mutated	2 (10)	6 (16)	0.699
ZRSR2	Wild type	8 (62)	33 (97)	0.004
	Mutated	5 (39)	1 (3)
EZH2	Wild type	12 (92)	30 (88)	1.000
	Mutated	1 (8)	4 (12)

Variables with significant differences between the two groups (P<0.05) are highlighted in bold. AML: acute myeloid leukemia; BM: bone marrow; CPSS: CMML specific prognostic scoring system; FAB: French-American-British; Int: intermediate; WHO: World Health Organization.

Probes differentially methylated in CMML identify two clusters with clinical relevance

To explore associations among CMML samples, 2-dimensional hierarchical cluster analysis was performed, considering all CMML patients, using the most variable probes that displayed an absolute difference >0.4 (n = 333; P<0.05) between patients and controls (Supplementary Table 3). The analysis revealed 3 main clusters, one with the 10 control samples plus two patients, and then two other clusters with 23 and 39 patients, respectively (Figure 1). According to this methylation profile, we grouped the patients into two groups and then analyzed the differences between their biological and clinical characteristics (Figure 1, Table 2). Statistically significant differences were observed in the following parameters between Group 1 (n = 25) and Group 2 (n = 39): age, WHO classification, AML progression, cytogenetic risk groups, hemoglobin (Hb) levels, platelet count, CPSS groups, and SRSF2 mutations (Table 2). Group 1 was enriched in patients with advanced disease features, such as CMML-2, presence of anemia or thrombocytopenia, higher risk CPSS groups, and intermediate or high risk karyotypes. Therefore, we performed survival analysis, which revealed significant differences regarding the outcome of the patients in the two groups. Group 1 had an inferior overall survival (OS) and progression-free survival (PFS) compared to Group 2 [median OS of 13 months, confidence interval 95% (CI95%) 10–15, for Group 1 vs. median OS of 55 months, CI95% 13–96, for Group 2, P = 0.014; median PFS of 24 months, CI95% 20–27, for Group 1 vs. median PFS not reached, for Group 2, P = 0.017] (Figure 2A, Supplementary Table 5). Regarding molecular data, significant differences were observed in the following genes: ZRSR2 mutations were restricted to Group 1, while Group 2 was enriched in patients with mutations in TET2 and SRSF2 genes (Table 2, Figure 1).

Gene expression and prognostic value of the DNA methylation of specific genes

Due to the fact that assessing genome-wide methylation profiles has a limited applicability in a diagnostic setting, the methylation level of individual selected gene promoters, which, according to previous studies, may play a role in MDS and myeloid-related diseases, was correlated with OS and PFS: CDH1, ETS1, DAXX, FADD, DAPP1, AATK, CYFIP, TP53BP2, and AIM2, which were hypermethylated in CMML samples, and BCL2, which was hypomethylated. Only hypomethylation of BCL2 was predictive of poorer OS and PFS. Patients with demethylation of BCL2 had an inferior OS and PFS compared to patients that displayed BCL2 methylation levels comparable to the ones observed in healthy controls (median OS of 14 months, CI95% 6–21, for hypomethylation vs. median OS of 55 months, IC95% 9–100, for normal methylation, P = 0.002; median PFS of 23 months, CI95% 13–33, for hypomethylation vs. median PFS not reached, for normal methylation, P = 0.001) (Figure 2B, Supplementary Table 4, Supplementary Table 5). In order to assess if increased or reduced methylation levels of these genes resulted in gene expression changes, we analyzed their expression by RT-PCR in 20 CMML patients and 4 healthy controls (Supplementary Table 6). Regarding hypermethylated genes, gene expression of AATK, AIM2, CYFIP2, DAPP1, and FADD was significantly decreased in patients compared to controls. Conversely, CDH1 was significantly overexpressed in CMML samples, despite being hypermethylated, while no significant changes in gene expression were observed for genes DAXX, ETS1, and TP53BP2. On the other hand, BCL2 hypomethylation did not result in increased expression of the gene.

Multivariate survival analysis

In order to assess the independent prognostic role of DNA methylation data, we performed an adjusted multivariate survival analysis, including also clinical, biological, and genetic features that were statistically significant in the univariate OS and PFS analyses (Supplementary Table 5). None of the variables related to DNA methylation remained significant in these analyses. Independent prognostic variables for OS were MP-CMML FAB subtype/leukocyte count ≥13 × 10⁹/L, hemoglobin level >10 g/dL, and EZH2 mutations; for PFS, hemoglobin level >10 g/dL and presence of splenomegaly (Supplementary Table 6).

Differential methylation in patients with distinct cytogenetic and molecular features

Finally, differential methylation analysis was performed between patients that presented with distinct cytogenetic and molecular features (Figure 3). The most significant differences were observed between patients with different risk karyotypes and between patients with different mutational status of TET2 and EZH2, so we focused on these comparisons.

Figure 3. Differential methylation between CMML patients displaying distinct molecular and cytogenetic features. The bar plot shows the number of probes differentially methylated (P<0.05) between groups of patients characterized by different molecular features (presence vs. absence of specific gene mutations) and cytogenetic [intermediate/high risk karyotype (IHR-K) vs. low risk karyotype (LR-K)].

Higher methylation levels were detected in patients with intermediate/high risk karyotype (IHR-K) compared to the low risk karyotype (LR-K) group (Supplementary Figure 3A). A total of 240 and 12 probes were significantly hypermethylated and hypomethylated, respectively, in the group with IHR-K (>0.3 methylation difference, P<0.05). Interestingly, 70% (167/240) of the hypermethylated probes corresponded to CpGs located at promoter or enhancer regions. Overall, 122 different genes were significantly more methylated in patients with IHR-K and only 3 genes were less methylated. The combination of differential methylation analysis and DMBs analysis across the genome allowed us to identify regions and genes with significantly different levels of DNA methylation between these two groups and revealed an enrichment in transcription factor proteins. The region displaying the highest difference between the two groups of patients corresponded to WT1 gene, which was more methylated in patients with IHR-K (Supplementary Figure 4A). Another region identified corresponded to the homeobox HOXB gene cluster, located at chromosome 17q21.32 and consisting of 9 genes (Supplementary Figure 4B). Other hypermethylated transcription factors included PAX8, NR5A2, LMO3, PRDM16, THRB, EVX2, TFAP2B, TFAP2E, PDLIM4, and other homeobox genes (HOXA2, HOXA10, HOXC4, HOXC9, HOXC12, HOXD9, HOXD10).

Then, we focused on the differences between TET2 mutant and wild type patients. Cases with TET2 mutations were characterized by significantly higher methylation levels (Supplementary Figure 3B), with 667 probes hypermethylated and only 9 probes hypomethylated, affecting 518 and 8 genes, respectively (>0.3 methylation difference, P<0.05). GO analysis of the differentially methylated promoter regions revealed a significant enrichment in several GO terms, including immune response, cell activation and multiple hematological processes, such as leukocyte activation, proliferation or migration, monocyte and granulocyte chemotaxis, granulocyte activation, myeloid leukocyte migration, and myeloid leukocyte mediated immunity, among others (Supplementary Table 7). Regarding EZH2, although a high number of probes was differentially methylated between mutant and wild type patients (Figure 3), these differences were more subtle and only 44 and 46 probes were hyper- and hypo-methylated, respectively, when we applied the same cut-off value (>0.3 methylation difference, P<0.05). In contrast to TET2 mutants and patients with IHR-K, EZH2 mutations were not associated with global higher methylation levels (Supplementary Figure 3C).

Discussion

Cytogenetic abnormalities and recurrent genetic lesions may be insufficient to explain the biological diversity of CMML. The initiation and development of cancer involves both genetic and epigenetic alterations and increasing evidence shows that hypermethylation of genes occurs in myeloid malignancies and potentially contributes to the pathogenesis of these diseases [10]. Nevertheless, the role of DNA methylation in CMML remains unclear and little is known about the relevance of genome-wide methylation changes in the disease progression. In the present study, we performed an integrative genome-wide DNA methylation profile analysis in 64 CMML patients using the Infinium HumanMethylation450 BeadArray that interrogates the DNA methylation status of 450,000 CpGs throughout the genome. It should be noted that, while the Infinium HumanMethylation450 BeadArray used in this and many other studies covered much of the genome, including 99% of RefSeq genes, it is mostly focused on promoter CpG islands and their shores, and thus the effect of CpGs located at promoter-distal and intergenic regulatory regions could not be assessed.

Our results show significant higher levels of DNA methylation in >1,800 CpGs in patients compared to controls. Aberrant DNA methylation of promoters in CMML and MDS has previously been reported in specific genes [16,24–27]. Here, with the use of a wider approach, we found significant and widespread promoter hypermethylation in CMML, which affected 636 regulatory regions, including 296 promoter sites in 244 different genes. As promoter hypermethylation affected many known TSGs, such as BCL2, BRCA2, PRDM2, SEPT9, and VHL, likely reflecting their transcriptional inactivation, it is likely that this contributed to positive selection in neoplastic evolution. The extent of the hypermethylation further suggests that it also includes many passenger alterations, or so called epiphenomenon effect, that are not disease causative [26]. Of note, epigenetic deregulation was not limited to cancer-associated genes but appeared to be a more widespread phenomenon. Aberrant methylation was detected in several genes involved in the WNT signaling pathway, similarly to what has been reported in MDS, secondary AML [26,28], and several solid tumors [29–31]. Among these genes, aberrant methylation of CDH1 has been consistently reported in MDS patients and has been correlated with poor prognosis [15,25,32]. Activation of the WNT signaling pathway has been associated to cancer initiation, and WNT antagonists or inhibitors are known TSGs, suggesting that hypermethylation-induced silencing of these genes may play a role in the pathogenesis or progression of these malignancies [33]. Many genes involved in inflammation mediated by chemokine and cytokine signaling pathways, PDGF signaling and apoptosis were also affected by aberrant DNA methylation, some of which have been described in other myeloid malignancies, such as MDS and AML. Methylation-induced silencing of the apoptotic gene ETS1 was also reported by del Rey et al. in low risk MDS patients [34]. Conversely, del Rey et al. also reported methylation-induced silencing of BCL2 in low risk MDS patients, while in our cohort, BCL2 promoter was methylated in healthy controls and hypomethylated in a subset of CMML patients. In fact, a significant decrease in OS and PFS was observed in this last group. However, we did not observe a significant increase in BCL2 expression in CMML patients compared to controls, suggesting that it may be regulated by other mechanisms. We also detected DNA hypermethylation in several genes previously reported in MDS and AML, including BCL11A, BRCA2, CYFIP2, DAPP1, DAXX, FADD, KLF11, LCK, LRRF1P1, NOTCH4, PDCHB3, PDCHB7, SMPD3, TGFA, TNF, and ZFP161 [13,14,28,35,37]. Furthermore, we observed that hypermethylation of some of these genes (AATK, CYFIP2, DAPP1, FADD) led to a significant decrease in gene expression. All these results suggest that some of these genes may be involved in pathogenic steps that are common across the spectrum of myeloid malignancies and that may be a necessary, but not sufficient, event during the malignant transformation process of hematopoietic cells. However, specific patterns and distributions of epigenetic deregulation may characterize different forms of hematologic malignancies, potentially contributing to the biology of each disease [28]. Identifying targets for altered DNA methylation conserved across all myeloid neoplasms will help in the identification of key drivers of these diseases. The identification of methylation changes in common pathways provides the rationale for testing the inhibition of key signaling components of these pathways for novel therapeutic intervention strategies.

In our study, the effect of DNA hypermethylation was not completely stochastic and certain chromosomal regions presented more hypermethylation than others. Accordingly, it has been shown that DNA hypermethylation can extend over regions of the genome from 100 to >1,000 kb in size and is linked to long-range epigenetic silencing [38]. Identified by DMBs analysis, the region with most significant differences in DNA methylation between patients and controls was located at 5q31.3 and included the PCDH gene family, which is formed by 3 clusters of PCDH genes (α, β, and γ) arranged in tandem that codify for a superfamily of cell adhesion molecules. The expression of PCDH genes by hematopoietic stem cells could play an important role in the hematopoietic niche [39]. The genes within each α and γ cluster share the same last three 3′ exons, while genes in the β cluster are distinct and consist of only one exon [40]. Each PCDH gene in these three gene clusters also has a CpG island promoter. Global methylation silencing of clustered PCDH genes has been previously described in solid tumors [38,41]. Another significant region detected by this analysis was 17q24.3-q25.3. Interestingly, Bullinger et al. identified a CpG island-dense region on 17q25.3 that was predictive of poor outcome in AML patients [13].

Hierarchical cluster analysis of CMML samples using the most variable probes between patients and controls revealed two main clusters: Group 1 and Group 2. Group 1 was enriched in patients with advanced disease features, such as CMML-2, presence of anemia or thrombocytopenia, intermediate/high risk karyotypes and higher risk CPSS groups. Furthermore, patients from Group 1 had an inferior OS and PFS compared to Group 2, although these differences were not statistically significant in the multivariate analyses when clinical, biological, and genetic data were included. These results suggest that DNA methylation alteration has a major impact on the biology of the disease and that these altered patterns of DNA methylation may be an important determinant of the different clinical behavior of CMML patients. Moreover, differences in DNA methylation patterns at diagnosis may be useful biomarkers for prediction of clinical outcome.

The two identified clusters also correlated with cytogenetic and molecular characteristics. Group 2 was enriched in patients with normal karyotype. Correlation between genome-wide DNA methylation patterns and cytogenetic features has also been described in both AML and ALL, although clustering is never absolute [14,42]. Patients with intermediate/high risk karyotypes were characterized by higher levels of DNA methylation, supporting that aberrant methylation is a progressive process that increases in more advanced states of the disease, as previously hypothesized [26]. Another possibility is that aberrant DNA methylation could in some way be triggered by chromosome instability [28]. Comparison between patients with intermediate/high risk karyotypes (IHR-K) vs. low risk karyotypes (LR-K) revealed an abnormal methylation in several transcription factors in the IHR-K group. Both WT1 and its antisense WT1-AS1 were significantly hypermethylated in cases with IHR-K. WT1 can function as an oncogene or a TSG and plays a role in the pathogenesis of several hematological neoplasms [43]. Promoter methylation of WT1 and WT1-AS has been observed in hematological neoplasms and solid tumors [44–47]. Another differentially methylated region was also detected at 17q21.2-q21.33, which includes the HOXB gene cluster that codifies for a family of transcription factors that have been shown to play a role in both myeloid and lymphoid malignancies [48–51].

Regarding molecular data, we found that ZRSR2 mutations were restricted to Group 1, while Group 2 was enriched in patients with mutations in TET2 and SRSF2 genes. This interesting observation is limited by the fact that mutation data was only available for a subset of patients of which more fell into Group 2 than 1. Differences in DNA methylation between patients with different molecular features were also explored. As expected, since TET2 enzyme is involved in the methylation cycle, the highest number of differentially methylated probes was observed between patients with TET2 mutant and wild type cases. Differential methylation analysis revealed higher methylation levels in patients with TET2 mutations, in accordance with previous studies [19,21]. Focusing on specific genes, Yamazaki et al. reported hypermethylation of AIM2 and SP140 promoters in TET2 mutant CMML patients [21], both of which were also highly methylated in TET2 mutant cases in our cohort. Of note, AIM2 hypermethylation correlated with a significant decrease in gene expression. The role of this gene in hematological neoplasms is unknown, but it has been reported that decrease of AIM2 contributes to hepatocellular carcinoma metastasis [51]. Similarly, Pérez et al. detected hypermethylation in these two genes and other 11 genes in TET2 mutant CMML patients, 9 of which were also more methylated in our series (CCIN, CD48, GPR109B, LAX1, MS4A1, NKG7, SLC22A12, TRAF1, and VHL) [52]. Compared to these studies, our genome-wide approach allowed us to identify a total of 518 genes hypermethylated in TET2 patients. Of note, these were enriched in several GO terms related to hematological processes. While we have used a CpG island-focused technique, a recent study has pointed out that TET2 mutations might have a more pronounced effect on DNA methylation outside of CpG islands [20], which we could not assess due to the nature of the technique used in the present study.

In conclusion, our data show that epigenetic lesions are abundant in CMML, suggesting that some driver lesions in CMML could be of epigenetic nature. In addition, the clinical significance of DNA methylation patterns is underlined by the fact that it contributes to identifying groups of patients that share common biological and clinical features. Finally, some of these alterations could be targets for the development of novel therapeutic approaches.

Materials and methods

Patients and samples

A total of 64 patients with CMML and 10 controls (healthy hematopoietic stem cell donors) were included in the study. Patients were diagnosed at the Institut Català d'Oncologia (ICO) following the 2008 WHO criteria [1]. Patients were also classified according to the FAB classification and the prognostic score for each patient was calculated according to the CPSS [2,53]. The study was approved by the ICO-Hospital Germans Trias i Pujol Ethics Committee and an informed consent was obtained from each patient in accordance with the Declaration of Helsinki.

Samples

BM samples from CMML patients were collected at diagnosis in all cases and total BM cells were obtained after erythrocyte lysis. Total cells from CD34+ enriched mobilized PB from 10 healthy donors were obtained. Mobilized CD34+ cells were obtained from PB after donors’ treatment with G-CSF for five days. Genomic DNA was extracted in all cases with QiaAmp DNA Blood Mini kit (Qiagen, Hilden, Germany) and quantified using Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, CA, USA). Total RNA was extracted from RNA samples in 20 CMML patients and 4 controls using the QiaAmp RNA Blood Mini kit (Qiagen, Hilden, Germany). The purity and quantity of DNA and RNA was evaluated by spectrophotometry, using the Nano Drop 2000 (Thermofisher Scientific, Wilmington, DE).

Cytogenetics

Conventional G-banding cytogenetics was performed on BM samples at diagnosis following standard procedures [3]. Karyotypes were described according to the International System for Human Cytogenetic Nomenclature [54]. In one uninformative case with no metaphases, SNP-A analysis was performed as previously described [55].

Targeted deep sequencing

Targeted deep sequencing of a panel consisting of 83 myeloid-related genes had been previously performed on a subset of patients (n = 46), as previously published [23]. Libraries were prepared and enriched using the Kapa Library Preparation Kit (Kapa Biosystems, MA, USA) and the SeqCap EZ capture chemistry (Nimblegen, Roche, Basel, Switzerland), and then sequenced in an Illumina MiSeq following a 150 bp paired-end reads standard protocol. Sequencing data was analyzed using an in-house bioinformatics pipeline [23]. Briefly, reads were aligned against human genome build 19 (hg19) using BWA 0.7.12. Post-alignment including local indel realignment and base recalibration was performed using the tools in GATK 3.4.46 software package. Packages SAMtools 1.2 and VarScan 2.4.0 were used for variant calling and ANNOVAR for variant annotation. High-probability oncogenic mutations were called by eliminating sequencing and mapping errors and by discarding variants located in highly variable regions or with low coverage (<100x), as well as synonymous variants and polymorphisms described on the available databases.

DNA methylation microarrays

DNA methylation was assessed with the Infinium HumanMethylation450 BeadArray (Illumina, San Diego, CA, USA). First, sodium bisulfite conversion was performed on 1μg of DNA using the EZ DNA Methylation-Gold™ Kit (Zymo Research, Orange, CA, USA) following the manufacturer's protocol. Bisulfite converted DNA was then processed and hybridized to the array according to the manufacturer's recommendations. The HumanMethylation450 DNA Methylation array assesses the methylation status of 485,764 cytosines distributed over the whole genome, which correspond to CpG dinucleotides (99.3%) and CNG targets (0.7%) [56]. According to their associated RNA transcripts, 361,766 CpGs (74.4%) correspond to classic coding messenger RNA genes, 4,168 (0.85%) are linked to non-coding RNAs (microRNAs and long non-coding RNAs) and for 119,830 (24.6%) sites there are no annotated transcripts associated with the described CpG location [56]. Fluorescence signals generated for unmethylated and methylated cytosines were transformed into beta-value (β-value) ratio ranging from 0 (completely unmethylated site) to 1 (completely methylated site). Raw data are available at the Gene Expression Omnibus (GEO) database under the accession number GSE105420. Raw IDAT files from Infinium HumanMethylation450 experiments were imported into RnBeads for DNA methylation analysis, including quality control, data preprocessing and normalization [57]. The methylation β-values were normalized using the BMIQ normalization method [58]. RnBead was used to perform unsupervised principal component analysis (PCA) and to explore associations between traits [57]. Differential methylation analysis was performed to identify hypermethylated and hypomethylated probes between the compared groups. Aberrant methylation was defined as a mean absolute difference >0.3 (adjusted P<0.05) between the average β-values of the two groups, for each specific site. To calculate DMBs, for each comparison, the genome was divided in bins (0.1, 1, 10 Mb). For each bin, a score was calculated dividing the probes differentially methylated (we set the cut-off for the differential methylation using all the probes that showed an adjusted P≤0.01 and an absolute differential methylation >0.1), with the total probes present in the bin and then multiplied for the probes differentially methylated. The bins with the higher rank, that showed a score belonging at the higher 5% of the calculated scores, were further investigated. Hierarchical clustering was performed using the function hclust, with default parameters, from the Bioconductor Heat Plus package v.2.22. All statistical analyses were performed using RnBeads package version 1.8.0. Post-statistical analyses were done using in-house R scripts.

Gene expression analysis

A total of 150 ng of RNA were retrotranscribed using the Reverse Transcription Master Mix and cDNA products were pre-amplified with the PreAmp Master Mix (Fluidigm Co., San Francisco, CA, USA). In all cases protocols supplied by the manufacturer were followed. DELTAgene® assays (Fluidigm Co.) were designed for 10 selected genes (CDH1, BCL2, ETS1, DAXX, MFADD, DAPP1, AATK, CYFIP, TP53BP2, and AIM2) plus two housekeeping genes (GAPBH and GUSB) (Supplementary Table 8). Gene expression was assessed by Real-Time PCR using the microfluidic system Biomark HD and the FLEXsix™ Gene Expression IFCs, following the manufacturer's recommendations (Fluidigm Co.). Each measurement was performed in triplicate. The relative expression levels between CMML patients and healthy controls, normalized to GAPDH housekeeping gene, was determined using the 2^−ΔΔCt method [59].

Statistical and survival analyses

Baseline demographic and clinical–biological characteristics were described as median and range for continuous variables, while frequency and percentage were used for categorical variables. Comparisons of variables between patient subsets were performed by χ², Fisher's exact or Mann-Whitney U test, as appropriate. OS was defined as time from diagnosis to the last follow-up or death from any cause, while PFS was defined as time from diagnosis to progression to AML (presence of >20% of blasts in BM or PB) or death from CMML. Survival curves were calculated using the Kaplan-Meier method and log-rank test was used for comparisons between groups. Where disease-modifying treatment was administered, duration of follow-up was censored at the time of starting disease-modifying treatment (specifically, allogeneic stem cell transplantation, aggressive chemotherapy, or hypomethylating agents). Multivariate analysis was performed using Cox proportional-hazards regression model, considering Wald Backward as selection method. Two-sided P values<0.05 were considered as statistically significant. The statistical package SPSS, version 23.0 (SPSS Inc., Chicago, IL, USA) was used for all analyses.

Disclosure of potential conflicts of interest

The authors have no potential conflicts of interest to disclose.

Acknowledgements

The authors would like to thank Lauro Sumoy, Gabriel Rech, Sílvia Marcé and Diana Domínguez for their technical assistance.

Additional information

Funding

This work was supported in part by a grant from the Ministerio de Educación Cultura y Deporte (FPU13/03770); MINECO (BFU2015-66559-P); Instituto de Salud Carlos III, Ministerio de Economia y Competividad (PI 11/02519); Redes 2014 (SGR225, SGR35); and the CERCA Programme of the Generalitat de Catalunya. Research in the Buschbeck lab is further supported by MINECO-ISCIII (PIE16/00011); Deutsche José Carreras Leukaemie Stiftung (DJCLS R 14/16); AFM-Téléthon (AFM-18738); and the Marie Skłodowska Curie Training network ‘ChroMe’ (H2020-MSCA-ITN-2015-675610). With economical support from Fundació Internacional Josep Carreras, Obra Social “la Caixa” and Celgene Spain.