DNA Methylation Profiles Capture Clinical and Molecular Heterogeneity of MDS and Can be Harnessed for the Development of Robust Biomarkers Predictive of Response to Azacitidine

While significant progress has been made to understand the genetic landscape of MDS, less is known about the epigenetic makeup of this disease and how this may impact biology and response to azacitidine (AZA). Thus, we performed genomic, epigenomic and transcriptomic analysis on CD34+ cells from a multicenter cohort of 94 intermediate or higher risk MDS patients treated with AZA who had documented responses per IWG 2006. DNA methylation (DNAme) by ERRBS, RNA-seq, mutational profiling, and detailed clinical, cytogenetic and laboratory data were documented.

We first sought to identify aberrant DNAme and expression (GE) profiles in MDS by comparing them to normal CD34+ controls (NBM; n=22 ERRBS & 14 RNA-seq). Despite prior reports of profound promoter hypermethylation in MDS compared to NBM, we uncovered extensive hypomethylation in MDS, with 9,297 hypomethylated and 5,892 hypermethylated differentially methylated regions (DMR). Hypomethylated DMR preferentially targeted enhancers (54% of DMR vs 34% of background; p < 0.001) at key biological pathways such as WNT, mTOR, Hedgehog and Hippo signaling pathways; hypermethylated DMR were enriched in NF-kappa B and NOTCH signaling pathways. In contrast, GE analysis identified only 66 differentially expressed genes between MDS and NBM (36 up and 30 down). Gene set enrichment analysis (GSEA, , FDR < 0.05) showed that genes related to cell cycle (NES=2.5), cancer (NES=2.4) and inflammation pathways (NES=3) were enriched in MDS.

Next we used unsupervised clustering to determine whether DNAme captures disease heterogeneity, correlating the results with genetic and clinical data. Seven clusters were identified that are enriched for cases harboring specific molecular features. Cluster I (n=13) was enriched for cases with mutations in either ASXL1 (8/13, p=0.01) or RUNX1 (7/13, p=0.02) with co-occurrence of splicing factor mutations (SRSF2, U2AF1 or ZRSR2; 9/13 cases). Cluster III included 7/10 cases harboring mutations in genes related to DNAme (DNMT3A, TET2, or IDH1/2, p=0.029) with DNMT3A mutated in 5/7. Cluster IV (n=17) was not significantly enriched for any single mutation, though it primarily consisted of cases carrying splicing factor mutations (11/17 cases). Cluster V (n=7) was enriched for cases with either BCOR (4/7; p=0.007) or STAG2 (4/7; p=0.016) mutations with co-occurrence of SRSF2/U2AF1 mutations (5/7 cases). Out of 9 cases with SF3B1 mutations in the cohort, 4 were in Cluster VI (n=7, p=0.0012) while Cluster VII (n=18) included 9/17 TP53-mutant cases found in the entire cohort (p=0.0006). Notably, Cluster II (n=14) displayed an overall paucity of mutations (mean 1.7 vs. 2.2-3.7 in other clusters; p=0.007), higher hemoglobin levels compared to the rest of the cohort (mean:10.9 vs 9.14 g/dL; p=0.001) and normal cytogenetics. Within cluster II, a subcluster of 6 cases were enriched for DDX41 mutations (p<0.0001). Notably, 11/14 cases in cluster II responded to AZA (≥HI, p=0.0029).

We hypothesized that this epigenetic information at diagnosis may help identify AZA-sensitive patients. A direct comparison between sensitive and resistant cases at diagnosis identified 1,980 hypermethylated and 2,280 hypomethylated DMR in AZA-sensitive cases. In contrast, no significant differences in GE were observed. Using a machine-learning approach (K neighbors), we harnessed these differences to develop a molecular biomarker predictive of response to AZA. After randomly dividing our cohort, we used 61 cases for feature selection and training, and identified 22 predictive DMRs. Testing of the predictor in a 27-patient validation cohort yielded an AUC score of 0.85, indicating excellent discrimination. To further improve this performance, we integrated GE data by incorporating information from the strongest enriched gene set from GSEA. This integrative biomarker, with 21 DNAme + 21 GE features, had improved predictive performance, with AUC score=0.92.

Taken together, our results demonstrate that aberrant DNAme in MDS is not distributed randomly but rather is highly correlated with disease phenotypes, capturing clinically relevant heterogeneity, beyond what is identified by other methodologies. This epigenetic information can be harnessed for the development of robust biomarkers predictive of AZA response. Importantly, integrative approaches combining GE and DNAme data can further improve the predictive performance of these biomarkers.