Specific molecular signatures predict decitabine response in chronic myelomonocytic leukemia

Chronic myelomonocytic leukemia (CMML) is a myelodysplastic syndrome/myeloproliferative neoplasm (MDS/MPN) overlap syndrome (1) that was historically classified within MDS (2) until 2001 (3). CMML shares many characteristics with MDS, including dysplasia in one or more myeloid cell lineages and increased risk of transformation into acute myeloid leukemia (AML). However, a distinguishing feature of CMML is the presence of persistent peripheral monocytosis (>1 × 10⁹/l). CMML can be subdivided into 2 subtypes on the basis of blast count: CMML1, with less than 10% bone marrow (BM) blasts, and CMML2, which has between 10% and 19% blasts.

Substantial epigenetic abnormalities have been described in both MDS and MDS/MPN. Mutations in epigenome-modifying enzymes are highly prevalent in these disorders, including those responsible for DNA methylation and demethylation — DNA methyltransferase 3A (DNMT3A) (4) and ten-eleven translocation 2 (TET2) (5, 6), respectively — as well as those involved in histone-modifying complexes — additional sex combs-like 1 (ASXL1) (7) and enhancer of zeste homolog 2 (EZH2) (8–11). Although the precise mechanisms through which these mutations drive the aberrant epigenetic changes observed in MDS are still not completely understood, it has been shown that MDS and MDS/MPN are characterized by a DNA hypermethylation that increases with disease severity (12, 13).

MDS and MDS/MPN are resistant to conventional chemotherapies; however, epigenome-modifying drugs can be used successfully as therapeutics to treat these disorders. In particular, the nucleoside analogs azacytidine (AZA) and decitabine (DAC) are commonly used to treat MDS and CMML (14, 15). Both AZA and DAC are DNA methyltransferase inhibitors (DMTis), and while their precise mechanism of action in treating MDS and MDS/MPN remains a point of controversy, they are known to be incorporated into DNA during the S phase, where they covalently trap DNA methyltransferases and target them for proteasomal degradation (16, 17). DMTis can also cause DNA damage (18), and because AZA is mostly incorporated into RNA, it may have additional effects on RNA processing and translation (19). Despite the utility of DAC and AZA, only a subset of MDS and CMML patients respond to them. Only approximately 50% of patients treated with DMTis show a hematological improvement (HI) or better that is associated with a survival benefit (20). Furthermore, as many as 6 months of treatment may be required for the therapeutic benefit of DMTis to become apparent, thus forcing half of the patients to undergo long periods of treatment before they can be deemed resistant to this therapy. Currently, there are very few means of predicting response versus resistance, and even this is exclusive to AZA (21). Additionally, few alternative treatments exist for patients who fail to respond to DMTis, and their prognosis is extremely poor. Therefore, it is critical that we better understand the molecular profiles associated with sensitivity and resistance to DMTis in order to improve risk stratification strategies as well as shed light on the mechanisms of resistance.

While some studies have suggested that reversal of methylation and/or transcript reexpression of certain loci was associated with clinical response to DMTis (22–28), epigenetic studies to date have failed to identify any strong correlation between response to these agents and the presence of specific baseline DNA methylation profiles (23, 26, 27, 29, 30). We hypothesized that this lack of correlation was due to the promoter-centric nature of assays used over the past decade and that methylation differences associated with potential for therapeutic response were likely present in these patients upon diagnosis at promoter-distal and intergenic regulatory regions. In this study, we report, for the first time to our knowledge, the identification of DNA methylation and expression differences in diagnostic BM specimens from a cohort of CMML patients treated with DAC. These differences, detected through the use of genome-wide next-generation sequencing assays, reveal underlying biological differences between these 2 groups of patients and point to a novel mechanism of resistance to DMTis.

Somatic mutations do not correlate with response to DAC in CMML. Somatic mutations in epigenome-modifying enzymes and other genes are prevalent in MDS and CMML (4–6, 31–35). Recently, it has been reported that mutations in TET2 and DNMT3A are associated with improved response to DMTi therapy in MDS and related disorders (36–38). Despite this, the presence of these mutations did not translate to an improved overall survival rate in any of these studies, indicating that therapeutic response and survival benefit are likely influenced by multiple different factors. Moreover, these findings have not been recapitulated in CMML exclusively (39). To determine whether particular genetic or epigenetic abnormalities are associated with DMTi sensitivity or resistance in this disease, we studied a cohort of primary CMML cases. BM mononuclear cells (BM MNCs) were collected from 40 patients with de novo CMML at the time of their diagnosis. All patients included in this study were enrolled in a clinical trial conducted by the FISM and received single-agent treatment with DAC as frontline therapy (20 mg/m²/day for 5 days), and response was evaluated after 6 cycles of treatment. Responsive patients (n = 20) were defined as those who achieved either complete remission, marrow complete remission, partial remission, or HI, as defined by the 2006 International Working Group (IWG) response criteria for myelodysplasia (40). Patients with either stable disease or progressive disease were considered to have primary resistance to DAC (n = 20). As shown in Table 1, there were no significant differences in terms of age, gender, BM monocytosis, blast percentage, cytogenetics, or presence of either splenomegaly or extramedullary lesions between responder and nonresponder patients. Using MiSeq to sequence DNA isolated from the diagnostic BM MNCs, we performed targeted resequencing of the following panel of genes mutated at frequencies greater than 5% in CMML: SRSF2, TET2, ASXL1, NRAS, DNMT3A, RUNX1, U2AF1, TP53, JAK2, KIT, KRAS, SF3B1, EZH2, IDH1, and IDH2. As with previous reports, SRSF2, TET2, and ASXL1 were the most frequently mutated genes in this cohort of patients (6, 32, 34, 35, 41–44). However, no somatic mutation was significantly correlated with response to DAC in our cohort (Fisher’s exact test, P = NS for all mutations ) (Figure 1A and Table 2).

Figure 1

Somatic mutations in CMML do not correlate with DAC response or specific epigenetic clusters. Mutational status of a panel of 15 genes frequently mutated in CMML according to (A) therapeutic response to DAC or (B) DNA methylation hierarchical clustering.

Table 1

Clinical characteristics of the FISM CMML patient cohort treated with DAC

Table 2

