Identifying Clones in Myelodysplastic Syndromes by Using Single‐Cell RNA Sequencing

Myelodysplastic syndromes (MDS) are heterogeneous myeloid lineage neoplasms characterized by cytopenia, recurrent infections, and excessive bleeding [1, 2]. MDS prevalence increases with age [1], and the incidence in men is almost twice that in women [2]. MDS are clonal disorders arising from a rare hematopoietic stem cell (HSC) population [3, 4] defined phenotypically as (Lin−)CD34+CD38−CD90+CD45RA−CD49f+ [5-7]. Although it is known that HSCs play crucial roles in both normal hematopoiesis [5] and MDS [6], a comprehensive understanding of gene expression patterns in MDS-derived HSCs is lacking. We hypothesized that some of the aberrant HSCs of MDS patients may differ in transcriptional signature from control HSCs. To test our hypothesis, we recruited 11 subjects suspected of MDS. As part of the routine clinical MDS workup, bone marrow (BM) aspirations were collected. From these, 1150 Lin−CD34+CD38−CD90+CD45RA−CD49f+ HSCs were sorted and subjected to single-cell RNA sequencing (scRNA-seq). Karyotype analysis was performed for six subjects. Further methodological details regarding the study population, cell sorting, library preparation, and scRNA-seq analysis are given in the Supporting Information S1. Consistent with previous studies [1, 2], the median age in the MDS+ group was 78 years, and there was a male predominance, compared to a median age of 76.5 and a female predominance in the MDS– controls (Table 1). The groups did not differ significantly with respect to median blood counts, incidence of moderate/severe anemia (hemoglobin < 10 mg/dL), or incidence of bicytopenia. In MDS– subjects, no other BM-related pathologies were detected. Most sequenced HSCs (1101/1150, 95.74%) passed quality control (Supporting Information S1: Figure S1). Most of these 1101 HSCs comprised one large cluster (cluster 0, Figure 1A) that did not express a unique transcriptomic signature and did not enable detection of abnormal clones (Supporting Information S1: Figure S2A). The 69 cells (6.2%) that comprised cluster 1 highly expressed S and G2/M cell-cycle markers (Supporting Information S1: Figure S2B) and we thus concluded that these cells were proliferating. Importantly, the fraction of proliferating cells was not significantly affected by MDS status [5.23% (37/707) of MDS+-derived cells versus 8.12% (32/394) of MDS–-derived cells; p_val = 0.06, χ2 = 3.59]. Erythroid markers were highly expressed in the six cells (0.54%) comprising cluster 2, suggesting erythroid differentiation, and those cells were thus removed from the analysis. The above results, taken together, show that transcriptional state alone did not enable detection of potentially pathological HSC clones in MDS. Gene expression patterns of HSCs derived from MDS patients and controls. (A) UMAP projection displaying 1101 HSCs derived from seven MDS+ patients and four MDS– controls. (B) Heatmap depicting the smoothed expression level of 14 chrY genes (X-axis) in 733 cells derived from six MDS+ males and one MDS– male control. Genes are ordered by their chromosomal location, and cells (Y-axis) are sorted by subject. (C) Histogram depicting the number of cells as a function of the mean normalized expression level of chrY genes. Distribution curves for each subject are indicated. Colors of subjects are as in (B). The red distribution is subject b (MDS–), used as reference. Vertical red line marks mean chrY expression value of 0.92. (D) UMAP projection (same as in A) colored by chrY status in cells derived from subjects c and d. Red—del(Y)-HSC (N = 177); turquoise—cells without chrY deletion (N = 70). (E) Volcano plot depicting differentially expressed genes (DEGs) between del(Y)-HSCs and chromosomally intact HSCs. The Wilcoxon test Bonferroni corrected adjusted p-value is presented. The horizontal dashed line indicates an adjusted_p_value cutoff of 0.05. Vertical lines indicate a log2(average fold_change) (FC) of ± 1. Genes indicated are significantly differentially expressed (adjusted_p_value < 0.05) and have a log2(average FC) > |1|. The genes KDM5D, mapped to chrY, and ANXA1 mapped to chr9, are shown in red boxes. (F) Bar plot depicting results of pathway enrichment analysis of DEGs between del(Y)-HSCs and chromosomally intact HSCs from subjects c and d. The top four hallmark gene sets are shown on the Y-axis, and the change in the activity of each pathway using the qval scale is shown on the X-axis (for further information please refer to Supporting Information S1). The number on each bar is the FC of the specific pathway. All presented pathways have an adjusted_p_val < 0.05. Panels were grouped using Inkscape. We then tested whether inferred copy number aberrations (CNAs) could be used to define clonal HSCs. No autosomal chromosomal changes were detected in MDS+ subjects (Supporting Information S1: Figure S2C), in accordance with clinical karyotype analysis data. Intriguingly, focusing on the Y chromosome, we showed that 239 of the 733 male-derived HSCs exhibited lower expression levels of chromosome Y (chrY) genes, defined as avg . chrY . expression $\text{avg}.\,\text{chrY}.\,\text{expression}$ < 0.92. This finding suggested chrY deletion [del(Y)-HSC, Figure 1B]. Two MDS subjects (c and d) displayed a higher proportion of del(Y)-HSCs than control subjects (Figure 1C), suggesting that those subjects carried mosaic loss of chrY (mLOY). Of the HSCs derived from subject c, 44% (50/113) were del(Y)-HSCs, consistent with the clinical karyotype that identified 4/10 (40%) of BM cells bearing the chrY deletion. Notably, HSCs obtained from subjects c and d were not distinctly separated according to their chrY status, suggesting minor transcriptional differences in chrY deletion (Figure 1D). To reveal genes potentially relevant to the pathogenesis of MDS, we then searched for differentially expressed genes between 177 inferred del(Y)-HSCs and 77 chrY intact HSCs derived from subjects c and d. Of the 41 differentially expressed non-HLA genes [Bonferroni adjusted p_val < 0.05 and log2(average fold_change) >|1|; Supporting Information S1: Table S1], 35 were downregulated in the del(Y)-HSCs (Figure 1E). As expected, all eight chrY genes were significantly downregulated in del(Y)-HSCs. Those genes included KDM5D (of chrY), which codes for lysine-specific demethylase 5D protein, and was previously reported to be involved in leukemogenesis [8]. Notably, the ANXA1 gene (of chr9), encoding for annexin A1, an immunomodulatory protein, was downregulated in del(Y)-HSCs. Potential involvement of ANXA1 in the pathogenesis of MDS is consistent with a previous study that reported reduced ANXA1 in MDS, especially in those patients progressing to leukemia [9]. Pathway analysis showed that the HALLMARK_APOPTOSIS pathway, which includes ANXA1, was downregulated in del(Y)-HSCs versus chromosomally intact cells derived from the same subjects (adj_p_val < 0.0001, fold_change = −2.41, qval = 2.788, Figure 1F). Although it has previously been reported that upregulation of BCL2 reduces apoptosis in MDS [4], BCL2 was not differentially expressed in our analysis. Taken together, our findings identified aberrant HSC clones bearing a chrY deletion that showed decreased apoptosis with ANXA1 downregulation. Due to chrY involvement, we performed differential expression analysis between sexes. We identified 17 downregulated and 30 upregulated genes in female- and male-derived HSCs, respectively [Bonferroni adjusted p_val < 0.05 and log2(average fold_change) > |1|; Supporting Information S1: Table S2]. The four genes that mapped to chromosome X and the 13 that mapped to chrY were upregulated in female- and male-derived HSCs, respectively. In summary, we investigated the role of rare HSCs in MDS patients at the single-cell resolution. We found chromosomal abnormalities, specifically chrY deletion, in accordance with previous findings [6]. Two of the genes that were downregulated in MDS del(Y)-HSCs of two male patients, KDM5D and ANXA1, were previously reported to be involved in MDS and leukemia [8, 10]. Additionally, MDS del(Y)-HSCs downregulated apoptosis without a statistically significant upregulation of the BCL2 gene. Despite the small numbers of participants and sequenced cells, and the sex imbalance between MDS and control groups, our results suggest that in del(Y)-HSCs, KDM5D and ANXA1 are markers for and potential drivers of MDS. U.G., R.G., R.P. and T.S. conceived the study. U.G. and I.L. selected and recruited patients. L.O. and R.G. handled specimen enrichment, fluorescence-activated cell sorting and single cell sorting. H.N.-G., R.P. and T.S. analyzed the data. R.P. and T.S. wrote the manuscript. All authors reviewed the results and approved the final version of the manuscript. This study was supported by a multidisciplinary internal grant from Ben-Gurion University of the Negev to U.G., R.G., and T.S., and the Israel Science Foundation (ISF) grants no. ISF883/21 and ISF458/21 to R.G. and T.S., respectively. The funders had no role in the design, data collection, data analysis, and reporting of this study. This study was approved by the Institutional Helsinki Committee of Soroka University Medical Center (Institutional approval No. SOR-0115-17). All study participants signed an informed consent form. The authors declare no conflicts of interest. The data that support the findings of this study are openly available in NCBI Gene Expression Omnibus (GEO) at https://www.ncbi.nlm.nih.gov/geo/, reference number GSE293014. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.