Single-cell DNA methylation analysis uncovers epigenetic pathways in the transformation of MDS to AML

Acute myeloid leukemia (AML) can arise from an antecedent myelodysplastic neoplasm (MDS) [ 1 ] despite the treatment with hypomethylating agents (HMA). This secondary AML (sAML) is associated with resistance to therapies and shortened survival [ 1 ]. Both MDS and AML undergo aberrations in DNA methylation [ 2 , 3 , 4 ], but the MDS to AML transformation has been mostly studied at the genetic level [ 5 ]. We wondered if uncovering the single-cell DNA methylation (scDNA-Met) ecosystem could yield relevant clues to understand the biology of these aggressive sAMLs in a similar manner that has been shown for hematopoietic stem cell differentiation during aging [ 6 ]. We dissected in a longitudinal manner the epigenetic phylogeny of patients with MDS at diagnosis (Dx), at the end of HMA therapy (Post-HMA), and upon the rise of sAML using Smart-RRBS (reduced representation bisulfite sequencing) [ 7 , 8 ], which allows single-cell DNA methylation resolution, with complementary studies of RNA, DNA mutational status, and protein expression also at the single-cell level. Since the origin of MDS and AML is mostly confined within the hematopoietic stem and progenitor cells (HSPCs) compartment [ 9 , 10 ], we have herein selected by fluorescence-activated cell sorting HSPCs defined as Lin- CD34+ CD38− expressing at least one of the leukemic stem cell markers (CD45RA or CD123; that is, Lin- CD34+ CD38− CD45RA+ CD123− or Lin- CD34+ CD38− CD45RA− CD123+, Supplementary Methods ), previously used to enrich for MDS and AML stem cells and from now on termed malignant stem cells (mSCs) [ 5 ], from 15 bone marrow (BM) longitudinal samples from 5 patients that transformed from MDS to AML after the treatment with HMA (azacitidine) in monotherapy to delineate the epigenetic single-cell architecture by Smart-RBBS [ 7 , 8 ] (Fig. 1A and Supplementary Methods ). As controls, HSPCs from one age-matched healthy BM (66 years), which do not express either CD45RA or CD123 (Lin- CD34+ CD38− CD45RA− CD123−), were also isolated ( Supplementary Methods ) [ 5 ]. The MDS clinicopathologic characteristics are shown in Supplementary Tables S1 and S2 . HMA treatment induced a complete clinical response in three patients (#3, #26, and #28), whereas two patients (#6 and #8) were classified as non-responders. We profiled single-cell DNA methylation of the three matched samples from the five patients studied and the healthy BM to measure intra-sample epimutation rate (PDR, proportion of discordant reads) variability, reflecting the degree to which individual cells experience epigenetic alterations (epimutations) throughout disease progression. This epimutation rate was obtained by assessing the methylation concordance/discordance of adjacent CpGs within the same sequencing read and then calculating the proportion of discordant reads per cell (Fig. 1A, B and Supplementary Methods ) [ 7 ]. We also measured absolute cell-to-cell PDR differences (ΔPDR) that compute PDR differences between cells from the same sample ( Supplementary Methods ), capturing the variation in the epimutation rate among single cells within a sample. This implies that an increase in ΔPDR values indicates greater epigenetic heterogeneity (Fig. 1A, B and Supplementary Methods ) [ 7 ]. We observed that there was an overall progression towards a higher epimutation rate (PDR) from healthy BM to the final transformation to AML (Fig. 1C, D and case-by-case analyses in Supplementary Fig. S1 ), suggesting an increase in DNA methylation pattern diversity accompanied by higher intercellular heterogeneity in the road toward the most aggressive myeloid malignancy. In the same line, we noticed that sAML samples showed higher PDR levels than the corresponding MDS ones (Fig. 1C, D ). This increase was supported by all HMA-responders, while non-responders showed variable patterns (Supplementary Fig. S1 ). Interestingly, when we subdivided cells according to HMA response, responder patients did not exhibit a higher PDR at MDS diagnosis compared to normal BM (Fig. 1D , top left and Supplementary Fig. S1A ), whereas non-responders showed a higher epimutation rate (Fig. 1D , top right, and Supplementary Fig. S1A ). Also, HMA-responders did not show any significant differences in ΔPDR between MDS and sAML samples, whereas a higher value was found for non-responders (Fig. 1D , bottom). These results suggest a higher intracellular epigenetic variability between responders and non-responders. Fig. 1: Longitudinal changes in epimutation rates (PDR) and epigenetic heterogeneity (ΔPDR) during MDS-to-AML progression. Full size image A Schematic overview of the study design. Image created with BioRender.com. B Epimutation rate calculation as the proportion of discordant reads (PDR, top), adapted from Gaiti et al. [ 7 ], and brief explanation about the relationship between ΔPDR and clonal heterogeneity (bottom). C Top, boxplot showing the epimutation rate (PDR) per cell across the four conditions: healthy bone marrow (hBM), at MDS diagnosis (Dx), post-treatment with a hypomethylating agent (Post-HMA), and progression to acute myeloid leukemia (AML). Bottom, boxplot showing the absolute cell-to-cell difference in epimutation rate (ΔPDR) for each condition. D Top, boxplots for responder and non-responder patients showing the PDR and ΔPDR for each cell and condition (hBM, Dx, Post-HMA, and AML). Bottom, boxplots showing ΔPDR per condition. Comparisons were performed using pairwise Wilcoxon tests, and the depicted P values were adjusted for multiple testing. We also found that post-HMA MDS samples presented higher PDR and ΔPDR compared to the Dx ones for all cases (Fig. 1C ), as well as when we separated patients into responders and non-responders (Fig. 1D and Supplementary Fig. S1 ). Thus, our results suggest that the demethylating drug is associated with increased epigenetic heterogeneity, but it is not related to clinical efficacy. However, it is noteworthy that the gene NPHP4 was differentially methylated between responders and non-responders, in line with previous findings that correlate the methylation status of this gene with the clinical response to HMAs [ 3 ]. In turn, comparing ΔPDR levels between responder and non-responder patients upon HMA therapy, we detected a significant increase in responders (Supplementary Fig. S2 ), suggesting that the more effective the HMA is, the greater the epigenetic heterogeneity. Consequently, even when PDR values are similar on average, the absolute difference in PDR values between cells (ΔPDR) is higher in responder patients. Importantly, we observed that epimutation rate was similar for post-HMA vs sAML regardless of the HMA response (Fig. 1D ), but only in all non-responder patients (Supplementary Fig. S1 ). Furthermore, the ΔPDR was reduced in AML samples compared to post-HMA ones, considering all cases and when divided by HMA response (Fig. 1C, D ). Nevertheless, this decrease in ΔPDR levels was only observed in three patients (including responders, #3 and #28, and non-responders, #8) after evaluating each patient individually (Supplementary Fig. S1 ), depicting different mechanisms driving MDS-to-AML transition at the single-cell epigenetic level. Thus, while differences are observed at the group level, individual cases highlight the extent of inter-patient heterogeneity. Projecting single-cell DNA methylation data onto lineage trees revealed distinct epimutation trajectories (Fig. 2A, B ). Two HMA-responder patients (#3 and #28) exhibited DNA methylation clades corresponding to the MDS and sAML diagnosis samples (Fig. 2A ), which suggests that there are differential DNA methylation patterns associated with each disease stage. However, the heterogeneous distribution of cells displayed by patient #26, who also responded to HMAs (Fig. 2A ), confirmed the epigenetic variability of our patient cohort and emphasized the need to identify the key factors that explain this diversity in larger cohorts. Strikingly, we observed the presence of single mSCs from the post-HMA stage in the clade enriched in sAML cells in refractory patients (Fig. 2B ), showing similar DNA methylation patterns. This could mean that these post-HMA mSCs may be selected during treatment, mediating therapy resistance and disease progression. Further investigations are thus required to determine the molecular mechanisms leading to the epigenetic evolution of these mSCs upon therapy and their role in promoting the MDS-to-AML transformation. Fig. 2: Epigenetic phylogenies and cellular composition dynamics in HMA-responder and non-responder patients. Full size image Epigenetic phylogenies corresponding to the five patients studied (responders, A , and non-responders, B ), where cells are colored according to the time point: at MDS diagnosis (Dx), post-treatment with a hypomethylating agent (Post-HMA), and progression to acute myeloid leukemia (AML). The color code shown in ( A ) applies to all phylogenies. The scale bar indicates the amount of epigenetic change it represents. Consequently, the length of each branch is proportional to the amount of epigenetic change indicated by the scale bar. C UMAPs representing the distributions of BM populations defined by scRNA-seq at MDS diagnosis and upon HMA treatment. Dashed circles represent HSC/MPP, CMP and GMP cell populations. CMP common myeloid progenitors, CNLOH copy-neutral loss of heterozygosity, DC dendritic cells, double-negative thymocytes CD4− CD8− double-negative progenitors in the thymus, ETP early T-cell precursors, GMP granulocyte-macrophage progenitors, HSC/MPP hematopoietic stem cells/multipotent progenitors, ILC3 Group 3 innate lymphoid cells, MAIT mucosal-associated invariant T cells, MK megakaryocytes, pDC plasmacytoid dendritic cells. To complement the results described above, three additional single-cell multi-omics layers (RNA expression, targeted DNA sequencing and protein) were analyzed in these patients (Fig. 1A ). Among the median of 3.81 million CpG sites (Interquartile Range, IQR: 2.61-5.21 million) analyzed per sample, the genomic distribution of the identified DNA methylation sites (among promoters, enhancers, gene bodies and intergenic regions) is shown in Supplementary Fig. S3 . Taking advantage of the fact that the Smart-RRBS technique enables to match DNA methylation with expression data for the same cells [ 7 , 8 ], we weighted paired gene promoter methylation and gene expression values according to single-cell coverage using a pseudo-bulk strategy ( Supplementary Methods ). We identified 27 genes in which DNA methylation correlated with expression levels (adjusted P value < 0.05) (Dataset S1 ). Interestingly, among the top candidates, we found STAT5B , which plays a critical role in myelopoiesis [ 11 ]. Interrogating scRNA-seq data from 104,919 cells by the Chromium 10x technology assigned 35 cell populations (Fig. 2C ). Case-by-case data are depicted in Supplementary Fig. S4 . We observed that HMA-responders underwent a depletion in the mSC content, mainly located in the stem cell (HSC/MPP) and myeloid progenitor (CMP and GMP) compartments, at the post-HMA time point (Fig. 2C ), which is in line with previous findings that demonstrated the reduction of the number of total and mutated cells within the HSPC compartment in responders [ 12 ]. However, it should be noted that it was primarily driven by a single patient (#3, Supplementary Fig. S4 ). On the other hand, we did not find any changes in mSC levels in non-responders at the same time point (Fig. 2C and Supplementary Fig. S4 ). Since we found a strong correlation between single-cell expression data derived from Smart-RBBS and 10x data (Spearman ρ = 0.864, adjusted P value < 0.001), we used this last approach for the subsequent analysis, as it allows a high-throughput profiling of thousands of cells, can identify rare cell types and provides a more robust clustering [ 13 ]. To go in-depth on the differential expression between responder and non-responder patients, both before and after HMA treatment, we focused on differentially regulated biological signaling pathways. In responders, HMA therapy was associated with the upregulation of cell-cycle checkpoints and DNA repair pathways, probably due to the accumulation of DNA damage caused by this treatment, and downregulation of FLT3 and MAPK-related signaling, frequently linked to poor outcomes in AML (Dataset S2A ). In addition, at MDS diagnosis, responders showed enrichment of signaling and immunity-related pathways (Dataset S2B ). We further evaluated the expression of two different signatures in the studied stem cell and myeloid progenitor compartments: one linked to AML-derived HSC/Progenitor-like cells at single-cell level [ 14 ] and another one associated with poor outcomes in myeloid malignancies at bulk RNA-seq resolution [ 15 ] (Dataset S3 ). We observed that both signatures were significantly underrepresented after HMA treatment in responders (Supplementary Fig. S5 ). These data match the observed depletion of CD34+ HSPCs in MDS cases that respond upon HMA treatment [ 4 , 12 ]. We also took advantage of the Tapestri technology to perform simultaneous scDNA-seq to study 53 genes mutated in myeloid malignancies (Dataset S4A ) and scProt-seq of 42 cell-surface proteins (Dataset S4B ) [ 12 ] in 17,849 cells from the above-described triplets of samples. Case-by-case and mutational data are described in Supplementary Fig S6 and Dataset S4C . We identified distinct genetic variants in 12 genes evenly distributed across the main altered biological processes in myeloid neoplasms: epigenetic/chromatin modifiers ( EZH2 , TET2 , IDH1 , ASXL1 ), transcription ( RUNX1 , BCOR ), RNA splicing ( SF3B3 , U2AF1 ) and signal transduction ( NRAS , PTPN11 , FLT3 , NPM1 ). Among these, initiating events involved alterations in epigenetic regulators (80% of the cases), and ASXL1 and TET2 were the genes harboring the highest number of variants (3 and 4, respectively) (Supplementary Fig. S6 ). In general, these initiating events were followed by mutations that affected transcription and RNA splicing genes, whereas late events were associated with signaling variants (primarily in NRAS and FLT3 ), adhering to overarching “rules” [ 10 ]. We did not observe any association between the type of epigenetic gene (or genes) mutated and the PDR and ΔPDR values determined by our single-cell DNA methylation approach. The mutational architecture analyses of these patients related to the HMA response are shown in Supplementary Fig. S6 . HMA-responders experienced a remarkable reduction of mutant mSC clones at the post-HMA sample vs the Dx time point, in line with previous reports [ 12 ], especially in patients #3 and #28 (Supplementary Fig. S6 ), whereas the derived sAML samples underwent a persistence of mutant clones, expansion of initially minor clones at Dx, and/or the irruption of new mutant mSC clones, depending on the patient (Supplementary Fig. S6 ). One interesting case was patient #3 (Supplementary Fig. S6 ), who acquired a FLT3 mutation during AML progression. Since FLT3-related pathways were transcriptionally downregulated in responders (Dataset S2A ), the emergence of the FLT3 mutations may represent a genetic mechanism enabling mSCs to overcome non-genetic constraints for transformation. Finally, once the HSPC compartment (CD34+ CD117+ cells) and mSCs (CD34+ CD45RA+ and CD34+ CD123+ cells) were defined using scProt-seq data (Supplementary Fig. S7A ), we observed that both the amount of mSCs and their number of mutations decreased considerably after HMA therapy, being restored after AML transformation, in responders. Nevertheless, the number of mSCs with more than 3-4 mutations remained constant in non-responders upon treatment and malignant evolution (Supplementary Fig. S7B ). It is also worth noting that scDNA-seq data reveal substantial genomic heterogeneity within both responders and non-responders, likely further contributing to variability in methylation and transcriptional responses. Overall, our data suggest that both genetic and epigenetic paths foster the transformation from MDS to AML through aberrant DNA methylation profiles [ 2 , 3 , 4 ] and complex clonal architectures [ 9 , 10 , 12 ]. These findings support the concept that multiple mSC populations co-evolve in parallel during MDS-to-AML progression [ 10 ], shaping a dynamic and competitive landscape rather than a linear evolutionary path. While we describe several of these interactions, the marked heterogeneity, captured by single-cell technologies, highlights that current models likely illustrate only part of the underlying complexity. The small sample size is a limitation of our study and a comprehensive validation in larger, deeply annotated cohorts to firmly establish the biological and clinical consequences for diagnosis, prognostic refinement and patient stratification is required. In this context, using novel cutting-edge single-cell multi-omics approaches would help us to understand how genetic and epigenetic alterations co-evolve, converge and establish interdependent evolutionary trajectories, ultimately uncovering actionable vulnerabilities. A more granular characterization of epigenetic heterogeneity in MDS, both at diagnosis and following HMA exposure, holds translational potential. Identifying and therapeutically targeting mSC subsets that share convergent DNA methylation states, potentially shaped or selected by treatment, could enable precision strategies aimed not merely at disease control but at intercepting clonal evolution and preventing AML transformation altogether.