Somatic mutations of the FISM cohort did not correlate with response

We have previously shown, as have others, that distinct DNA methylation profiles in AML and acute lymphoid leukemia (ALL) are strongly correlated with the presence of specific molecular and cytogenetic subtypes (12, 45–48). To determine whether similarly distinct methylation patterns in CMML can be linked to the presence of specific somatic mutations, we examined DNA methylation patterns in the same specimens through enhanced reduced representation bisulfite sequencing (ERRBS) (45), a deep-sequencing method that captures and accurately quantifies DNA methylation at approximately 3 million CpG sites. ERRBS data were available for 39 of the 40 patients (19 nonresponders and 20 responders). The percentage of methylation measured by ERRBS was highly concordant with the findings of the quantitative single-locus DNA methylation validation assay MassARRAY Epi­TYPER (ref. 49 and Supplemental Figure 1; supplemental material available online with this article; doi:10.1172/JCI78752DS1). Unsupervised clustering analysis of the patients based on their DNA methylation patterns did not reveal a correlation between gene mutations and particular methylation clusters (Figure 1B). In addition, there was no significant difference in the observed patient survival time between the 2 top-level methylation clusters (log-rank test, P = 0.33).

Next, we performed supervised analyses comparing TET2, ASXL1, DNMT3A, and SRSF2 WT and mutant cases to identify the differentially methylated regions (DMRs) associated with each of these mutations. As expected, given its role in de novo DNA methylation, we identified a predominantly hypomethylated profile associated with DNMT3A mutations (total DMRs: 243; hypomethylated DMRs [hypo-DMRs]: 197; hypermethylated DMRs [hyper-DMRs]: 46) that was targeted mainly at intergenic and intronic regions (Figure 2A). By contrast, TET2 loss-of-function mutations were associated with the presence of hypermethylation compared with that seen in TET2 WT cases (total DMRs: 188; hypo-DMRs: 48; hyper-DMRs: 140) (Figure 2B). Mutations in ASXL1, another epigenetic modifier, were associated with a specific signature consisting of equal proportions of hyper- and hypo-DMRs (total DMRs: 144, hypo-DMRs: 82, hyper-DMRs: 62). Both hyper- and hypo-DMRs in ASXL1-mutant CMML cases were strongly depleted at promoter regions (hyper-DMRs 3% vs. background 21%, P = 6.79 × 10^–5; hypo-DMRs 5% vs. background 21%, P = 4.30 × 10^–5) and significantly enriched at intergenic regions (hypo-DMRs 54% vs. background 38%, P = 2.84 × 10^–3) (Figure 2C). Notably, mutations in the splicing factor SRSF2 were linked to the strongest DNA methylation differences, with a total of 724 DMRs (hypo-DMRs: 383; hyper-DMRs: 341). In this case, hypermethylated DMRs were strongly enriched at promoter regions (hyper-DMRs 31% vs. background 21%, P = 1.44 × 10^–5) and depleted at introns (hyper-DMRs 19% vs. background 33%, P = 1.50 × 10^–8) (Figure 2D). While SRSF2 itself does not have any direct epigenetic function, it is likely that mutations in this gene lead to mis-splicing and the consequent deregulation of other epigenome-modifying genes, resulting in this strong epigenetic signature. Additionally, the observed survival time was not significantly different between the patients with or without individual DNMT3A, TET2, ASXL1, and SRSF2 mutations (log-rank test, P = 0.61, 0.067, 0.93, and 0.58, respectively).

Figure 2

Distinct DNA methylation profiles are associated with recurrent somatic mutations in DNMT3A, TET2, ASXL1, and SRSF2. Volcano plots illustrating the methylation differences between DNMT3A-mutant (n = 5) (A), TET2-mutant (n = 17) (B), ASXL1-mutant (n = 15) (C), or SRSF2-mutant (n = 21) (D) samples versus WT patients (n = 39 for the number of mutated samples). DMRs are indicated by red dots (beta-binomial test, FDR <0.1 and absolute methylation different ≥25%). Pie charts illustrate the relative proportion of CpG tiles and DMRs annotated to the RefSeq promoter, exonic, intronic, and intergenic regions.

A specific epigenetic profile distinguishes DAC-resistant CMML patients at diagnosis. Previous efforts by many groups, including ours, have failed to identify baseline epigenetic differences between DMTi-sensitive and -resistant patients (12, 27, 30). However, all of these studies were performed using platforms that examined DNA methylation within CpG islands and gene promoters. A growing body of recent evidence suggests that DNA methylation and other epigenetic modifications at enhancers and other distal regulatory regions play a key role in transcriptional regulation and that these regions are often located at a significant distance from the transcription start site of the target gene (50). Therefore, we hypothesized that key epigenetic differences may exist between DAC-sensitive and -resistant patients at diagnosis that are located distally from promoters, targeting enhancers and other distal regulatory regions.

For this purpose, we used the ERRBS assay, a deep-sequencing–based method that targets not only promoter regions but also intronic, exonic, and distal intergenic regions (45). Using the MethylSig package, we performed a direct comparison between the diagnostic DNA methylation profiles of DAC-sensitive and DAC-resistant patients (51). We identified 167 DMRs that displayed a methylation difference of 25% or more between responders and nonresponders and that were statistically significant at an FDR of less than 0.1. Among these DMRs were regions displaying higher methylation in responders, as well as regions of lower methylation as compared with those in nonresponders (Figure 3A and Supplemental Table 1). Hierarchical clustering of our cohort using these DMRs was sufficient to achieve a perfect segregation of DAC-sensitive and -resistant patients (Figure 3B). These findings indicate that numerous epigenetic differences exist at the time of diagnosis that correlate with a patient’s likelihood of responding to DAC treatment.

Figure 3

Baseline DNA methylation differences distinguish DAC responders and nonresponders at the time of diagnosis. (A) Volcano plot illustrating methylation differences between 20 DAC-sensitive and 19 DAC-resistant patients. Mean methylation difference between the 2 groups is represented on the x axis and statistical significance (–log₁₀ P value) on the y axis. Beta-binomial test identified 167 DMRs, which are indicated by red dots (FDR <0.1 and absolute methylation difference ≥25%). (B) Hiearchical clustering of the patients using the 167 DMRs illustrates the power of these genomic regions in segregating the patients into nonresponder (blue) and responder (red) groups.

