Analysis of Myelodysplastic Syndrome Stem Cells at Single Cell Resolution during DNA Methyltransferase Inhibitor Therapy

The myelodysplastic syndromes (MDS) arise in and are maintained by hematopoietic stem cells (HSCs). Serial sampling of patients treated with DNA methyltransferase inhibitors (DNMTIs) and lenalidomide has demonstrated that disease HSCs (MDS HSCs) persist at significant levels even in patients achieving complete clinical and cytogenetic responses. As MDS HSCs are the functional unit of clonal selection both during therapy and subsequent disease progression, we hypothesized that the molecular heterogeneity of MDS HSCs may underlie therapeutic resistance. We therefore sought to perform single cell RNA-sequencing (RNA-seq) on MDS HSCs from patients with known responses to therapy, with the intention of identifying novel therapeutic vulnerabilities. To characterize MDS HSC heterogeneity, we FACS-purified HSCs (Lin-CD34+CD38-CD90+CD45RA-) from paired bone marrow (BM) specimens taken from four MDS patients before and after two to four 28-day cycles of the DNMTI decitabine, as well as two patients who were not treated due to stable disease, and two normal age matched controls. Specimens from both responding and non-responding patients were included. We captured and sequenced a total of 869 single cells from 14 samples, sequencing to an average depth of 4.8 million reads. In a subset of samples (n=7) we also performed bulk RNA-seq (average 1500 cells) for comparison. The sequencing data was of high quality, with an average of 80% mapped reads. We confirmed our ability to accurately quantify transcript levels using ERCC spike-in controls, observing a linear correlation between expected concentration and observed FPKM (fragments per kilobase per million). Single cell RNA-seq revealed vast intratumoral heterogeneity in MDS HSCs that was otherwise missed by bulk RNA-seq, as evidenced by the presence of transcripts variably expressed among cells from the same specimen (Fig. 1A). Despite this intratumoral heterogeneity, single cell transcriptomes were able to completely separate individual MDS patients using principal components analysis and hierarchical clustering, consistent with the known heterogeneity of MDS. MDS HSCs further clustered separately from normal age-matched HSCs, with the top 10% of genes contributing to this separation enriched for Gene Ontology (GO) categories including pathways implicated in MDS biology such as mRNA splicing, nonsense mediated decay, and P53 mediated DNA Damage Response (all P Unsupervised hierarchical clustering of all pre-treatment MDS HSCs revealed clustering of cells from responders separately from non-responders (Fig. 1B). Differential gene expression analysis identified a cluster of genes (FDR Finally, within each sample we measured the spread of gene expression using dispersion (log[variance/mean]) within bins based on expression levels, defining variable genes as those with a dispersion >1.75 at a mean FPKM >2. The highest number of variable genes were in normal HSCs (mean=141), with the next highest in responders prior to treatment (mean=80), and the least number of variable genes in non-responders prior to treatment (mean=9.5). We speculate that the low number of variable genes in non-responders reflects a higher degree of clonal dominance. All post-treatment MDS HSCs demonstrated a relatively low number of variable genes (mean=25), suggesting that therapy induces clonal selection. In sum, our data illustrate the robustness of single cell RNA-seq to define the intrinsic variability of individual MDS HSCs, implicating perturbed ribosomal biogenesis and transcriptional variability as novel predictors of response to therapy. As we expand our data set with additional patients, we expect to identify additional pathways that mediate therapeutic response and resistance, as well as mutations and variably expressed genes that are selected for during therapy and drive disease progression. Disclosures No relevant conflicts of interest to declare.