Single-cell multiomics analysis of chronic myeloid leukemia links cellular heterogeneity to therapy response
eLife assessment
This study presents fundamental insights into the heterogeneity of chronic myeloid leukemia (CML) stem cells and their response to tyrosine kinase inhibitor therapy, shedding light on potential mechanisms underlying treatment failure. The study's robust methodology, supported by validation with bulk RNA-seq data and surface marker analysis, provides compelling evidence for the identified associations between cellular composition and treatment outcome. These findings contribute to our understanding of CML pathogenesis and may inform the development of more targeted therapeutic strategies.
https://doi.org/10.7554/eLife.92074.3.sa0Fundamental: Findings that substantially advance our understanding of major research questions
- Landmark
- Fundamental
- Important
- Valuable
- Useful
Compelling: Evidence that features methods, data and analyses more rigorous than the current state-of-the-art
- Exceptional
- Compelling
- Convincing
- Solid
- Incomplete
- Inadequate
During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments
Abstract
The advent of tyrosine kinase inhibitors (TKIs) as treatment of chronic myeloid leukemia (CML) is a paradigm in molecularly targeted cancer therapy. Nonetheless, TKI-insensitive leukemia stem cells (LSCs) persist in most patients even after years of treatment and are imperative for disease progression as well as recurrence during treatment-free remission (TFR). Here, we have generated high-resolution single-cell multiomics maps from CML patients at diagnosis, retrospectively stratified by BCR::ABL1IS (%) following 12 months of TKI therapy. Simultaneous measurement of global gene expression profiles together with >40 surface markers from the same cells revealed that each patient harbored a unique composition of stem and progenitor cells at diagnosis. The patients with treatment failure after 12 months of therapy had a markedly higher abundance of molecularly defined primitive cells at diagnosis compared to the optimal responders. The multiomic feature landscape enabled visualization of the primitive fraction as a mixture of molecularly distinct BCR::ABL1+ LSCs and BCR::ABL1-hematopoietic stem cells (HSCs) in variable ratio across patients, and guided their prospective isolation by a combination of CD26 and CD35 cell surface markers. We for the first time show that BCR::ABL1+ LSCs and BCR::ABL1- HSCs can be distinctly separated as CD26+CD35- and CD26-CD35+, respectively. In addition, we found the ratio of LSC/HSC to be higher in patients with prospective treatment failure compared to optimal responders, at diagnosis as well as following 3 months of TKI therapy. Collectively, this data builds a framework for understanding therapy response and adapting treatment by devising strategies to extinguish or suppress TKI-insensitive LSCs.
Introduction
The persistence, burden, and evolvability of LSCs during therapy provide common threads to address the major challenges remaining in CML (Cortes et al., 2021; Michor et al., 2005). For example, the acquisition of the BCR::ABL1 kinase domain (KD), and additional somatic mutations within LSCs could either necessitate a switch to another TKI or precipitate a blast crisis (Bolton-Gillespie et al., 2013; Giustacchini et al., 2017; Nieborowska-Skorska et al., 2012; O’Hare et al., 2007; Sorel et al., 2004). Moreover, even in the absence of KD mutations, the LSC burden at diagnosis is linked to cytogenetic and molecular response in the chronic phase (Khorashad et al., 2006; Mustjoki et al., 2013), and their presence is associated with recurrence during TFR (Pagani et al., 2020), thus compelling life-long treatment in the majority of patients.
The de facto standard for disease monitoring, BCR::ABL1 transcript/gDNA quantification using qPCR, Sanger, and/or Next Generation Sequencing (NGS), however, predictably only captures cells positive for BCR::ABL1 expression and translocation and does not inform the number and identity of fully leukemogenic cells (Radich et al., 2019). This is pertinent given that a combination of cell sorting and BCR::ABL1 quantification suggested that the presence of BCR::ABL1 signal in the stem cell compartment rather than the lymphocytes correlates with TFR loss (Kinstrie et al., 2020; Pagani et al., 2020). Using HSC specific markers, however, is confounded by the observations that in remission, patients in chronic phase restore Philadelphia chromosome-negative (Ph-) hematopoiesis (Bergamaschi et al., 1994; Carella et al., 1993; Claxton et al., 1992; Coulombel et al., 1983; Deininger, 2003), implying the co-existence of BCR::ABL1- HSCs and BCR::ABL1+ LSCs. However, because their proportion varies across individuals and treatment stages (Mustjoki et al., 2013), and both LSCs and HSCs reside within similar immunophenotypic compartments (Lin-CD34+CD38-/low) (Eisterer et al., 2005; Petzer et al., 1996), discrimination between these populations – and thereby reliably estimating LSCs – remains difficult. As a result, understanding CML LSC identity, heterogeneity, and vulnerability to TKI therapy remain outstanding challenges in CML.
We recently implemented BCR::ABL1 targeted single-cell RT-qPCR in combination with index sorting for surface markers to dissect the LSC compartment at diagnosis and following 3 months of TKI therapy (Warfvinge et al., 2017). This demonstrated that at diagnosis, Lin-CD34+CD38-/low BCR::ABL1+ cells could be divided into seven molecularly distinct subpopulations, each reflecting unique lineage and cell state signatures. Importantly, the proportions of BCR::ABL1+ LSC subpopulations changed dynamically upon therapy and only a few surface markers efficiently discriminated between BCR::ABL1+ vs. BCR::ABL1- cells. Of these, despite being generally reduced in BCR::ABL1+ cells upon therapy, CD26 was more frequently detected in the subpopulation most persistent during TKI treatment. The substantial treatment-induced changes observed within the stem cell population suggested that BCR::ABL1+ LSCs are themselves heterogeneous in terms of sensitivity to TKIs allowing us to define BCR::ABL1+Lin-CD34+CD38-/lowCD45RA-cKIT-CD26+ as the most TKI-insensitive cells. Apart from being detected at diagnosis, these cells are strikingly enriched after TKI therapy (Warfvinge et al., 2017), a finding since confirmed in a transgenic mouse model (Shah et al., 2023), and in clinical observations documenting the long-term persistence of CD26+ cells in CML patients (Pacelli et al., 2023). Together, these findings argue that functional heterogeneity within LSCs cannot be predicted solely by surface markers but is intimately linked to their cell state and gene expression signature, thereby motivating a strategy to simultaneously capture their lineage potency/affiliation, BCR::ABL1 status, and molecular program in addition to the surface markers.
Although recent single-cell omics studies have begun to explore the cellular landscape in CML (Giustacchini et al., 2017; Krishnan et al., 2023; Patel et al., 2022; Zhang et al., 2020), these have largely been limited to measuring single modalities at a time, and thus, multiomic investigation of BCR::ABL1+ LSCs and BCR::ABL1- stem cells from the same patient are currently lacking. Therefore, we used cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq Stoeckius et al., 2017) to quantify the global gene expression program, and highly multiplexed cell-surface protein profile simultaneously within the same cells from nine CML patients. The single-cell multiomics maps enabled us to visualize the cellular heterogeneity of each patient by delineating the abundance of molecularly defined leukemic stem versus lineage progenitors, and link cell composition at diagnosis to TKI therapy response. Moreover, the maps highlighted the co-existence of molecularly distinct BCR::ABL1+ vs. BCR::ABL1- stem cells, guiding their prospective isolation by novel combinations of CD26 and CD35 markers, dissection of imbued molecular programs, division kinetics and changes in heterogeneity following TKI therapy.
Results
Single-cell multiomic CITE-seq analysis defines the heterogeneity of CML stem cell fractions
To characterize the multimodal heterogeneity within the Lin-CD34+ and Lin-CD34+CD38-/low CML bone marrow (BM) compartment at diagnosis and estimate how the heterogeneity contributes to therapy response, we applied CITE-seq on nine patients retrospectively stratified according to molecular response following 12 months of treatment. To allow for accurate specification of LSPC heterogeneity without confounding effects generated from cell-type-irrelevant transcriptional signals related to BCR-ABL transformation per se, or inter-patient variability due to batch effects, we used Scarf projection (Dhapola et al., 2022) to map each individual CML cell to a Lin-CD34+ healthy bone marrow reference with a predefined heterogeneity. Each CML cell was then annotated by label transfer according to the annotation of the reference cells with the highest molecular similarity (Figure 1A).
The CITE-seq analysis combines single-cell RNA-seq with simultaneous staining for a panel of >40 antibodies conjugated to uniquely barcoded oligonucleotides (henceforth Antibody derived tags, ADTs) aimed to cover the spectrum of normal hematopoietic stem (CD34+CD38-/low fraction) and lineage progenitors (CD34+ enriched fraction, HSPCs) as well as surface markers reported to enrich for CML stem cells (Herrmann et al., 2014; Järås et al., 2010; Kinstrie et al., 2020; Landberg et al., 2018; Sadovnik et al., 2016; Warfvinge et al., 2017). Upon sequencing, the scRNA-seq and ADT-seq libraries provide a concurrent readout of gene expression and abundance of surface proteins, respectively. This multiomic approach surpasses scRNA-seq in cellular subtype characterization and is not constrained by the conundrum of spectral overlap of fluorophores in FACS Stoeckius et al., 2017, and availability of metal isotopes in mass cytometry as noted previously (Gullaksen et al., 2019). Our sample set of nine CML patients are from various clinical studies with attendant information on BCR::ABL1 transcript level as per the International Scale (IS %) available after 12 months of TKI therapy. These were retrospectively classified as optimal, treatment failures, and warning cases in accordance with the guidelines from European LeukemiaNet (ELN) (Hochhaus et al., 2020 Supplementary file 1).
To build a normal reference, we used our recently described CITE-seq profiles of Lin-CD34+ cells from normal BM from two age-matched healthy donors (henceforth, nBM) (Sommarin et al., 2021). The normal reference of 4696 Lin-CD34+ cells was visualized in a Uniform Manifold Approximation and Projection (UMAP, Figure 1B). Although single-cell omics datasets generate expression values for several thousand genes (dimensions), not all are equally informative. UMAP is a non-linear algorithm for dimensionality reduction and preserves the global data structure and enables visualization of the dataset in two dimensions. The nBM cells were subsequently clustered and annotated based on marker genes as well as information from antibody-derived tags (ADTs) against surface markers (Figure 1C–D, Supplementary file 2). We identified eleven cell clusters consisting of primitive cells, lineage progenitors and cycling cells in the UMAP, with near-to-equal distribution across the two nBM samples (Figure 1—figure supplement 1). The primitive cluster was marked by expression of stem-like gene expression e.g., CRHBP, HLF as well as a CD90+CD38-/low surface profile. Closest to the primitive cluster were MPP1 and 2 without frank manifestation of lineage-specific markers. The MPP1 cluster, in turn, was connected to the clusters with gain of myeloid-like, and lymphoid-like gene expression programs. This included myelo/lymphoid (My/Ly), myeloid (My: MPO, CD33), and lymphoid/dendritic/monocyte (Ly/pDc/mono) clusters (Ly86, CD135). The MPP2 cluster, in comparison, was connected to the megakaryocyte and erythroid lineage progenitors such as megakaryocyte-erythroid (MEP), megakaryocyte (MkP: ITGA2B, CD41), erythroid progenitors (ErP: HBA1, CD71, CD105) and basophil/mast cell clusters (Baso/mc: HDC, MS4A2-3, ITGB7).
Following CITE-seq of 14,274 Lin-CD34+ cells from nine CML patients, the cells were projected onto the nBM reference and annotated by label-transfer. Subsequently, we generated UMAPs for each patient and color-coded cells according to the transferred labels enabling us to visualize and highlight cellular heterogeneity on an individual basis (Figure 2A, Figure 3—figure supplement 1). A comparison of abundance of clusters showed that CML patients, in general, were enriched for basophil/mast cell, MEP, MkP clusters and depleted for primitive, MPP1, 2, and ly/pDC/Mono clusters (Figure 2B). Inspection of profiles per patient showed that while most of the clusters were consistently enriched or depleted, ErP progenitors and cycling cells showed highly individual patterns (Figure 3—figure supplement 2B). Furthermore, visual inspection of the individual UMAPs revealed an intriguing and distinct separation of cells within the primitive cluster in most of the patients (Figure 2A). Taken together, CITE-seq and label transfer through nBM reference cell projection could successfully define the molecular heterogeneity within the Lin-CD34+ LSPC at diagnosis and revealed an average increase in abundance of MEP, MkP, and basophil/mast cells as compared to nBM.
Identification of gene signatures unique to CML cells and shared across CML stem and progenitor populations
The identification of nBM and CML primitive populations and lineage progenitors prompted us to systematically compare pair-wise gene expression differences between all eleven clusters in CML and their nBM counterparts (Figure 2C, Supplementary file 5). For each pair-wise comparison, we highlighted only the genes that were differentially expressed between the pair but not shared with any other cluster thereby defining strict cluster-specific changes. Although this revealed unique changes across all cluster pairs, the most striking changes were seen between primitive clusters from CML versus nBM with 384 differentially expressed genes. Moreover, erythroid progenitors also displayed large changes followed by ly/pDC/mono and ErP progenitors.
In contrast, when considering genes differentially expressed in CML relative to nBM, and crucially shared among all the cluster comparisons, we identified a set of 71 genes that comprised a pan-cml signature. This approach is preferrable to a ‘bulk’ comparison of all CML CD34+ cells lumped together versus nBM cells, where the most numerically abundant cluster, and the most highly expressed genes therein, are likely to contribute disproportionately (Figure 2D). Of the 71 genes, 50 genes were found to be consistently up-regulated in CML stem and progenitors. Notably, the individual genes from this set were expressed at varying levels in individual clusters (Figure 3—figure supplement 2A, Supplementary file 7), reflecting a heterogenous origin of the CML signature. This included a constellation of genes such as surface markers, CD81 (Quagliano et al., 2020), nucleoside diphosphate kinase with transcriptional regulatory activity, NME2 (Kar et al., 2012; Thakur et al., 2009; Thakur et al., 2014; Thakur et al., 2011; Tschiedel et al., 2012; Tschiedel et al., 2008; Yadav et al., 2014) previously implicated in CML, genes of histone family, and components of SWI/SNF chromatin remodeling complex, BRD7 (Kaeser et al., 2008), as well as genes of histone methyltransferase family, KMT3A and 5 A (Husmann and Gozani, 2019). The analysis also identified 21 genes specifically and consistently upregulated in all clusters across nBM versus CML, and among others included DPY30 gene implicated in the maintenance of adult HSCs (Yang et al., 2016).
Progenitor heterogeneity at diagnosis is linked to TKI therapy outcome
To assess how cellular heterogeneity within the Lin-CD34+ compartment at diagnosis is related to therapy response, the patients were stratified according to BCR::ABL1IS (%) following 12 months of TKI therapy as per ELN recommendations (Supplementary file 1). The stem and progenitor heterogeneity at diagnosis was then evaluated relative to the treatment outcome by assessing cell composition of the individual CML patients 1–9 classified as optimal (CML1-4), warning (CML 5–7), and treatment failure (CML 8–9).
By enumerating the proportion of molecularly defined clusters of individual patients, we found that all clusters from the normal bone marrow reference samples were represented in each patient. However, there was remarkable heterogeneity between patients in terms of relative abundance of constituent clusters (Figure 3A, Figure 3—figure supplement 1B). Importantly, this inter-patient heterogeneity followed patterns connected to therapy outcome. A focused comparison between optimal responders (CML 1–4) and treatment failure cases (CML 8–9) showed distinctive enrichment of specific cell clusters (Figure 3B), where failure cases displayed a higher abundance of the molecularly defined primitive cluster, multipotent progenitors 1 and 2, and myelo-lymphoid (My/Ly) and Ly/pDC/monocyte clusters at diagnosis. In contrast, the optimal responders had a higher burden of megakaryocyte-erythroid progenitors (MEP), megakaryocyte progenitors (MkP), Erythroid progenitor (ErP), Baso/mc, and cycling clusters at diagnosis. Among others, the most conspicuous changes between optimal responders and treatment failures were observed in the primitive cluster (fourfold higher in failures), and MEP and erythroid cluster (∼threefold higher in optimal). The warning cases (CML 5–7) also showed heterogeneity in terms of cell composition (Figure 3A); for instance, while CML 5 had a relatively higher content of primitive cells akin to treatment failures, CML 6 and 7 had a profile resembling optimal responders. Thus, these observations indicate that the prospective treatment outcome is reflected by the composition of stem and progenitor cell types at diagnosis.
Large-scale deconvolution of an independent dataset shows that cell composition is linked to cytogenetic response to TKI therapy
To further assess whether the hematopoietic stem and progenitor composition of CML patients at diagnosis is connected to their therapy response, we utilized an independent publicly available bulk gene expression dataset from CML patients with known therapy response (n=59 McWeeney et al., 2010). This was done with an aim to computationally deconvolute the individual profiles into constituent cell types by employing the same nBM used to define the heterogeneity in the CITE-seq data as a molecular reference, subsequently estimate their relative abundances, and relate to TKI response (Figure 3C). We used CIBERSORTx (Newman et al., 2019), a widely used algorithm, which has also been recently used to deconvolute bulk transcriptomes from acute myeloid leukemia (AML) (Zeng et al., 2022). As the gene expression profiles were measured at diagnosis from CD34+ enriched cells from either bone marrow or peripheral blood, this represented a cellular fraction broadly similar to ours for comparison. Importantly, the CML patient cohort was stratified by the percentage of Ph+ cells after 12 months of Imatinib treatment as responders (0% Ph+ metaphases, complete cytogenetic response, CcyR) or non-responders (>65% Ph+ metaphases, lack of even a minor cytogenetic response) as per the original study (McWeeney et al., 2010). We adhered to the patient annotation (responder or non-responder) provided by the original authors for consistency (see methods for more details).
The deconvolution results demonstrated a heterogeneous cell composition across patients, and overall comparison of imatinib non-responders (n=18 out of 59) versus responders (n=41 out of 59) showed a >threefold statistically significant enrichment of primitive cells in non-responders (Figure 3D, Figure 3—figure supplement 2B). To assess whether the source of origin of CD34+ cells in the public dataset could bias our results, we compared the abundances of imputed cell types from non-responders and responders based on their extraction from either bone marrow, or peripheral blood separately. This once again revealed a statistically enriched presence of primitive cells in non-responders versus the responders in the bone marrow (∼threefold higher). Moreover, and unlike the findings from overall comparison, myeloid cells were also found to be statistically enriched within responders (∼2.5 fold higher Figure 3E). Although MPP2 and MkP clusters showed increased abundance in non-responders and responders, respectively, the difference did not reach statistical significance. Notably, the primitive cells were rarely detected in peripheral blood in contrast to the bone marrow, confirming the notion that cell type fractions differ in different tissues (Figure 3D), and the bone marrow is predictably a more reliable and abundant source of true primitive cells in CML. Taken together, our analyses revealed distinct cellular heterogeneity across CML patients at diagnosis and highlighted that cell composition is associated with therapy response. Specifically, a higher burden of cells with a stem-like signature at diagnosis is a strong indicator of treatment failure.
Identification of BCR::ABL1+ and BCR::ABL1- primitive cells and their surface markers by combining single-cell gene expression with ADTs
Interestingly, according to the ADT data, the established CML stem cell markers CD25, CD26, and CD93 could only capture a fraction of the primitive cells (Figure 3—figure supplement 2C). To further investigate the molecular and immunophenotypic heterogeneity within the primitive cluster we performed CITE-seq on the stem cell enriched Lin-CD34+CD38-/low compartment from the same patients (n=8; 6779 cells Figure 4—figure supplement 1A and Figure 4—figure supplement 2A). A stepwise analysis (Figure 4A) was performed to resolve the BCR::ABL1 status of the individual primitive cells and determine their immunophenotype. First (1), the primitive cells from the Lin-CD34+CD38-/low fraction of all CML patients were merged and classified by BCR::ABL1 status using publicly available BCR::ABL1+ gene signatures. Next (2) CITE-seq data from both the Lin-CD34+ and the Lin-CD34+CD38-/low sorted population from the same patient was merged and visualized in a joint UMAP, while (3) the BCR::ABL1 status of the Lin-CD34+CD38-/low primitive cells was linked to the merged UMAPs by cell ID. Finally (4) the log2 fold change in ADT expression between BCR::ABL1+ and BCR::ABL1- primitive cells were calculated and the remaining cells of the UMAP were color-coded by the cluster identity generated by label-transferred projection to nBM.
The BCR::ABL1 status of the primitive cells was determined using two previously established gene expression signatures comparing Lin-CD34+CD38-/low BCR::ABL1+ LSCs to normal Lin-CD34+CD38-/low HSCs, and Lin-CD34+CD38-/low BCR::ABL1- non-leukemic stem cells (Giustacchini et al., 2017). These signatures have been derived by analyzing cells with SMART-seq for global transcriptome, and primers specific to the BCR::ABL1 fusion gene, simultaneously and, therefore, represent useful anchors to query and label cells as either BCR::ABL1+ or BCR::ABL1- (See methods for detailed analysis).
In the UMAP generated from all Lin-CD34+CD38-/low primitive cells, we identified three coarse clusters, where two consisted of primitive CML cells with representation from all patients (cluster 2–3), and one represented Lin-CD34+CD38-/low primitive cells from nBM which were included as negative control (cluster 1 Figure 4—figure supplement 1B–C). When we queried each cluster for enrichment of either LSC or HSC-like gene signatures the analysis clearly demonstrated that the primitive cells from CML samples could be divided based on the presence of BCR:ABL1, where CML cluster 3 displayed an evident enrichment of BCR::ABL1+ LSC specific signatures as compared to CML cluster 2 consisting of BCR::ABL1- non-leukemic stem cells (Figure 4—figure supplement 1D–E). Interestingly, both CML clusters were also found to share gene characteristics with normal stem cells (Figure 4—figure supplement 1D). Accordingly, we annotated Lin-CD34+CD38-/low primitive cells from individual patients as either BCR::ABL1+ or BCR::ABL1- and could observe that the distribution varied across patients, suggesting a variable load of leukemic stem and normal stem-like cells within the primitive fraction (Figure 4—figure supplement 1F, Figure 4—figure supplement 2A).
Remarkably, subsequent matching of BCR::ABL1+ and BCR::ABL1- cell IDs to the primitive cells in the joint Lin-CD34+ and Lin-CD34+CD38-/low UMAPs (indicated by red and black cell color, respectively), indeed revealed that BCR::ABL1+ LSCs and BCR::ABL1- HSCs were distinctly and consistently separated in the UMAP for each patient (Figure 4B, Figure 4—figure supplement 3). BCR::ABL1+ LSCs in all patients were positioned immediately juxtaposing the downstream progenitor populations forming a continuum of differentiation hierarchy characteristically observed in single-cell analysis of HSPCs (Pellin et al., 2019; Safi et al., 2022; Velten et al., 2017). In contrast, BCR::ABL1- HSCs formed an isolated and distinct cluster clearly detached from the rest of the cells in the majority of patients, indicating that during the chronic phase of CML, active hematopoiesis is dominated by the BCR::ABL1+ LSCs while BCR::ABL1- HSCs reside in the bone marrow albeit in an inactive state with reduced contribution to hematopoiesis as previously suggested (Chen et al., 2023; Coulombel et al., 1983). Moreover, using the ADT data we could validate the molecularly defined BCR::ABL1+ LSCs also by surface expression, with CD26, CD25, and CD93 displaying a 2–8 fold elevated expression in BCR::ABL1+ primitive cells as compared to BCR::ABL1-. A noteworthy patient-specific variability in surface protein expression on the most primitive cells was also established (Figure 4B).
To gain an understanding of the nature of molecular program of BCR::ABL1+ vs. BCR::ABL1- primitive cells, we analyzed differentially expressed genes between the pooled primitive BCR::ABL1+ and BCR::ABL1- populations from all CML patients. This identified 244 genes specifically upregulated in BCR::ABL1+ and 180 genes upregulated in BCR::ABL1- cells (Figure 4C, Supplementary file 6). Interestingly, among genes up-regulated in BCR::ABL1+ primitive cells included the expression of established CML markers; IL2RA (CD25), DPP4 (CD26), and IL1RAP as well as LEPR (CD295) a receptor recently described as a novel marker for Lin-CD34+CD38-/low LSCs (Landberg et al., 2018). In contrast, BCR::ABL1- cells displayed an up-regulated expression of the stem cell-related gene CRHBP (Barbieri et al., 2022; He et al., 2005) as well as CXCR4 important for homing and maintenance of normal HSCs (Sugiyama et al., 2006). Upon applying the gene signatures to Lin-CD34+ cells, we could clearly identify as well as partition primitive cells into BCR::ABL1+ vs. BCR::ABL1- (Figure 4D and Figure 4—figure supplement 3). This is especially relevant as the Lin-CD34+ fraction contains a higher proportion of lineage progenitors but a smaller fraction of primitive cells making their identification challenging. Taken together, this multistep analysis approach generated unique detailed UMAPs allowing for inspection of BCR::ABL1+ LSCs and BCR::ABL1- HSCs, their relation to other progenitor cells, as well as further validate BCR::ABL1+ LSCs identity by surface protein expression.
A unique combination of CD26 and CD35 surface markers captures molecularly defined BCR::ABL1+ LSCs and BCR::ABL1- HSCs
A long-sought aim in CML research has been to generate protocols for effective discrimination and isolation of CML LSCs from residual HSCs within individual patient bone marrow. The single-cell multiomic feature space of CITE-seq together with the identification of molecular distinct BCR::ABL1+ LSCs and BCR::ABL1- HSCs provides a unique opportunity to identify efficient protocols for their separation.
To specifically define the surface markers for their ability to capture BCR::ABL1+ vs. BCR::ABL1- primitive cells, we compared significantly up-regulated ADTs for either Lin-CD34+CD38-/low BCR::ABL1+ LSCs or BCR::ABL1- HSCs cells across patients and identified markers with consistent expression across the data set. We found 28 unique markers to be significantly up-regulated for either BCR::ABL1+ or BCR::ABL1- cells (black or red colored bars in Figure 4B, Figure 5A). Of these 28 markers, 17 ADTs marked BCR::ABL1+ LSCs and were detected in at least one patient. However, only CD26, and CD25 were explicit and consistent throughout the entire cohort. CD93 was present on BCR::ABL1+ LSCs in 7/8 patients and other previously suggested LSC markers like CD33 and CD56 specifically labeled BCR::ABL1+ primitive cells in only a fraction of the patients. The documented CML LSC marker IL1RAP was not in the ADT panel and hence could not be probed. In contrast, 11 ADTs were specifically detected in BCR::ABL1- primitive cells from at least one patient. Interestingly, CD35 with ~fourfold elevated expression as compared to BCR::ABL1+ LSCs was the only marker found to be consistently specific for the BCR::ABL1- HSC population across all patients. CD62L and CD10 were specifically present on BCR::ABL1- HSCs in 6/8 and 5/8 patients, respectively. Intriguingly, we recently showed that CD35 captures all HSCs in healthy human hematopoiesis and that stem cell activity as measured by xenotransplantation assays and epigenetic profiling is restricted to CD35+ HSCs (Sommarin et al., 2021). By applying the BCR::ABL1 signatures also to the Lin-CD34+ primitive cells (Figure 4—figure supplement 4), we could visualize the relative expression within the Lin-CD34+ UMAPs in relation to the cells BCR::ABL1 status and further confirm that the primitive fraction, previously captured by the established CML stem cell markers CD25, CD26, and CD93, consisted of specifically BCR::ABL1+ LSCs (Figure 5B–C, Figure 5—figure supplement 1 and Figure 5—figure supplement 2A). In addition, we observed that our previously described surface marker CD117, whose absence characterized LSCs persisting 3 months of TKI therapy, displayed a highly variably expression at diagnosis both marking a proportion of LSCs and HSCs across patients. Interestingly, despite that relative expression of CD35+ was lower in the primitive cells compared to MEP/ERP cell clusters the CD35+ gate captured ~50% of the BCR::ABL1- primitive cluster.
Given the consistency of detection across patients, we specifically focused on the ability of ADTs against CD25, CD26, and CD35 to separate between BCR::ABL1+ LSCs and BCR::ABL1- HSCs within the stem cell-enriched Lin-CD34+CD38-/low compartment. CD26 showed the highest capture of the BCR::ABL1+ primitive cluster (>95%) with a consistent capture percentage across patients (Figure 5D). This was followed by CD25 which captured >80% of the BCR::ABL1+ primitive cluster on average but displayed higher variability in capture efficiency compared to CD26 across patients. In comparison, the percentage of BCR::ABL1- primitive cluster captured by CD26 and CD25 was considerably lower (~5%). In contrast, the CD35 + gate showed a remarkable specific capture of BCR::ABL1- primitive cluster (~50%) (Figure 5D).
To further confirm whether CD26 and CD35 together could efficiently separate BCR::ABL1+ and BCR::ABL1- primitive cells, we gated for different combinations of CD26 and CD35 within the Lin-CD34+CD38-/low compartment (Figure 5E, Figure 5—figure supplement 2B). Predictably, whereas a combination of CD26+CD35- ADTs captured BCR::ABL1+ primitive LSC cluster at high purity, the CD26-CD35+ combination marked BCR::ABL1- primitive HSC cluster with no contamination of BCR::ABL1+ LSCs (Figure 5E). Interestingly, CD35 expression divided the CD26-, BCR::ABL1- primitive cells in one CD35 +and one CD35- fraction. This is in line with our recent observation that human BM immunophenotypic stem cells can be divided by CD35, where all HSC activity is captured within the CD35+ fraction (Sommarin et al., 2021). Accordingly, ADT combination of CD26-CD35- captured BCR::ABL1- primitive cells likely representing non-leukemic MPPs downstream of normal CD26-CD35+ stem cells. Notably, only a minority of cells were positive for ADTs for the CD26+CD35+ combination as evidenced by an overall low capture of cells in general. Thus, this shows the importance of positive, HSC-specific markers for efficient isolation of CML LSCs, and detection of residual HSCs in leukemia. Taken together, these observations suggested that although patients are heterogeneous for a given marker, BCR::ABL1+ LSCs can be consistently captured by the CD26+CD35- combination while BCR::ABL1- HSCs expressed a CD26-CD35+ immunophenotype within Lin-CD34+CD38-/low compartment in the bone marrow of CML patients.
A comparative analysis of CD26+CD35- versus CD26-CD35+ cells for BCR::ABL1 expression and response to TKI therapy
To explore whether the combination of CD26 and CD35 could indeed purify BCR::ABL1+ vs. BCR::ABL1- in Flow Cytometry protocols, we sorted CD26+CD35- and CD26-CD35+ cells from the Lin-CD34+CD38-/lowcompartment from patients both before and following 3 months of Bosutinib therapy and employed quantitative real-time qPCR analysis using Taqman probes against the BCR::ABL1 fusion gene (Figure 6A). The analysis showed a strong signal within CD26+CD35- cells for BCR::ABL1 expression while CD26-CD35+ cells were BCR::ABL1low/negative, validating the capacity of the combination of CD26 and CD35 to separate BCR::ABL1+ and BCR::ABL1- stem cells at diagnosis as well as after second generation TKI therapy. In addition, real-time qPCR analysis confirmed that Lin-CD34+CD38-/lowCD26+CD35- cells from CML patients at diagnosis displayed a strong signal for GAS2 (Figure 6B), one of the top genes from the CITE-seq BCR::ABL1+ LSC signature (Figure 4C). Furthermore, CD26+CD35- cells maintained GAS2 expression following 3 months of Bosutinib therapy (Figure 6B). GAS2 has previously been linked to CML disease progression, cell cycle control, apoptosis, and response to Imatinib (Janssen et al., 2005; Radich et al., 2006; Zhou et al., 2014).
To assess the level of quiescence between the BCR::ABL1+ LSCs and BCR::ABL1- HSCs, we sorted Lin-CD34+CD38-/lowCD45RA-CD26+CD35- and CD26-CD35+ single cells from two CML patients at diagnosis into individual wells and examined their division kinetics for 140 hr. However, LSC and HSC populations were indistinguishably deeply quiescent with average times to first division of 80–100 hr in culture (Figure 5—figure supplement 2D–C). Moreover, cell cycle analysis was performed by gating the corresponding cell populations in the Lin-CD34+CD38-/low CITE-seq data using ADTs. In accordance with the in vitro study, there was no distinct difference in cell cycle status between the populations (Figure 5—figure supplement 2E–F).
The higher burden of primitive cells in treatment failures observed above prompted us to query to which extent patients differed in the load of isolatable BCR::ABL1+ vs. BCR::ABL1- at diagnosis and their sustenance in response to TKI therapy. Intriguingly, we noted a substantially higher ratio of CD26+CD35- LSCs over CD26-CD35+ HSCs at diagnosis in prospective treatment failures compared to the bone marrow of optimal responders, both using CITE-seq ADTs as well as in FACS analysis for CD26 and CD35 (Figure 6C–D, Figure 5—figure supplement 2B). Importantly, following three months of TKI therapy, the percentage of BCR::ABL1+ CD26+CD35- cells decreased in both groups of patients; however, there was a striking difference in the relative proportion of CD26+CD35- LSCs/CD26-CD35+ HSCs in both groups. The stem cell compartment of optimal responders was dominated by healthy HSCs at three months following treatment. In contrast, the stem cell population of failure patients still comprised of CD26+CD35- LSCs without re-establishment of HSCs and restoration of normal hematopoiesis (Figure 6E). Taken together, we here for the first time present a protocol for separation and isolation of both BCR::ABL1+ LSCs and BCR::ABL1- HSCs from the same CML patients, two populations evidently relevant to therapy response.
Discussion
As new and more potent TKIs continue to be developed (Braun et al., 2020), the persistence of fully leukemogenic cells even in TFR necessitates life-long therapy in most patients. The ensuing treatment-emergent adverse events, and financial toxicity (Cortes et al., 2021; Lipton et al., 2022; Zafar, 2016) motivate a search for organizing principles to predict therapy response, dampen disease progression, and achieve durable TFR in a larger fraction of CML patients. One recent development has been to invoke the presence of somatic mutations other than BCR::ABL1 to explain primary resistance to TKIs, and disease progression (Branford et al., 2019). A complementary perspective is to consider leukemic cell state and heterogeneity as a fundamental determinant of response to TKI to augment the mutational foundation, as has been recently demonstrated in AML (Zeng et al., 2022).
Our single-cell multiomics maps show clear differences in overall cell composition within stem and progenitor compartments in leukemia patients at diagnosis versus nBM. Importantly, treatment failures and optimal responders displayed distinctive enrichment of specific cell clusters. Using our 11-molecularly defined clusters as anchors, we deconvoluted the bulk gene expression profiles from n=59 CML patients to infer constituent cell populations and found that a more profound stemness/primitive signature in the BM was consistently associated with inferior therapy response. Indeed, there is growing recognition of the burden of primitive cells, and the ratio of BCR::ABL1+ vs. BCR::ABL1- within the primitive compartment as potentially clinically relevant features such as hemoglobin, blast count, overall survival, progression-free survival, and therapy response (Fathy El-Metwaly et al., 2021; Krishnan et al., 2023; Mustjoki et al., 2013; Thielen et al., 2016). Apart from the differential burden, the primitive cells from optimal and sub-optimal responders could also be qualitatively different e.g., in terms of gene expression profile as suggested recently (Krishnan et al., 2023). The relevance of BCR::ABL1+ LSC burden to first and second generation TKI therapy outcome stands in contrast to the leukemic progenitors (Ph+ CD34+CD38+) which did not show such a correlation (Mustjoki et al., 2013). Moreover, longstanding observations have revealed that CML patients at diagnosis contain both BCR::ABL1+ LSCs vs. BCR::ABL1- HSC within the Lin-CD34+CD38-/low immunophenotypic compartment, and BCR::ABL1+ LSCs suppress residual normal HSCs (Chen et al., 2023; Coulombel et al., 1983). Predictably, a lower fraction of BCR::ABL1- normal stem cells at diagnosis and during therapy relates to hematological toxicity with delayed or compromised restoration of Phneg hematopoiesis upon successful response to TKIs (Janssen et al., 2012). Here, we for the first time present an effective means to separate CML LSCs from HSCs within the stem cell compartment of patients. Prospective optimal responders to TKI treatment had a higher content of BCR::ABL1-CD26-CD35+ cells at diagnosis, and 3 months following therapy versus the treatment failures. We have recently demonstrated that Lin-CD34+CD38-/lowCD45RA-CD35+ cells are at the top of human hematopoietic hierarchy and display chromatin accessibility and transcriptional enhancers in line with their multilineage long-term reconstitution ability (Sommarin et al., 2021).
With the high-throughput and resolution now provided by single-cell methods, BCR::ABL1+ cells can be distinguished from BCR::ABL1- cells residing within the same immunophenotypic compartment as described above, and importantly, heterogeneity within BCR::ABL1+ LSC subpopulations in terms of cell surface markers, molecular signature, and TKI response can be measured. The use of CITE-seq enabled the generation of patient-specific maps providing a panoramic yet exquisitely detailed view of cellular heterogeneity. Intriguingly, one of the emerging principles from the single-cell studies has been that not all LSC subpopulations are equally sensitive to BCR::ABL1 inhibition; while a majority of the LSCs are depleted, at least a fraction survives (Giustacchini et al., 2017; Warfvinge et al., 2017). This is corroborated by an overall reduction in BCR::ABL1+ Lin-CD34+CD38-/low cells but a striking enrichment of cKIT-CD26+ subpopulation displaying primitive and quiescent molecular program upon commencing TKI therapy (Warfvinge et al., 2017). Using CITE-seq, we confirmed the existence of these cells already at diagnosis and showed that BCR::ABL1+ vs. BCR::ABL1- can be efficiently discriminated as CD26+CD35- and CD26-CD35+, respectively within the Lin-CD34+CD38-/low fraction, and thus facilitating improved assessment of LSC burden.
The inspection of single-cell maps raises questions both regarding the origin of the interpatient heterogeneity in LSC burden as well as how these cells differentially manage to survive TKI therapy. Given the heterogeneity within normal HSCs with respect to lineage bias and clonal output (Haas et al., 2018), it is tempting to speculate that stem/progenitor cell acquiring the BCR::ABL1 oncogenic hit perhaps might be different across CML patients ensuring disparity in LSC load and characteristics. An equally plausible alternative is variability in the bone marrow niche that might favor differential abundance of stem/lineage-biased progenitors in patients (Baryawno et al., 2017). Once established, at least a subset of the LSCs is likely to be independent of BCR::ABL1 kinase signaling for their survival (Bhatia et al., 2003; Corbin et al., 2011; Hamilton et al., 2012). However, whether their self-sustenance involves switching to non-kinase activity of BCR::ABL1 protein, or achieving total BCR::ABL1 independence, and whether such mechanisms are inherently active in treatment naïve cells or become active only after exposure to TKIs is poorly understood (Zhao and Deininger, 2020).
We foresee future investigations to provide several key insights. First, accurate and direct detection of fusion transcripts such as BCR::ABL1 in single cells as part of scRNA-seq remains a major constraint, however, new approaches should allow direct measurement of BCR::ABL1 on a massive scale. Second, CITE-seq, a single-cell multiome approach, like any other high-throughput single-cell method comes with methodological, technical, and biological challenges. Cell sample handling, library preparation, sequencing quality and depth, drop-out of lowly expressed genes, and taking a molecular ‘snapshot’ of cells in time are all merely some factors that could introduce variability and noise to the data. A simpler, faster, and cheaper FACS panel detecting CD26+ and CD35+ surface markers can possibly work as an applicable clinical tool to quantify leukemic and non-leukemic stem cells in CML patients. Future prospective studies of this FACS panel coupled to clinical trials in large patient cohorts will establish the diagnostic and prognostic value of the relative abundance of these populations. It will be importance to evaluate if CD26-CD35+ cells are critical for restoring normal hematopoiesis once the TKI therapy diminishes the leukemic load, and whether patients with low counts of CD26-CD35+ cells at diagnosis have a relatively higher risk of developing hematologic toxicity such as cytopenia during therapy. The relative levels of CD26- and CD35+ stem cells for TFR will also be an important area of investigation.
Third, joint longitudinal analysis of leukemic as well as immune compartments is likely to be informative, especially for TFR. The observations that patients with activated immune signature at diagnosis are more likely to optimally respond to TKIs (Radich et al., 2023) serve as apt reminder for collective analyses. Fourth, our capacity to generate single-cell datasets, ironically outpaces our ability to extract information. Development of tools that allow biologists and clinicians to analyze large-scale datasets without requiring dedicated bioinformatics infrastructure can overcome the challenge (Dhapola et al., 2022). Finally, the clonal relationships and lineage output of candidate BCR::ABL1+ LSC subpopulations remain unknown. By coupling sc-multiomics with barcoding analysis and lineage tracing (Wagner and Klein, 2020), it should be possible to evaluate, for instance, whether the BCR::ABL1+CD26+ cells detected in patients in TFR are a sub-clone that persists from diagnosis through therapy, and whether it is responsible for recurrence in TFR.
The notion of LSCs has been invoked for a long time in the genetically defined and molecularly targeted paradigm of CML (Holyoake and Vetrie, 2017). Still, these entities remain essentially peripheral to the disease management (Zhao and Deininger, 2021). We posit that the heterogeneity of LSCs is a barrier towards their efficient measurement and safe purging. Our high-resolution single-cell multiomics maps suggest how to probe and deconstruct heterogeneity of CML, thereby permitting inference of leukemic vs. non-leukemic cells, estimation of BCR::ABL1+ LSCs, enumeration of their molecular features, and prospective isolation. Understanding how the cellular heterogeneity and plasticity emerges in the absence of extensive genetic variability would inform if the fully leukemogenic residual cells could either be safely eliminated pharmacologically or kept perpetually suppressed by an empowered immune system to avoid recurrence (Hsieh et al., 2021; Zhao and Deininger, 2023). Alternatively, longitudinal sampling and sc-omics combined with barcoding may also reveal if the stem cell subpopulations and states are interconvertible (Lenaerts et al., 2010; Marjanovic et al., 2013), and, therefore, appear as potentially inexhaustible pool that can only be slowly eroded by several years of TKI treatment in a subset of patients. A more formidable barrier is the unstated but still presumed equivalence of immunophenotype and function. Rather than describing LSCs solely by surface markers, we propose that treating leukemic cellular clones as the fundamental units of selection and evolution during therapy would have a more meaningful impact in predicting response to existing therapy and switch to another TKI. The next line of advances will require assessing LSCs as pool of clones defined by their ability to contribute to primary and secondary resistance in patients on therapy, and recurrence in TFR without recourse to animal models.
Materials and methods
Patient samples
Request a detailed protocolThe bone marrow was obtained from patients enrolled in either BosuPeg clinical trial (NCT03831776), NordCML006 (NCT00852566), BFORE (NCT02130557), or from Skåne University Hospital, Lund after informed consent and processed according to the guidelines approved by regional research ethics committees of sites involved in the trial. The information on age, gender, TKIs, BCR::ABL1IS %, and therapy response as per ELN guidelines are provided as Supplementary file 1. MNCs or CD34+ cells were enriched from the BM aspirates using magnetic microbeads and subsequently cryopreserved. The collection, processing, and generation of CITE-seq of bone marrow from age-matched healthy donors used in this study has been described previously (Sommarin et al., 2021).
CITE-seq sample preparation and facs sort
Request a detailed protocolSample preparation was performed according to Stoeckius et al., 2017 with minor adaptations. In brief, CML MNC, or CD34 enriched BM samples was thawed, FCR blocked 1/5 (130-046-703, Miltenyi) for 10 min, and washed before staining. Samples were stained with a PeCy5 conjugated lineage cocktail (CD2, CD3, CD14, CD16, CD19, CD235a), CD34 FITC, and a subset of the CITE-seq panel for 30 min on ice. Samples were washed and split into two tubes and stained with either CD38-PeCy7 or CD38-CITE-seq antibody followed by individual hashtags and the rest of the CITE-seq panel for 30 min on ice (Supplementary file 2 and Supplementary file 3). Cells were washed and resuspended in PBS, 2% FBS, 1/100 7AAD (Hyclone, 559925, BD Bioscience). Two populations (Lin-CD34+ and Lin-CD34+CD38-/low) were sorted per sample using FACSAriaII/Aria III (BD Bioscience) into 300 μl PBS, 0.04% BSA. Sorted cells were centrifuged at 300 × g for 10 min and volume was adjusted to 45 μl.
CITE-seq library preparation and sequencing
Request a detailed protocolCITE-seq library preparation and sequencing Cells were processed through Chromium Controller (10 x Genomics) with a total of 20–30.000 cells loaded per lane. Single-cell cDNA, HTO (Hashtags), and ADT (CITE-seq antibodies) libraries was prepared using Chromium Single-Cell 3’ V3 as per manufacturer instructions (CG000183 Rev C) as reported previously (Stoeckius et al., 2017). In brief, to increase HTO and ADT library yield, HTO and ADT-specific primers were spiked-in during the cDNA amplification step. cDNA and HTO/ADT libraries were isolated using 2 X SPRI beads per the manufacturers protocol. HTO and ADT libraries were subjected to adaptor ligation and sample index in separate PCRs using KAPA Hifi PCR Readymix (10 and 11 PCR cycles, respectively). The cDNA library was subjected to fragmentation, end repair, A-tailing, adaptor ligation, and sample index PCR (12 cycles). cDNA, HTO, and ADT libraries were sequenced together (65% cDNA, 30% ADT, 5% HTO) on Illumina sequencers with the following read length configuration: Read1=28, i7=8, i5=0, Read2=91. The raw data was processed using Cell Ranger 3.0.2 with GRCh38 as a reference genome.
CITE-seq analysis
Cell-gene matrices produced by Cell Ranger were analyzed using Scarf (v.0.18.12 Dhapola et al., 2022). In brief, cells with low or high gene count (<1000 or>9000 genes) and cells with high percentage mitochondrial gene UMIs (>11% or >6% for CML samples and healthy BM reference, respectively) were excluded from the analysis. The matrices were demultiplexed using Otsu thresholding (automatized pixel-based thresholding of HTO count per cell, background buffer = 0.1, override value = 0.5 Otsu thresholding was run twice for CML4, CML9. For nBM the thresholds were adjusted to HTO_1=32, HTO_2=18, HTO_3=17). HTO UMAPs were generated through Scarf (dims = 0) and clustered with Paris clustering (n_clusters = 2). If any RNA Leiden cluster overlapped with HTO paris cluster >70% and conversely, it was considered doublets and excluded from analysis.
Post-demultiplexing, for each sample, 2000 highly variable genes (HVGs) were identified (min_cells = 50) with the ‘mark_hvgs’ function in Scarf and used for PCA (principal component analysis) (CML8; HVGs = 1000, min_cells = 30). KNN (K-Nearest Neighbors) graphs were calculated based on the top 20 principal components (n_neighbours = 11) and utilized to build UMAP embedding of the cells (min_dist = 1, spread = 2, n_epochs = 2000). Clustering was performed using Lieden clustering (resolution set to 1 and 0.76 for CML and nBM samples, respectively).
Projection of Lin-CD34+ and Lin-CD34+CD38-/low CML subpopulations.
Mapping reference
Request a detailed protocolThe nBM is merged CITE-seq data of two Lin-CD34+ FACS sorted bone marrow samples (Sommarin et al., 2021) reanalyzed as described above’. Cluster identity was determined by (a) the top 10 genes in ‘run_ marker_search’ function in Scarf. In brief, each genes expression value is ranked per cell and a gene’s mean rank per cluster is divided by the sum of mean to determine the cluster specificity, (b) CLR (centered-log ratio) normed ADT expression across clusters was used to determine the immunophenotypic identity the cells, and (c) cell cycle scoring using the ‘run_cell_cycle_scoring’ function in Scarf was used to estimate cell cycle phase.
Label transfer through reference-based cell projection
Request a detailed protocolCML subpopulations were mapped to Lin-CD34+ nBM using Scarf (v.0.18.12). In short, cells were projected using ‘K-Nearest Neighbours (KNN) mapping through Scarf’s ‘run_mapping’ function (neighbors = 5). Cell cluster identity for CML cells was determined using the ‘get_target_classes’ function with the threshold set to 0.5 (>50% of the total weight score from the five top-matched reference cells must be cluster-specific to assign identity). Mapping scores (how frequently the reference cells end up as one of the top five neighbours of the projected CML cells) for the individual reference cells was calculated per projected sample using ‘get_mapping score.’.To visualize the mapping score and compare the mapping across samples, cells size of reference cells was set proportional to the score in the healthy bone marrow reference UMAP.
Following the mapping of Lin-CD34+CD38-/low CML cells (n=8), the primitive clusters from all patients were merged in Scarf. The ‘make_graph’ function in Scarf was used to identify the primitive CMLs cells top 500 HVGs (min_cells = 30, gender-specific genes excluded). Lin-CD34+CD38-/low healthy bone marrow cells (n=1) were projected onto the healthy reference as described above and the primitive cluster was added to the primitive CML cell clusters in Scarf. The UMAP was generated using the 500 CML-specific HVGs (k=11, PC = 20, n_epoch = 500) and with PCA fitted to the primitive CML cells only using Scarf’s ‘make graph’ function. Subsequently, cells were clustered using the Leiden clustering method (resolution = 0.1).
BCR::ABL1 status of the primitive clusters was determined using two DEG signatures: (1) Normal HSC vs BCR::ABL1+, and (2) BCR::ABL1- vs BCR::ABL1+ where Lin-CD34+CD38-/lowBCR::ABL1+ cells have been validated through single-cell RT-qPCR (Giustacchini et al., 2017). Up- and down-regulated genes were defined as log2 fold change >1 and <0, respectively. The expression of signature genes was normalized in Scarf, Z-scored to provide an average value count of the signature per cell and visualized in the UMAP using Scarf.
To visualize the BCR::ABL1+ primitive cells alongside the full CML heterogeneity, Lin-CD34+ and Lin-CD34+CD38-/low CML CITE-seq data from the same patient sample was merged in Scarf (The sub populations were FACS sorted at the time point and sequenced together). UMAPs of the merged data sets were generated as previously described (2000 HVGs, neighbors = 11, PC = 20). Cell indexes from BCR::ABL1+ and BCR::ABL1- primitive CML clusters was matched to the Lin-CD34+CD38-/low cell data.
ADT expression of BCR::ABL1+ and BCR::ABL1- cells were retrieved from the Zarr file in Scarf. The data thereafter was CLR (centered-log ratio) normalized and the log2 fold change was calculated between the two groups of cells using mean expression. p-values were calculated with the Mann–Whitney U test and they were corrected for multiple hypotheses testing using the Bonferroni method.
Analysis of differentially expressed genes
Request a detailed protocolDESeq2 was used to identify differentially expressed genes. Significantly up- or down-regulated was selected based on adjusted p-value <0.01 and log2 fold change >1 or < –1. Subsequently, volcano plots showing log2 fold change and -Log10 p-values were generated. Mitochondrial/ribosomal genes, as well as genes with a total count below 10 was excluded from the analysis.
To compare gene expression between clusters 1–11 from Lin-CD34+ CML cells, n=9 and their normal counterparts, their corresponding normal bone marrow reference clusters (n=2) were randomly pseudo-bulked into three replicates. Cluster-specific signatures genes were defined as significant up- or down-regulated genes unique to the cluster comparison, and not shared with not any other clusters. CML signature genes were defined as significantly up- or down- regulated genes present in all 11 clusters. Heatmap was generated from the DEseq2 data output (VST), showing the log2 mean expression of all samples per cluster. BCR::ABL1+ and BCR::ABL1- clusters >25 cells was used for DEG analysis of Lin-CD34+CD38-/low primitive cells (n=7).
ADT gating
Request a detailed protocolADT expression was CLR (centered-log ratio) normalized in Scarf. The normalized data was transformed into its antilog, multiplied with a scaling factor of 1000, and exported as FCS files. The FCS files was analyzed using FlowJo software. Gated cells ADT scale values were exported as CSV files and matched to cells ADT expression in CITE-seq data to determine original cell identity.
Gene expression deconvolution
Request a detailed protocolThe gene expression counts from CITE-seq molecularly defined populations served as input for CIBERSORTx (Newman et al., 2019) at default settings run on its docker version. Deconvolution of publicly available CD34+ bulk gene expression microarray dataset GSE14671 McWeeney et al., 2010 was performed after robust multichip average (RMA) normalization (Irizarry et al., 2003). Upon estimation of relative abundances of subpopulations from CD34+ cells, the values were normalized to 1, the value for each subpopulation represented its percentage/fraction within the query sample. Individual patients were annotated as non-responders or responders to Imatinib therapy as defined in the original study McWeeney et al., 2010; non-responders with >65% Ph + metaphases (the patients did not achieve even a minor cytogenetic response as described by the original authors), responders with 0% with Ph + metaphases (achieved CCyR) after 12 months of Imatinib therapy. Notably, the threshold for defining therapy response was proposed by the original authors, and the patients were labeled as either responder or non-responder (McWeeney et al., 2010). We adhered to the patient labels provided by the original authors as described in GSE14671.
Flow cytometry for phenotyping
Request a detailed protocolLymphoprep kits were used to separate mononuclear cells (MNCs), and magnetic microbeads were used to enrich CD34+ cells (Miltenyi). Antibodies against lineage-specific markers and those for specific populations mentioned in Supplementary file 3 were used to stain the cells. Sorting and analysis were done using FACS ARIAII/III or LSR FORTESSA (BD Biosciences); data analysis was done with FlowJo (Tree Star).
qPCR for BCR::ABL1
Request a detailed protocolRNA from FACS-sorted cells was purified using Single-Cell RNA Purification Kit (Norgen, cat# 51800) followed by cDNA synthesis using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, cat# 4368814). Following pre-amplification of cDNA (SsoAdvanced PreAmp Supermix BioRad, cat# 172–5160) for 12 cycles, qPCR was performed with TaqMan gene expression master mix (Applied Biosystems, cat# 4369016) and TaqMan probes for BCR::ABL1 (Assay ID: Hs03024541_ft), GAS2 (Assay ID: Hs00169477_m1 and GAPDH Assay ID: Hs02758991_ft).
Time to first division in vitro
Request a detailed protocolLin-CD34+CD38-/lowCD45RA-CD26+CD35- or CD26-CD35+ BM single cells from two CML patients were sorted into individual wells of 96-well U-shaped-bottom TPP plates using the FACS Aria III flow cytometer (Becton Dickinson). Cells were sorted directly into 200 μl Stem Span Serum-Free Expansion Medium (SFEM, StemCell Technologies), supplemented with 100 units/ml penicillin, 0.1 mg/ml streptomycin (Hyclone) and the following cytokines from PeproTech: SCF (100 ng/ml), Flt3L (100 ng/ml), TPO (50 ng/ml), and IL7 (10 ng/ml). Cells were visualized and counted in each well thrice a day using an inverted microscope for >140 hr. Dead cells as well as cells that did not divide during this period were excluded from the analysis. The Curve fit-sigmoid approach in GraphPad Prism was used to analyze the data.
Data availability
The data from sequencing have been deposited in GEO under accession ID: GSE236233. Previously published data used within this study is available at: GSE173076 (nBM), and GSE14671 (deconvolution). Data and code are available under the following OSF repository: https://osf.io/ns2j8/. The CITE-seq data has been published in an analyzed, ineractive format under the following CellHub repository: https://scarfweb.nygen.io/cellhub/2b2fec2ceb52.
-
NCBI Gene Expression OmnibusID GSE236233. Single cell multi-omics analysis of chronic myeloid leukemia links cellular heterogeneity to therapy response.
-
NCBI Gene Expression OmnibusID GSE173076. Single-Cell Multiomics Reveals Distinct Cell States at the Top of the Human Hematopoietic Hierarchy [CITE-seq].
-
NCBI Gene Expression OmnibusID GSE14671. Expression signature to predict major cytogenetic response in chronic phase CML patients treated with imatinib.
References
-
Polyclonal hematopoiesis in interferon-induced cytogenetic remissions of chronic myelogenous leukemiaBlood 79:997–1002.
-
Human chronic myeloid leukemia stem cells are insensitive to imatinib despite inhibition of BCR-ABL activityThe Journal of Clinical Investigation 121:396–409.https://doi.org/10.1172/JCI35721
-
Cytogenetic studies in patients on imatinibSeminars in Hematology 40:50–55.https://doi.org/10.1053/shem.2003.50043
-
CD34+/CD38- stem cell burden could predict chronic myeloid leukemia patients’ outcomeAsian Pacific Journal of Cancer Prevention 22:3237–3243.https://doi.org/10.31557/APJCP.2021.22.10.3237
-
Titrating complex mass cytometry panelsCytometry. Part A 95:792–796.https://doi.org/10.1002/cyto.a.23751
-
Histone lysine methyltransferases in biology and diseaseNature Structural & Molecular Biology 26:880–889.https://doi.org/10.1038/s41594-019-0298-7
-
BRD7, a novel PBAF-specific SWI/SNF subunit, is required for target gene activation and repression in embryonic stem cellsThe Journal of Biological Chemistry 283:32254–32263.https://doi.org/10.1074/jbc.M806061200
-
Cell plasticity and heterogeneity in cancerClinical Chemistry 59:168–179.https://doi.org/10.1373/clinchem.2012.184655
-
Characterization of primitive subpopulations of normal and leukemic cells present in the blood of patients with newly diagnosed as well as established chronic myeloid leukemiaBlood 88:2162–2171.
-
Identification of CD25 as STAT5-dependent growth regulator of leukemic stem cells in Ph+ CMLClinical Cancer Research 22:2051–2061.https://doi.org/10.1158/1078-0432.CCR-15-0767
-
Evidence of ABL-kinase domain mutations in highly purified primitive stem cell populations of patients with chronic myelogenous leukemiaBiochemical and Biophysical Research Communications 323:728–730.https://doi.org/10.1016/j.bbrc.2004.08.169
-
Simultaneous epitope and transcriptome measurement in single cellsNature Methods 14:865–868.https://doi.org/10.1038/nmeth.4380
-
Mechanisms of non-metastatic 2 (NME2)-mediated control of metastasis across tumor typesNaunyn-Schmiedeberg’s Archives of Pharmacology 384:397–406.https://doi.org/10.1007/s00210-011-0631-0
-
Human haematopoietic stem cell lineage commitment is a continuous processNature Cell Biology 19:271–281.https://doi.org/10.1038/ncb3493
-
Lineage tracing meets single-cell omics: opportunities and challengesNature Reviews. Genetics 21:410–427.https://doi.org/10.1038/s41576-020-0223-2
-
Dpy30 is critical for maintaining the identity and function of adult hematopoietic stem cellsThe Journal of Experimental Medicine 213:2349–2364.https://doi.org/10.1084/jem.20160185
-
Financial toxicity of cancer care: It’s time to interveneJournal of the National Cancer Institute 108:djv370.https://doi.org/10.1093/jnci/djv370
-
Declaration of Bcr-Abl1 independenceLeukemia 34:2827–2836.https://doi.org/10.1038/s41375-020-01037-9
-
Eradicating residual chronic myeloid leukaemia: basic research lost in translationThe Lancet. Haematology 8:e101–e104.https://doi.org/10.1016/S2352-3026(21)00001-6
Article and author information
Author details
Funding
Swedish Cancer Foundation
- Göran Karlsson
Ragnar Söderbergs stiftelse
- Göran Karlsson
Knut och Alice Wallenbergs Stiftelse
- Göran Karlsson
Vetenskapsrådet
- Göran Karlsson
Barncancerfonden
- Göran Karlsson
Incyte Biosciences Nordic AB
- Rebecca Warfvinge
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to acknowledge patients, study nurses, and other personnel in the clinical centers for their participation in this project. We also thank Clinical Genomics Lund, SciLifeLab and Center for Translational Genomics (CTG), Lund University, for providing expertise and service with sequencing and analysis, and the Lund Stem Cell Center FACS Facility for expert Flow Cytometry technical support. This work was funded by grants from the Swedish Cancer Society, the Ragnar Söderberg Foundation, the Knut and Alice Wallenberg Foundation, the Swedish Research Council, the Swedish Childhood Cancer fund, and a grant from Incyte Biosciences Nordic AB.
Ethics
Clinical trial registration NCT03831776, NCT00852566, NCT02130557.
The bone marrow was obtained from patients enrolled in either BosuPeg clinical trial (NCT03831776), NordCML006 (NCT00852566), BFORE (NCT02130557), or from Sk;ne University Hospital, Lund after informed consent and processed according to the guidelines approved by regional research ethics committees of sites involved in the trial.
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
- Version of Record published:
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.92074. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, Warfvinge et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,461
- views
-
- 116
- downloads
-
- 2
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Cancer Biology
- Genetics and Genomics
The expression of eukaryotic genes relies on the precise 3'-terminal cleavage and polyadenylation of newly synthesized pre-mRNA transcripts. Defects in these processes have been associated with various diseases, including cancer. While cancer-focused sequencing studies have identified numerous driver mutations in protein-coding sequences, noncoding drivers – particularly those affecting the cis-elements required for pre-mRNA cleavage and polyadenylation – have received less attention. Here, we systematically analysed somatic mutations affecting 3'UTR polyadenylation signals in human cancers using the Pan-Cancer Analysis of Whole Genomes (PCAWG) dataset. We found a striking enrichment of cancer-specific somatic mutations that disrupt strong and evolutionarily conserved cleavage and polyadenylation signals within tumour suppressor genes. Further bioinformatics and experimental analyses conducted as a part of our study suggest that these mutations have a profound capacity to downregulate the expression of tumour suppressor genes. Thus, this work uncovers a novel class of noncoding somatic mutations with significant potential to drive cancer progression.
-
- Cancer Biology
- Immunology and Inflammation
The immunosuppressive microenvironment in pancreatic ductal adenocarcinoma (PDAC) prevents tumor control and strategies to restore anti-cancer immunity (i.e. by increasing CD8 T-cell activity) have had limited success. Here, we demonstrate how inducing localized physical damage using ionizing radiation (IR) unmasks the benefit of immunotherapy by increasing tissue-resident natural killer (trNK) cells that support CD8 T activity. Our data confirms that targeting mouse orthotopic PDAC tumors with IR together with CCR5 inhibition and PD1 blockade reduces E-cadherin positive tumor cells by recruiting a hypoactive NKG2D-ve NK population, phenotypically reminiscent of trNK cells, that supports CD8 T-cell involvement. We show an equivalent population in human single-cell RNA sequencing (scRNA-seq) PDAC cohorts that represents immunomodulatory trNK cells that could similarly support CD8 T-cell levels in a cDC1-dependent manner. Importantly, a trNK signature associates with survival in PDAC and other solid malignancies revealing a potential beneficial role for trNK in improving adaptive anti-tumor responses and supporting CCR5 inhibitor (CCR5i)/αPD1 and IR-induced damage as a novel therapeutic approach.