Response-associated DMRs localize preferentially to distal regulatory regions. Next, we sought to determine whether DMRs were distributed evenly across the genome or whether they were enriched at specific genomic regions. For this, we analyzed both the genomic distribution of DMRs as well as their association with known regulatory regions. Notably, our analysis of the distribution of DMRs relative to coding regions revealed that DMRs were significantly depleted at promoter regions (DMRs 10% vs. background 21%, binomial test P = 6.70 × 10^–5), with a concurrent enrichment at intronic regions, thus confirming our initial hypothesis. This distribution was not the same across hyper- and hypo-DMRs. While all DMRs were depleted at promoter regions, hyper-DMRs were significantly enriched at introns (hyper-DMRs 49% vs. background 33%, binomial test P = 1.29 × 10^–3), while hypo-DMRs were enriched at intergenic regions (hypo-DMRs 49% vs. 38% background, binomial test P = 0.03) (Figure 4A).

Figure 4

DMRs are enriched at distal intergenic regions and enhancers. (A) Pie charts illustrate the relative proportion of CpG tiles and DMRs annotated to RefSeq promoter, exonic, intronic, and intergenic regions. (B) Pie charts illustrate the relative proportion of CpG tiles and DMRs annotated to CpG islands, CpG shores, and regions beyond CpG shores. (C) Pie charts illustrate the relative proportion of CpG tiles and DMRs annotated to enhancers within gene bodies, enhancers within intergenic regions, and nonenhancer regions.

Next, we sought to determine the association of DMRs with regulatory regions. For this purpose, we analyzed their relative enrichment at CpG island and enhancer regions. Analysis of CpG islands and CpG shores demonstrated that DMRs were also significantly depleted at CpG islands (DMRs 14% vs. background 25%, binomial test P = 2.8 × 10^–4), with enrichment at CpG shores (DMRs 22% vs. background 15%, binomial test P = 8.79 × 10^–3). This pattern was conserved across both hyper- and hypo-DMRs (Figure 4B).

Recently, DNA methylation at enhancers was reported to strongly correlate with aberrant gene expression observed in cancer cells (52). We hypothesized that differential DNA methylation at enhancers, rather than at promoters, may be better correlated with differential response to DAC in CMML. Enrichment analysis of all DMRs relative to intragenic and intergenic enhancers revealed that DMRs were enriched for intragenic enhancers (DMRs 25% vs. background 18%, binomial test P = 0.01). When this analysis was stratified into hyper- and hypo-DMRs, it became apparent that hyper-DMRs showed the strongest enrichment at enhancer regions and, in particular, at enhancers located within gene bodies (hyper-DMRs 32% vs. background 18%, binomial test P = 8.14 × 10^–4). Conversely, hypo-DMRs were not significantly enriched at enhancer regions and were similarly distributed within gene body and intergenic enhancers (Figure 4C).

Finally, we asked whether the DMRs associated with DAC response were specifically enriched within relevant biological pathways. The 167 DMRs were annotated to known genes, and pathway enrichment analysis was performed against the KEGG pathway database. The MAPK signaling pathway, which plays a key role in the cell cycle, apoptosis, cell proliferation, and differentiation, was significantly enriched in DMR-associated genes (hypergeometric test P = 7.68 × 10^–3, FDR = 0.084) (Supplemental Figure 2A). There were 7 DMRs that were annotated to MAPK pathway genes, including STMN1, CACNAE1, PRKCB, MAPT, NFATC1, CRKL, and MKNK2 (Supplemental Table 2). Three of these DMRS — those annotated to STMN1, CACNAE1, and MAPT — were hypermethylated in DAC nonresponders, while MKNK2-, NFATC1-, CRKL-, and PRKCB-associated DMRs were hypermethylated in DAC responders (summarized in Supplemental Table 2). To further validate epigenetic deregulation of the MAPK signaling pathway in these patients, we performed MassARRAY EpiTYPER analysis of 3 of the affected MAPK genes in the pathway in a subset of samples (Supplemental Figure 2B). This analysis confirmed the increased methylation in the STMN1 and CACNAE1 DMRs in nonresponder patients, as well as validated the increased methylation of the NFATC1 DMR in responders.

DNA methylation differences can be harnessed for therapeutic response prediction. Given that our data identified, for the first time to our knowledge, the existence of baseline DNA methylation differences between DAC responders and nonresponders prior to DAC treatment, we hypothesized that these unique methylation profiles could be harnessed to predict at the time of diagnosis which patients would be sensitive or resistant to treatment. To test this, we used the percentage of cytosine methylation at each genomic location among patients in the FISM cohort (cohort 1) as potential predictors and applied a machine-learning approach, support vector machine (SVM) (53), to build a classifier (see details in Methods). Twenty-one 25-bp tile regions were identified by feature selection as the predictors with the highest predictability in the SVM classifier (Figure 5A, Supplemental Figure 3A, and Supplemental Table 3). Unsupervised analysis using only the methylation levels at the 21 selected tile regions revealed that they were sufficient to almost separate the 39 samples by response (Figure 5B and Supplemental Figure 3, B and C). There was no defined clustering of the patients according to their specific degree of response as shown by multidimensional scaling (MDS) analysis (Supplemental Figure 3C), which is concordant with the fact that the classifier was built to identify an all-or-nothing response versus no response and not to distinguish between types of responses. Ten-fold cross-validation was performed using the cases from cohort 1 to evaluate the predictive performance of the classifier, and the reported area under the receiver operating characteristic curve (ROC-AUC) was 0.99, indicating a strong predictive accuracy for the classifier model (Supplemental Figure 3D). In order to further assess the robustness of the SVM classifier built with the 21 selected features, we performed 3 different random splits of the same cohort 1 into training and test sets. We trained the classifier on each of the 3 sets of randomly selected samples and predicted the responses for the remaining samples in the cohort. The classifier was able to accurately predict response to DAC in 18 of 19 (accuracy = 94.74%) (Table 3), 13 of 14 (accuracy = 92.86%), and 9 of 9 (accuracy = 100%) patients, respectively (Supplemental Figure 4A).

Figure 5

Methylation profiles can be harnessed to classify patients according to DAC response at diagnosis. (A) Heatmap of 21 CpG tiles selected as the SVM classifier predictors. DAC-sensitive patients are indicated with the red bar and nonresponders with the blue bar. (B) Correspondence analysis (COA) using only the 21 CpG tiles included in the classifier could segregate the majority of the CMML cohort according to DAC response (responders are represented by red dots and nonresponders by blue dots).

Table 3

Prediction performance of the SVM classifier trained on 20 randomly selected samples and applied to the remaining 19 samples in the FISM cohort (accuracy = 94.74%)

Since validation in an independent cohort of patients is the gold standard for biomarker development, we identified a second cohort of patients in which to test the performance of our SVM classifier. Twenty-eight additional diagnostic CMML specimens from patients enrolled in a clinical trial from the Groupe Francophone des Myelodysplasies (GFM), all of whom had been treated with the same DAC regimen of 20 mg/m²/day for 5 days, were collected and subjected to ERRBS (Table 4 and Supplemental Table 4). Specimens from this second cohort (cohort 2) of 12 responder and 16 nonresponder patients consisted of sorted monocytes from peripheral blood (PB). The SVM classifier that had been developed using cohort 1 was applied blindly to these samples, without any prior knowledge of the therapeutic response labels for this second cohort. Due to the stochastic nature of ERRBS, the CpG coverage is never identical across all samples, thus leading to missing values for some regions of interest. In effect, only 6 of the 21 features were present in all 28 samples in cohort 2. Therefore, using only these 6 features, we first trained our SVM classifier on the 39 samples of the FISM cohort (cohort 1) and then applied the trained classifier on the GFM cohort (cohort 2). As shown in Table 5 and Supplemental Figure 4B, despite this limitation, the 6-feature classifier was still capable of correctly predicting response in 20 of 28 patients in the GFM cohort (accuracy = 71% and AUC = 0.82). Next, in order to increase the number of features being tested while still retaining a large enough cohort in which to test the predictive accuracy, we used 14 of the 21 features of the SVM classifier to predict response for 19 patients in the GFM cohort. Once again, we used only these 14 features to train the model on cohort 1, which consisted of the initial 39 patients, and then blindly applied the model to the 19 test samples from the GFM cohort (cohort 2). This modified classifier with 14 features was capable of accurately predicting therapeutic outcome for 15 of the 19 patients, which represents an accuracy of 79% and an AUC of 0.83. (Table 5 and Supplemental Figure 4B). Finally, we determined that of the original 21 features, 16 was the maximum number of features shared by at least 15 of the cohort-2 patients. We trained the model on cohort 1 using only these 16 shared features and then predicted response for the 15 patients in the independent cohort 2, achieving an accuracy of 87% with an AUC of 0.94 (Table 5 and Supplemental Figure 4B). These findings demonstrate that the SVM classifier developed using the original FISM cohort is general enough to be applied to and accurately predict the therapeutic outcome of fully independent samples (i.e., GFM cohort 2), which is a critical step in the development of a biomarker. Moreover, this robustness was maintained even across different cell types (BM MNCs in cohort 1 vs. PB monocytes in the validation cohort 2), further underscoring the power of the classifier to predict outcome in an independent cohort. While further validation in larger cohorts will be required to fully assess the accuracy of the features reported here, and additional studies of larger cohorts might help refine the selection of features to include those with the strongest accuracy over a large number of patients, our findings demonstrate that the epigenetic differences between responders and nonresponders at diagnosis have the potential to be harnessed as classifiers to predict clinical response to DAC.

Table 4

Clinical characteristics of the GFM CMML cohort treated with DAC

Table 5

Summary of the prediction performance of the independent validation cohort (GFM) in 3 scenarios using an increasing number of shared features of the 21 features preselected from the FISM cohort

DAC sensitivity can be linked to a specific transcriptional program at diagnosis. While it has been previously shown that reduced expression of uridine-cytidine kinase, an enzyme involved in nucleoside metabolism, is associated with response to AZA in MDS (54), we did not find that differential expression of this or other DMTi-metabolizing enzymes was associated with response to DAC in CMML (data not shown). Therefore, we sought to determine whether other transcriptional differences between DAC responders and nonresponders are indicative of response and can provide insight on functional pathways that contribute to DAC resistance. We performed RNA-sequencing (RNA-seq) on samples from 14 patients (8 responders and 6 nonresponders) in the cohort of CMML patients for whom we had high-quality RNA. Prior to performing differential analysis, we validated the ability of our RNA-seq approach to accurately detect quantitative variability by performing quantitative reverse transcriptase PCR (qRT-PCR) on RNAs from 13 of the 14 patients and determining the degree of agreement between the 2 methods (r = 0.85, R² value = 0.73, P < 0.0001) (Supplemental Figure 5A). As shown in Figure 6A, a direct comparison of the 2 groups of patients identified 601 genes with an absolute log₂ fold change greater than 1 and a P value of less than 0.05. Notably, this gene signature consisted of a majority of genes overexpressed in DAC-sensitive patients (405 upregulated genes), with only a small proportion of genes downregulated in these patients (Supplemental Table 5).

Figure 6

A specific transcriptional program is associated with response to DAC. (A) Heatmap illustrates gene expression differences between 8 DAC-sensitive (indicated by the red bar) and 6 DAC-resistant patients (indicated by the blue bar). Genes represented in the heatmap were identified by a GLM likelihood ratio test (P < 0.05 and absolute log₂ fold change >1). (B) Enrichment plots for GSEA using the expression difference–ranked gene list showing enrichment for cell cycle–related gene sets. NES, normalized enrichment score. (C) Box plots showing gene expression differences for CXCL4, CXCL7, and ITGB3 (red box plots denote responders; blue box plots denote nonresponders). P values were obtained from a GLM likelihood ratio test.

In order to identify biological differences that might explain the difference between these patients in their therapeutic response to DAC, we performed gene set enrichment analysis (GSEA) (55). Gene sets enriched in DAC-sensitive patients at an FDR of less than 0.1 were involved in proliferation, cell cycle activity, and DNA replication (Figure 6B). Likewise, genes reported as being downregulated in quiescent versus dividing CD34⁺ cells (56) were found to be upregulated in DAC responders. This enrichment of gene sets involved in the cell cycle and in DNA replication in DAC-sensitive patients is consistent with the requirement for DAC incorporation into the DNA during the S phase.

Primary resistance to DAC is associated with overexpression of ITGβ3 and the chemokines CXCL4 and CXCL7. As mentioned above, only a small fraction of genes were found to have at least a 2-fold overexpression in DAC-resistant patients. Among these, 3 genes that have previously been implicated in chemoresistance and leukemogenesis were overexpressed in nonresponders: CXCL4 (also known as PF4), CXCL7 (also known as PPBP), and integrin β3 (ITGB3) (Figure 6C). Thus, we hypothesized that overexpression of these genes might be a potential mechanism through which CMML acquires resistance to DAC. First, as shown in Figure 7A, we validated the overexpression of these genes in DAC-resistant patients by qRT-PCR. Notably, there was a statistically significant linear correlation between the levels of CXCL4 and CXCL7 expression by both RNA-seq (r = 0.9350, R² = 0.87, P < 0.0001) and qRT-PCR (r = 0.9865, R² = 0.9731, P < 0.0001), suggesting that these factors act in concert in the BM microenvironment (Figure 7B). While both chemokines were originally thought to be produced exclusively by megakaryocytes, there is evidence that monocytes (57, 58) and other cells within the BM also produce CXCL4 and CXCL7 (refs. 59, 60, and Supplemental Figure 5, B and C). To further confirm the overexpression of these chemokines in nonresponder patients as well as to determine the cellular source and localization of the proteins in the BM, IHC was performed on a subset of paraffin-embedded BM biopsies taken at diagnosis from responders and nonresponders. As shown in Figure 7, C and D, CXCL4 was primarily localized to megakaryocytes, while CXCL7 staining was stronger in an MNC population compatible with a monocytic origin. Importantly, there was increased CXCL4 and CXCL7 staining in BM from nonresponder patients as compared with that in BM from responders, confirming the presence of CXCL4 and CXCL7 proteins in the BM microenvironment, which, like mRNA levels, are increased in DAC-resistant patients.

Figure 7

CXCL4 and CXCL7 are upregulated in the BM of nonresponders. (A) qRT-PCR showing validation of overexpression of CXCL4, CXCL7, and ITGB3 in nonresponders; each point represents the mean of triplicate wells for each patient sample; the line and error bars indicate the group mean and SD, respectively. (B) Pearson’s correlation analysis of expression levels of CXCL7 and CXCL4 by RNA-seq and qRT-PCR. (C and D) Representative IHC images for CXCL4 (C) and CXCL7 (D) in diagnostic BM biopsies in DAC responders and nonresponders. Original magnification, ×40 (C and D, left panels), ×63 (C and D, right panels). Representative images from duplicate experiments are shown.

Previous studies have implicated serum levels of CXCL4 and CXCL7 as potential prognostic markers in MDS (61, 62). To determine whether serum levels of CXCL4 and CXCL7 could potentially serve as biomarkers for DAC response, we quantified the serum concentrations of these chemokines by ELISAs in 35 of 40 CMML patients (Supplemental Figure 6). There was no significant difference in serum CXCL4 or CXCL7 levels between responders and nonresponders. In addition, we found no significant correlation between BM mRNA levels and serum protein levels for these 2 chemokines, indicating that serum levels of these chemokines are not reflective of mRNA expression in the BM and mirroring previous observations documented for other chemokines in the BM and serum of AML patients (63, 64).

CXCL4 and CXCL7 abrogate the effect of DAC on hematopoietic cells. It has been previously reported that both CXCL4 and CXCL7 can reduce the chemosensitivity of BM cells to 5-fluorouracil in vitro (65), and CXCL4 has been implicated in cell cycle arrest (66) and quiescence (67, 68), which might be a mechanism through which it acts to prevent sufficient incorporation of DAC into cells of nonresponders. Therefore, we hypothesized that an overabundance of CXCL4 and CXCL7 in the BM microenvironment acts to overcome the effects of DAC. To test this, we cultured primary human CD34⁺ cells for 3 days in vitro with CXCL4 (50 ng/ml), CXCL7 (50 ng/ml), or a combination of both chemokines in either the presence or absence of low-dose DAC (10 nM) and then plated them in methyl­cellulose to test their clonogenic potential. The chemokines and low-dose DAC did not affect cell proliferation during the in vitro liquid culture period (Supplemental Figure 7A). Moreover, as previously reported, low-dose DAC did not reduce cell viability or induce apoptosis after 3 days in culture (Supplemental Figure 7, B and C, and ref. 69). However, 3 days of treatment with 10 nM DAC significantly reduced colony formation. Addition of either CXCL4 or CXCL7 alone did not have a significant impact on DAC-induced colony inhibition. However, concomitant treatment of CD34⁺ cells with CXCL4 and CXCL7 completely abolished the suppressive effect of DAC on colony formation (Figure 8A).

Figure 8

CXCL4 and CXCL7 promote resistance to DAC in CD34⁺ and primary CMML specimens. (A) Colony formation was inhibited by DAC but restored with the combination of CXCL4 and CXCL7. CD34⁺ cells were treated with 1 dose of CXCL4, CXCL7, or both (50 ng/ml each) or with vehicle (PBS containing 0.1% BSA) and daily 10-nM doses of DAC for 3 days. After 3 days of in vitro treatment with DAC, cells were plated in methylcellulose and incubated for 12 to 15 days before colonies were counted. Data represent the mean ± SD. Treatment with 10 nM DAC significantly decreased colony formation but failed to do so in the presence of CXCL7 and CXCL4 together. Shown in the 3 panels are the results of 3 independent experiments. Error bars represent the SD. (B) CXCL4 and CXCL7 abrogated the effect of DAC on the viability of primary CMML MNCs. CMML MNCs were treated in vitro for 72 hours with 10 nM DAC alone or in the presence of 50 ng/ml CXCL4, CXCL7, or both. Data represent the mean ± SD. Treatment with DAC alone significantly reduced the viability of these cells, but this effect was lost when CXCL4 or CXCL7 was added to the culture. All data represent independent experiments performed in 3 different CMML patients. Error bars represent the SD. *P < 0.05 and **P < 0.01 by unpaired 2-tailed Student’s t test.

Finally, we tested the ability of CXCL4 and CXCL7 to induce resistance in primary CMML cells. BM MNCs from diagnostic specimens collected from 3 patients were placed in liquid culture and treated for 72 hours with 10 nM DAC in the presence or absence of 50 ng/ml CXCL4, CXCL7, or a combination of both. Viability was assessed after 72 hours. Unlike normal CD34⁺ cells, which did not show diminished viability with 10 nM DAC (Supplemental Figure 5B), treatment of primary CMML cells with low-dose DAC led to a significant decrease in viability in all 3 patients (P < 0.01). However, concomitant treatment of CMML cells with CXCL4, CXCL7, or their combination abrogated the effect of DAC on all 3 patients (Figure 8B). Combined, these data support the hypothesis that the presence of excess CXCL4 and CXCL7 in the marrow microenvironment contributes to induction of DAC resistance in CMML cells.

While DMTis remain the only FDA-approved therapy for the majority of MDS and nonproliferative CMML patients, prognosis following DMTi treatment failure is extremely poor, with median survival for these patients barely reaching 6 months and approximately 50% of patients never even achieving a response in the first place (20, 70). This relatively low rate of therapeutic response is further complicated by the slow kinetics of DMTis, which may take as long as 6 to 12 months to show efficacy, thus committing the majority of patients to receive a drug to which they will ultimately be deemed resistant. Therefore, we set out to study the epigenetic and transcriptional characteristics associated with response to DAC in a cohort of CMML patients in order to identify molecular features that allow risk stratification at the time of diagnosis and, additionally, to explain the mechanisms behind the primary resistance to this agent. To better understand the molecular and mechanistic basis for DMTi response and effectively risk-stratify patients at diagnosis, we performed next-generation sequencing assays to study both the epigenome and the transcriptome of a uniformly treated cohort of CMML patients who differed in their response to DAC. The use of this improved technology, with extended genomic coverage and better dynamic range, allowed us to detect, for the first time to our knowledge, the presence of DNA methylation and gene expression differences present at the time of diagnosis that distinguish DMTi-sensitive and -resistant patients. The enrichment of these DMRs at distal enhancers, as well as the depletion of promoter-associated DMRs identified in this baseline epigenetic signature, underscores the importance of analyzing DNA methylation changes beyond promoter regions and explains the lack of statistically significant differential methylation observed in previous studies that were confined solely to promoter methylation analysis (12, 27, 30).

Moreover, our observation that the genomic locations predominantly affected by differential DNA methylation are distal regulatory regions adds more data to the strong evidence that emphasizes the critical role of long-range epigenetic gene regulation. Techniques to examine 3D chromatin architecture, such as chromosome conformation capture (3C) (71) and its subsequent iterations 4C (72, 73), 5C (74), and Hi-C (75), have indicated that gene regulation often occurs at very distant locations, in part through DNA looping at distal enhancers. In fact, only a small percentage (~7%) of gene-looping events have been reported to involve the nearest gene transcription start site (50). This argues for the critical role of distal, nonpromoter regulatory regions in controlling gene expression. If the differential methylation at nonpromoter regions does impact the expression of long-range target genes, this may explain why several previous studies have struggled to correlate differential DNA methylation with gene expression changes using nearest-gene annotations (30, 76).

We found that the MAPK pathway was significantly enriched in DMRs, with both gains and losses of methylation in responders and nonresponders within this pathway. These DMRs were localized to both intra- and intergenic genomic regions annotated for 7 genes involved in the MAPK pathway. While in-depth functional analysis of these DMRs will be required in additional experiments that are beyond the scope of our study, our findings support the results by others suggesting the importance of aberrant MAPK pathway signaling in contributing to MDS/MPN (77, 78), as well as to drug resistance and cell cycle progression in leukemic cells (79, 80). Furthermore, while it is known that multiple genes in the MAPK pathway can be mutated in CMML (81), our results indicate that the epigenetic alterations of genes in this pathway may also be present in CMML patients.

While previous reports on MDS and related malignancies have linked the presence of certain mutations — specifically, TET2 (36–38) and DNTM3A (37) — to an increased rate of response to DMTis, we could not find any correlation between the mutational status of these and other genes commonly mutated in CMML and response to DAC in our FISM cohort. This finding is in concordance with those of a previous report on CMML (39), which likewise failed to detect a correlation between response to DAC and mutational status, indicating that the impact of mutational status may be different in CMML patients compared with that in MDS patients or in mixed cohorts consisting of MDS patients as well as patients with other myeloid malignancies, including AML (37, 38) and MDS/MPN (37). Furthermore, the studies demonstrating better TET2- and DNMT3A-associated responses involved patients treated with AZA alone (38) or cohorts including both AZA- and DAC-treated patients (36, 37), which may also contribute to the differing result obtained in our study on patients who received DAC exclusively.

Conversely, DNA methylation status was indeed different at diagnosis between DAC-sensitive and DAC-resistant patients, and we demonstrate that these differences can risk-stratify patients at the time of diagnosis using an epigenetic classifier that exploits these identified methylation differences. Moreover, the SVM classifier developed in this study performed with 87% accuracy on an independent cohort, even when only a subset of the original features were included and 2 different cell types were used in the training and validation cohorts (BMN MNCs vs. PB monocytes). Thus, while the classifier reported here will require further extensive validation in larger, independent cohorts, the present study demonstrates not only that DNA methylation differences exist between patients with different responses to DAC but that these DNA methylation differences are sufficiently robust to be harnessed for use in the clinic as accurate classifiers. These classifiers have the potential to prevent patients who are unlikely to respond to DAC from receiving prolonged, unwarranted treatments with this drug and instead permit them to be quickly transitioned to alternative therapies.

In addition to epigenetic differences, our study also revealed baseline differences at the transcriptional level that correlated with response to DAC. Analysis of this response-associated signature demonstrated a strong enrichment for gene sets involved in cell cycle regulation among the genes upregulated in DMTi-sensitive patients. This finding is in line with the need for DAC to be incorporated into the DNA during cell cycle activity in order to exert its effects. By contrast, fewer genes were upregulated in resistant patients. Among these overexpressed genes, we found CXCL4 and CXCL7, two chemokines that have been previously implicated in mediating cell cycle arrest (66), quiescence (67, 68), and reduced chemosensitivity of BM cells to 5-fluorouracil in vitro (65). We therefore focused our efforts on studying the impact of these chemokines on response to DAC. In vitro treatment of both normal CD34⁺ cells or primary CMML MNCs with CXCL4 and CXCL7 blocked the effect of DAC on these cells, indicating that overexpression of these 2 genes may indeed lead to primary resistance to DAC and opening the possibility for future targeting of the downstream signaling cascades in order to overcome the effect of these chemokines.

Sample collection and processing

FISM cohort. BM specimens were collected before treatment from 40 patients with CMML. BM MNCs were isolated through Ficoll density centrifugation and viably frozen in 10% DMSO and 90% FBS. Patients with advanced CMML were enrolled in the nonrandomized clinical trial conducted by the FISM (NCT01251627; https://clinicaltrials.gov/) and were given DAC (20 mg/m²/day i.v.) for 5 days every 28 days for at least 6 cycles prior to being classified as responders or nonresponders, with response defined as HI or better according to IWG 2006 criteria (40). The clinical characteristics of the patients are summarized in Table 1. Genomic DNA and total RNA were isolated using the AllPrep DNA/RNA kit (QIAGEN) according to the manufacturer’s instructions.

GFM cohort. The patients were enrolled in the EudraCT 2008-000470-21 GFM trial (NCT01098084; https://www.clinicaltrials.gov/) and received DAC (20 mg/m²/day i.v.) for 5 days every 28 days for at least 3 cycles. Blood samples were collected using EDTA-containing tubes, mononucleated cells were isolated on Ficoll-Hypaque, and monocytes were enriched using the AutoMacs system (Miltenyi Biotec) through negative selection with microbeads conjugated with antibodies against CD3, CD7, CD16, CD19, CD56, CD123, and glycophorin A, then further enriched by positive selection with microbeads conjugated with a monoclonal mouse anti-human CD14 antibody (Miltenyi Biotec). Genomic DNA was extracted from the monocytes using the Norgen Biotek kit (Thorold) kit according to the manufacturer’s instructions. The clinical characteristics of the patients are summarized in Table 4.

Mutational sequencing

Target capture. Capture of the target regions (exons plus splice junctions) was carried out using a custom-designed HaloPlex Target Enrichment kit (Agilent Technologies) following the HaloPlex Target Enrichment System-Fast Protocol, version D.5.

Sequencing. DNA (500 ng) from each sample was quantified with a Qubit Fluorometer (Invitrogen) and used in the capture reaction. Each sample had a unique index. Libraries were quantified by Qubit, pooled, and run in an Illumina HiSeq 2500 rapid-run flow cell using the on-board cluster method for paired-end sequencing (2 × 100 bp reads).

Analysis. Sequencing results were demultiplexed and converted to a FASTQ format using Illumina BCL2FASTQ software. The reads were adapter and quality trimmed with Trimmomatic (82) and then aligned to the human genome (UCSC build hg19) using the Burrows-Wheeler Aligner (83). Further local indel realignment and base-quality score recalibration were performed using the Genome Analysis Toolkit (GATK) (84). Single-nucleotide variation and indel calls were generated with the GATK HaplotypeCaller. ANNOVAR (85) was used to annotate variants with functional consequence on genes as well as to identify the presence of these variants in dbSNP 137, the 1000 Genomes Project, ESP6500 (National Heart, Lung, and Blood Institute [NHLBI] GO Exome Sequencing Project), and COSMIC 67.

Genome-wide DNA methylation by ERRBS

High-molecular-weight genomic DNA (25 ng) was used to perform the ERRBS assay as previously described (45) and was sequenced on an Illumina HiSeq 2000. Reads were aligned against a bisulfite-converted human genome (hg18) using Bowtie and Bismark (86). Downstream analysis was performed using R statistical software (version 3.0.3) (87), Bioconductor 2.13 (88), and the MethylSig 0.1.3 package (51). Only genomic regions with coverage ranging from 10 to 500 times were used for the downstream analysis. DMRs were identified by first summarizing the methylation status of the genomic regions into 25-bp tiles and then identifying regions with an absolute methylation difference of 25% or more and an FDR of less than 10%. DMRs were annotated to the RefSeq genes (NCBI) using the following criteria: (a) DMRs overlapping with a gene were annotated to that gene; (b) intergenic DMRs were annotated to all neighboring genes within a 50-kb window; and (c) if no gene was detected within a 50-kb window, then the DMR was annotated to the nearest transcription start site (TSS).

Methylation classifier

SVM (53) was applied using R package e1071 (89) to classify the 2 groups of patients (responders and nonresponders), in which the percentage of methylation of the 25-bp tiles was used as a predictor. The probability mode and sigmoid kernel were used in the SVM function, otherwise the default parameters were applied. We performed 2-step feature selections for the SVM classifier: (a) 25-bp tiles were prefiltered by nominal P values of less than 0.05 and by an absolute methylation difference greater than 20%, calculated using the MethylSig package (51); (b) greedy forward-feature selection was applied on the remaining tiles. Briefly, we assessed and prioritized the predictability of each of the filtered tiles in the SVM model and then sequentially evaluated the combinatorial predictability of the tiles by adding 1 tile from the prioritized tiles to the classifier at a time. The final predictors of the SVM classifier were selected from the set of tiles that could optimally predict patient response. The predictability was assessed on the basis of 10-fold cross-validation. Specifically, we randomly partitioned the 39 samples for which ERRBS libraries were available into 10 complementary subsets, training the SVM model on 9 of the 10 subsets (called the training set) and predicting the classes (responder or nonresponder) on the 1 left-out subset (called the validation set or testing set). To reduce variability, 10 rounds of cross-validation were performed using different partitions, and the validation results were summarized over the rounds. During each round of validation, the probability of each sample being predicted as a responder was recorded, and then the ROC-AUC across 10 rounds was calculated with the R package ROCR (90), and this calculation was used as the assessment of the predictability. Complete code is provided in the Supplemental Methods.

EpiTYPER MassARRAY

Validation of CpG methylation of select genomic regions was performed by MALDI-TOF using EpiTYPER MassARRAY (Sequenom) (49) on bisulfite-converted genomic DNA from a subset of DAC responders and nonresponders. The primers used to amplify these genomic regions and the resultant amplicon sequences are listed in Supplemental Table 6.

RNA-seq

RNA-seq was performed on RNA samples from 14 patients (8 responders and 6 nonresponders) who had high-quality RNA (RNA integrity number >6 as determined by the Agilent 2100 Bioanalyzer). RNA-seq libraries were prepared using the Illumina TruSeq RNA Sample Prep Kit (version 2) according to the manufacturer’s instructions. A set of synthetic RNAs from the ERCC (91) at known concentrations were mixed with each of the cDNA libraries. Four separate samples were multiplexed into each lane and sequenced on a HiSeq 2000. The quality of reads obtained was evaluated using FastQC (Babraham Bioinformatics; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The sequenced libraries were aligned to the human genome (hg18) or to the ERCC spike-in reference sequence using TopHat, version 2.0.8 (92), with default parameters.

RNA-seq analysis

HTSeq (0.5.4p5) (93) was used to generate the count matrix with the following parameters: “htseq-count --mode=union --stranded=no” using the following 2 gene transfer format (GTF) annotation files, respectively: (a) the hg18 RefSeq gene GTF file downloaded from the UCSC genome browser for endogenous gene assembly; (b) the ERCC spike-in transcript GTF file downloaded from the official website (http://www.lifetechnologies.com/order/catalog/product/4456740) for ERCC spike-in assembly. The endogenous gene counts were normalized by ERCC spike-in library size, and the differential expression analysis was performed using the edgeR (version 3.4.2; Bioconductor) (94) generalized linear model (GLM). Genes with an absolute log₂ fold change greater than 1 and a P value of less than 0.05 were reported.

qRT-PCR

To validate the RNA-seq results, RNA from selected nonresponder and responder patients was reverse transcribed using the Verso cDNA synthesis kit (Thermo Scientific) with random hexamer primers, according to the manufacturer’s instructions. qRT-PCR was performed on the resulting cDNA in triplicate using intron-spanning and -flanking primer sets with Fast SYBR Green Master Mix and the StepOne Plus PCR System (Applied Biosystems) according to the manufacturer’s instructions. The primer sequences are listed in Supplemental Table 7.

ELISAs

ELISAs for CXCL4 and CXCL7/NAP2 on serum from the CMML patients were performed using the corresponding ELISA kits (RAB0402 and RAB0135) from Sigma-Aldrich according to the manufacturer’s directions. For CXCL4, the serum was diluted 1:500 in the sample dilution buffer provided in the kit.

IHC

For immunostaining, 3-μm-thick formalin-fixed, paraffin-embedded BM sections were deparaffinized in xylenes and hydrated in graded alcohols. Antigen retrieval was performed in EDTA (1 mM, pH 8.0) for two 15-minute cycles at maximum power in a microwave oven, and slides were then incubated with a CXCL4 antibody (1:300, catalog 500-P05; PeptroTech) or a CXCL7 antibody (1:50, catalog orb13423; Biorbyt). Immunostaining was performed with the BenchMark histo­stainer (Ventana Medical Systems, Roche) using a peroxidase detection kit with DAB substrate according to standard procedures. Sections were then counterstained with hematoxylin.

Cell culture and colony-forming assays

CD34⁺ cells were isolated from cryopreserved BM MNCs from femoral head specimens using the CD34 MicroBead Isolation Kit (Miltenyi Biotec) according to the manufacturer’s instructions. For CMML cells, the cryopreserved BM MNCs were rapidly thawed at 37°C and treated with DNAse to prevent cell clumping. Cells were plated in prestimulation media containing IMDM with 20% BIT (STEMCELL Technologies); IL-6 (20 ng/ml); SCF (100 ng/ml); TPO (100 ng/ml); and FLT3L (10 ng/ml) (PeproTech) and recovered overnight. The following day, either CXCL4 (50 ng/ml; PeproTech); CXCL7 (50 ng/ml; PeproTech); a combination of both chemokines (50 ng/ml each); or vehicle (PBS containing 0.1% BSA) was added as well as freshly prepared DAC (10 nM) (Sigma-Aldrich) or vehicle (water). DAC was replenished daily for a total of 3 days. Live cell numbers and viability were determined by trypan blue exclusion.

For colony assays, an equal number of live, treated CD34⁺ cells were plated in duplicate in H4435 Enriched MethoCult (STEMCELL Technologies). Colonies were counted after 12 to 15 days.

Apoptosis assays

Apoptosis was assessed using the Tali Apoptosis Kit with annexin V Alexa Fluor 288 and propidium iodide according to the manufacturer’s instructions and was measured on a Tali Image-Based Cytometer (all from Life Technologies).

Accession numbers

FISM cohort ERRBS and RNA-seq data are deposited in the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE61163). GFM cohort ERRBS data are also deposited in the GEO database (GEO GSE63787).

Statistics

For the analysis of clinical parameters, Fisher’s exact test was used for CMML type and sex; unpaired, 2-tailed Student’s t tests were used for clinical parameters with a normal distribution; Wilcoxon signed-rank tests were used when the samples were not normally distributed; and the log-rank test was used for survival. A P value of less than 0.05 was considered significant. Somatic mutations between nonresponders and responders was evaluated using Fisher’s exact test, and significance was considered at a P value of less than 0.05. For in vitro cell culture and colony-forming experiments, unpaired, 2-tailed Student’s t tests were used for comparisons, and significance was considered at a P value of less than 0.05. For correlation analysis between the RNA-seq and qPCR results, Pearson’s correlation was performed, and the r values and P values are indicated in the figures. The ERRBS and RNA-seq analyses were performed using a beta binomial test for differential methylation and a generalized linear model likelihood ratio test for differential gene expression. These methods were implemented through specific algorithms that are described in detail in their respective sections above.

Study approval

The current study was approved by the IRB of the University of Michigan Medical School and the ethics committee of the University of Florence, AOU Careggi-Firenze. The GFM clinical trial EudtraCT (2008-000470-21; https://eudract.ema.europa.eu/) was approved by the ethics committee of the Centre Hospitalier Universitaire de Dijon (Dijon, France). All samples were obtained from patients enrolled in clinical trials, and written, informed consent was obtained from these patients at the time of their enrollment in the study. Samples used in the current study were deidentified prior to use at the University of Michigan.