A 4-gene leukemic stem cell score can independently predict the prognosis of myelodysplastic syndrome patients

Myelodysplastic syndromes are clonal myeloid malignancies arising from hematopoietic stem cells (HSCs), which are characterized by bone marrow (BM) ineffective hematopoiesis, peripheral blood cytopenias,^1-5  and a propensity of transformation to acute myeloid leukemia (AML).^1-3  The clinical features and outcomes of MDS patients vary considerably, underscoring the importance of individualized management. The International Prognostic Scoring System (IPSS)⁶  and the revised IPSS (IPSS-R)⁷  have been widely used to risk-stratify MDS patients and guide the choice of treatment. However, the prognosis of patients may be different even in the same risk groups.^8,9  The search for new prognostic markers is needed for better risk classification of MDS patients.

Recently, Ng et al proposed a 17-gene leukemic stem cell (LSC) scoring system (LSC17 score) that could accurately predict the outcomes of AML patients.¹⁰  As the clonal origin of MDS and AML has been both demonstrated to lie within the stem cell compartment and tumor stemness is an established property pertaining to poor prognosis,^11-14  we surmise that LSC gene expression may harbor clinical significance in MDS patients. By using Z transformation and secondary multivariate Cox regression analysis, we identified expressions of 4 LSC genes were significantly correlated with overall survival (OS). We constructed a concise, integrated LSC signature–based 4-gene scoring system (LSC4) and found the LSC4 score was closely associated with clinical and biological features and could predict OS and leukemia-free survival (LFS) in MDS patients. The prognostic implication of the LSC4 score on OS remained significant in both IPSS-R lower-risk (very low, low, and intermediate risk) and higher-risk (high- and very-high-risk) groups. We also validated the prognostic significance of the LSC4 score in an independent internal cohort and 1 external cohort from GSE58831 in which microarray data were available.¹⁵  Furthermore, it was an independent risk factor for OS and LFS, irrespective of age, IPSS-R, and mutation status. We believe that after prospective validation with independent cohorts, this scoring system might provide directives to therapeutic decision in MDS patients, especially lower-risk-group patients in whom usually only supportive treatment is given. More aggressive therapy may be indicated for IPSS-R lower-risk patients with higher LSC4 scores.

We recruited 176 primary MDS patients diagnosed at the National Taiwan University Hospital (NTUH) from January 1992 to December 2010 who had cryopreserved BM samples for microarray analysis as a training cohort. The diagnosis was based on the 2016 World Health Organization (WHO) classification.¹⁶  Patients with antecedent chemotherapy or hematologic malignancies were excluded. Another independent set of 30 patients diagnosed with the same criteria from January 2011 to May 2012 were recruited as an internal validation cohort. For external validation, we collected publically available data from GSE58831¹⁵  in which 176 patients were annotated, but survival and microarray data for gene expression levels were both available in only 113 MDS patients. This study was approved by the institutional review board of the NTUH.

The median age of the 176 MDS patients in the training cohort was 68.7 years. Among the 164 patients who had karyotype data at diagnosis, 17.1% had IPSS-R very-high-risk MDS, 21.3% were high risk, 25% were intermediate risk, 32.9% were low risk, and 3.7% were very low risk (Table 1). Most patients (121 patients [68.8%]) received supportive care (only because a hypomethylation agent [HMA] was not reimbursed by National Health Insurance in Taiwan until 2013), and 55 received active treatment, including HMA, low-dose cytarabine, high-intensity chemotherapy, and/or HSC transplantation (HSCT). Nineteen (10.8%) patients, including 15 with MDS with excess blasts (EBs) and 4 with progressive disease during follow-up, underwent allogeneic HSCT. Approximately 10% patients in each group received supportive care due to the patient’s choice and/or comorbidity. During a median follow-up duration of 37.3 months (range, 0.1-130.9 months), 83 patients died of the disease, and 38 progressed to AML.

We profiled the global gene expression of BM mononuclear cells from the 206 patients by Affymetrix GeneChip Human Transcriptome Array 2.0 as described previously.¹²  The raw and normalized microarray data reported in this article have been deposited in the Gene Expression Omnibus with the accession number GSE97064. The expression levels of the previously identified 17 LSC genes,¹⁰  including MMRN1, DPYSL3, CDK6, LAPTM4B, NGFRAP1, CD34, AKR1C3, EMP1, SOCS2, NYNRIN, KIAA0125, GPR56, SMIM24, DNMT3B, CPXM1, ZBTB46, and ARHGAP22, were extracted for further analysis. Cytogenetic analyses were performed as described previously¹⁷  and interpreted according to the International System for Human Cytogenetic Nomenclature.¹⁸  We also analyzed the mutation statuses of 17 myeloid-relevant genes, including ASXL1, IDH1, IDH2, EZH2, TET2, and DNMT3A, genes related to the RNA-splicing machinery (including SF3B1, U2AF1, SRSF2, and ZRSR2), as well as FLT3/internal tandem duplication, NRAS, KRAS, RUNX1, MLL/partial tandem duplication, TP53, and SETBP1, by Sanger sequencing as previously described.^17,19-24 

We first conducted Z transformation for the expression of the 17 LSC genes at probe levels across the 176 MDS patients and set 0 as the mean and calculating unit standard deviation of each gene among the patients. We then used a multivariate approach to analyze the association between OS and the expression level of each LSC gene.

The LSC genes with significant association with OS were assigned for further multivariate Cox regression analysis to find genes independently associated with prognosis. These LSC genes were used to build the LSC prognostic scoring system. We performed the secondary multivariate Cox regression analysis to obtain the β values as the LSC genes’ weights in the scoring system. The prognostic LSC score was calculated as the sum of the normalized expression level of each component multiplied by its weight as follows: risk(j) = ΣLSC component LSC_i (j) × β_i, where j denotes the patient accession number, LSCi represents the normalized expression level of the LSC probe i after Z transformation, and β_i is the weight of the particular LSC probe i.

We used the Mann-Whitney U test to compare medians and continuous variables of distribution. The Fisher exact test or the χ² test was performed to examine the difference among discrete variables, including gender, WHO classification, cytogenetic changes, IPSS-R, and genetic alterations between patients with lower and higher LSC scores. OS was the duration from the date of initial diagnosis to the time of last follow-up or death of any cause, whichever occurred first. The training set was used to build the LSC scoring system, which was then applied to the validation set to confirm its significance. The survival prediction power of this LSC score was evaluated by both the log-rank test and the univariate Cox proportional hazards model. After >100 000 iterations, the prediction rate of our proposed LSC score was calculated as the fraction of random scoring systems that achieved P < .05. We plotted the survival curves with Kaplan-Meier analysis and calculated the statistical significance with the log-rank test. The Cox proportional hazards model was used in multivariate regression analysis. P < .05 was considered statistically significant. All statistical analyses were performed with BRB-ArrayTools (version 4.5.1; Biometric Research Branch, National Cancer Institute, Rockville, MD) and IBM SPSS Statistics 23 for Windows. Time-dependent receiver operating characteristic (ROC) curves analysis was performed using the R package time ROC. The Pearson’s correlation coefficient (PCC) was calculated in R language.

We first applied the LSC17 score constructed by Ng et al¹⁰  to 176 MDS patients in our training cohort and divided patients into higher- and lower-score groups with the median value as a cutoff level. We noticed that the LSC17 score could truly stratify the total MDS cohort into 2 risk groups with different OS and LFS, but further subgroup analyses showed that the scoring system had no prognostic significance in either lower- or higher-risk MDS subgroup according to IPSS-R or WHO classification (supplemental Table 1). Therefore, the prognostic-predicting capability of the LSC17 score is not as good in MDS as in AML, suggesting the more heterogeneous nature of MDS.

To construct a more simplified and powerful prognostic scoring system based on relevant LSC signature, we put the 17 LSC-related genes in a multivariate Cox model to identify the genes whose expression could independently predict OS (supplemental Table 2). We found the expression levels of LAPTM4B, NGFRAP1, CPXM1, CDK6, NYNRIN, and EMP1 genes were correlated with survival (P = 0.017, 0.03, 0.042, 0.05, 0.072, and 0.078, respectively). We then performed another round of Cox regression analysis for these 6 genes. The expression levels of LAPTM4B, NGFRAP1, EMP1, and CPXM1 remained significantly correlated with survival (P = .001, .027, .02, and .001, respectively). By integrating the β values as statistical weights, we constructed the LSC4 score, which was calculated with the following equation: [LAPTM4B] × 0.502 − [NGFRAP1] × 1.013 + [EMP1] X 0.181 + [CPXM1] × 0.381 (supplemental Table 3).

The LSC4 scores of the 176 MDS patients were calculated, and the relationship between LSC4 score and the expression of LSC-related genes were investigated (Figure 1; supplemental Figure 1). The LSC4 score had a strong correlation with CPXM1 and LAPTM4B expression (PCC = 0.81 and 0.78, respectively), as well as a moderate correlation with EMP1 (PCC = 0.55), but no correlation with NGFRAP1 (PCC = 0.02), consistent with coefficients of the LSC4 equation. The large range of expression correlations among the 4 selected genes (PCC = −0.02 to 0.60) suggests the complement of these genes for the prognostic prediction. Although some genes are highly coexpressed with CPXM1, such as NYNRIN (PCC = 0.79), CD34 (PCC = 0.68), DNMT3B (PCC = 0.60), and KIAA0125 (PCC = 0.48), only CPMX1 was selected as one of the predominant markers. It demonstrates that our procedure for marker selection can exclude redundant markers.

The 176 MDS patients in the training set were divided into 2 groups by the median value of the LSC4 scores. A comparison of clinical and laboratory features between the 2 groups is shown in Table 1. The high-score group had higher BM blast percentages at diagnosis (P < .001) compared with the low-score group. Patients with higher scores more frequently had MDS-EB by the 2016 WHO classification, including EB-1 and EB-2, but less MDS-SLD, MDS-MLD, MDS-RS-SLD, and MDS-RS-MLD compared with patients with lower scores. High-score patients were more frequently categorized into IPSS-R high- and very high-risk subgroups but less frequently to the IPSS-R low- and very-low-risk subgroups (Table 1).

Moreover, high-score patients had significantly higher incidence of poor-risk cytogenetics (21.2% vs 5.1%, P = .003) and complex karyotypes (≥3 abnormalities, 17.6% vs 3.8%, P = .005) (supplemental Table 4). Overall, 108 patients (61.4%) had at ≥1 gene mutation, 36 (40.9%) in the low-score subgroup and 72 (81.8%) in the high-score subgroup (P < .001). As listed in supplemental Table 5, the most common mutation in the high-score patients was ASXL1 mutation (36%), followed by RUNX1 (26.7%), SRSF2 (21.8%), and DNMT3A (18.4%) mutations. In contrast, SF3B1 mutation was the most frequent mutation (23.9%) in the low-score patients. High LSC4 score was closely associated with ASXL1, RUNX1, SRSF2, TP53, U2AF1, and ZRSR2 mutations but inversely associated with SF3B1 mutation (supplemental Table 5).

Patients with higher LSC4 scores had an inferior OS and LFS than those with lower scores (median, 14.6 months vs 83.6 months, P < .001; and 10.3 months vs 83.6 months, P < .001, respectively; Figure 2A-B). Subgroups analysis showed that the prognostic significance of LSC4 score for OS and LFS remained true in both IPSS-R lower-risk (very low, low, and intermediate risk) subgroup (Figure 3A-B) and IPSS-R higher-risk (high and very high risk) subgroup (Figure 3C-D). To further compare the power between LSC4 score and IPSS-R in prognostic prediction, the time-dependent ROC curves were estimated by the inverse probability of censoring weighting method. As illustrated in Figure 4, LSC4 score could be complementary to IPSS-R in prediction of both OS and LFS. Additionally, we performed bivariate analyses, including the LSC4 score and the IPSS-R, and LSC4 score and the presence of a complex karyotype, respectively. As a result, there was a moderate correlation between LSC4 score and IPSS-R (PCC = 0.536, P < .001) and a low correlation between LSC4 score and the presence of a complex karyotype (PCC = 0.221, P = .004).

The adverse implications of higher LSC4 scores on OS and LFS could also be demonstrated in the subgroups of patients with normal karyotype (n = 100; supplemental Figure 2A-B); patients without unfavorable cytogenetics such as complex karyotypes, monosomy 7, and del(7q) (n = 142, supplemental Figure 2C-D); and those with WHO lower-risk subtypes (MDS-SLD, MDS-MLD, MDS-RS-SLD, and MDS-RS-MLD; n = 100; supplemental Figure 3A-B) or WHO higher-risk subtypes (MDS-EB; n = 76) (supplemental Figure 3C-D).

We further analyzed separately the outcomes of MDS patients receiving different treatments. High-score patients consistently had a significantly inferior outcome in OS and LFS (supplemental Figure 4), regardless of whether they received palliative care (n = 121, 15.4 months vs 83.6 months, P < .001; and 13.1 months vs 83.6 months, P < .001, respectively) or active treatment (n = 55, 14.6 months vs 42.4 months, P < .001; 9 months vs 50.4 months, P < .001, respectively). For the patients treated with HMA (n = 17), those with a high LSC4 score also showed a significantly shorter OS and LFS (supplemental Figure 5; 14.6 months vs NR, P = .047, and 7 months vs not reached [NR], P = .02, respectively).

For multivariate analysis in the training cohort, we included parameters with P < .05 in univariate Cox regression analysis as covariates, including age, and mutations in ASXL1, TP53, SRSF2, and ZRSR2 (supplemental Table 6). Higher LSC4 score, either divided by a median (Table 2) or regarded as continuous values (supplemental Table 7), appeared to be an independent adverse prognostic factor for OS (P < .001 and P < .001, respectively) and LFS (P < .001 and P < .001, respectively). To verify the prognostication power of the LSC4 scoring system, we analyzed the expression levels of 17 LSC genes in an independent internal validation cohort of 30 MDS patients. Characteristics of patients in the training cohort and internal validation cohort were generally comparable, as shown in supplemental Table 8. Consistent with the findings in the training cohort, patients with higher LSC4 scores had a shorter OS and LFS (Figure 5A-B; 6.9 months vs NR, P = .005, and 3.9 months vs NR, P = .002, respectively) than those with lower scores in the validation cohort. Moreover, we validated the prognostic significance of this LSC4 scoring system in one external validation cohort of 113 MDS patients from GSE58831.¹⁵  A comparison of clinical and laboratory features between high- and low-score patients in this cohort is shown in supplemental Table 9. The high-score group had higher BM blast percentages at diagnosis (P = .004) and more frequently IPSS intermediate-2 subtype (P = .003) but lower absolute neutrophil counts (P = .005) and platelet counts (P = .029) compared with the low-score group. Low-score patients had a higher incidence of SF3B1 mutation than high-score patients, but there was no difference in cytogenetic changes and other molecular gene mutations between the 2 groups (supplemental Tables 10 and 11). Patients with higher LSC4 scores had significantly shorter OS (Figure 5C; 37.3 months vs NR, P = .01) than those with lower scores.

MDS is a heterogeneous disease with highly variable clinical presentations, survival, and rate of leukemia transformation among the patients. Traditional risk stratification systems like IPSS or IPSS-R have been widely adopted in the care of MDS patients for a long time. Compared with IPSS, the IPSS-R incorporated more parameters, including blast percentages, comprehensive cytogenetic classification, and the depth of cytopenias, with improved prediction power for prognosis assessment in MDS patients.⁷  However, the clinical outcomes of MDS patients may vary even in the same IPSS-R risk groups. It is warranted to search for more prognostic markers for better risk stratification.

In this study, we tried to integrate expression of genes related to stemness into the prognostic model. In a recent report, Ng et al proposed a LSC17 scoring system that clearly predicted the outcomes of AML patients.¹⁰  Since stemness of leukemia cells is a crucial factor of resistance to chemotherapy, it is reasonable that the expression profile of these LSC genes has significant impact on prognosis of AML patients.^{5,10,13,14,25,26}  MDS and AML are both myeloid malignancies rising from stem/progenitor cells, and they share some clinical and biological features.^11,27-31  The aberrant stem and progenitor cell populations have some cellular features in common with normal HSCs, including sustained self-renewal and proliferation capacity.²⁸  They are responsible for disease initiation, transformation, and relapse and are also more resistant to chemotherapy.²⁸  Multiple reports have shown functionally defined subsets of AML stem cells using in vitro and xenograft model systems.^14,32-34  In contrast, while MDS has also long been considered as a stem cell disorder, the characterization of LSCs in MDS is still not clear yet. We hypothesize that expression levels of LSC genes have a significant impact on the prognosis of MDS patients, although the roles of these LSC gene signatures in prognostication in MDS patients have not been clarified yet. We first applied the LSC17 score proposed by Ng et al¹⁰  to our MDS patient cohort but found that its correlation with clinical outcome was unsatisfactory. This might somewhat reflect the different nature of LSCs between MDS and AML.

In this study, we constructed a simple 4-LSC gene signature–based score for prediction of clinical outcomes, considering the expression levels and weights of only 4 LSC genes. We demonstrated that higher LSC4 scores predicted poorer prognosis for both OS and LFS. Patients with higher LSC4 score had worse OS and LFS than others in the same IPSS-R or 2016 WHO subgroup. These findings indicate the clinical heterogeneity in the same risk groups. Patients in each IPSS-R risk group could be further stratified into different prognostic groups based on this LSC4 scoring system. In addition, patients with higher LSC4 score consistently had poorer OS and LFS across subgroups, including patients receiving supportive care (P < .001 for both OS and LFS; supplemental Figure 4) or active treatment (P < .001 for both OS and LFS). These analyses indicated the prognostic power and clinical relevancy of the score, since more proactive follow-up and treatment strategies can be considered if a patient harbors a high LSC4 score, regardless of which risk group he or she was initially categorized in.

The prognostic significance of the LSC4 scoring system was confirmed in 2 independent validation cohorts. In bivariate analysis, there was a moderate correlation between LSC4 scores and IPSS-R, which is predictable, since both prognostic systems could well stratify patients’ outcome, yet the underlying mechanism warrants further studies. Meanwhile, a low correlation between LSC4 score and the presence of a complex karyotype further underlines the independent prognostic significance of LSC4 score. Although higher LSC scores were closely associated with high-risk mutations such as RUNX1, ASXL1, SRSF2, and TP53 mutations, multivariate analysis proved high LSC score to be an unfavorable prognostic factor independent of the IPSS-R and mutations.

To our knowledge, this is the first study integrating LSC gene signatures to predict the prognosis in MDS patients. This 4-LSC gene signature–based score is concise yet powerful and easier for clinical applications than using 17 genes. This LSC4 scoring system would be helpful for identifying those with poorer prognosis in the rather heterogeneous group of patients, especially those with IPSS-R lower-risk MDS; in patients with IPSS-R lower-risk but high LSC4 scores, more aggressive treatment, rather than traditionally palliative care, might be warranted.

From the 17 LSC genes, we identified that expression of LAPTM4B, NGFRAP1, EMP1, and CPXM1 is most relevant to the prognosis of MDS patients. LAPTM4B (the lysosome-associated protein transmembrane-4β gene) is considered as an oncogene, and its overexpression has been proven to promote the proliferation of various tumor cells, boost invasion and metastasis, resist apoptosis, initiate autophagy, and assist drug resistance.^35-38  Overexpression of LAPTM4B in breast cancer cells results in resistance to anthracycline, and its knockdown can sensitize the drug response.³⁹ NGFRAP1 (nerve growth factor receptor-associated protein 1) is an apoptosis-related gene, and its expression is downregulated in some solid organ malignancies and chronic lymphocytic leukemia.^40,41  Overexpression of NGFRAP1 was reported to be correlated with an adverse outcome in AML patients.⁴² EMP1 (epithelial membrane protein 1) encodes proteins associated with membrane blebbing, cell proliferation, and squamous cell differentiation.^43-45  Its high expression was an independent predictor of poor outcome in pediatric leukemia,⁴⁶  while its roles in solid cancers are still not fully understood.^47-50 CPXM1 (carboxypeptidase X, M14 family member 1) encodes a member of the carboxypeptidase family of proteins and is a positive regulator of adipogenesis, which may contribute to hyperplastic adipose tissue expansion via affecting extracellular matrix remodeling.⁵¹  Recently, a gene expression signature was developed to predict early molecular remission and long-term outcome in chronic myeloid leukemia patients and higher expression of CPXM1 in this signature was positively correlated with early molecular remission failure rate in chronic-phase chronic myeloid leukemia patients on frontline imatinib.⁵²  However, the role of CPXM1 in hematological and oncological malignancies is still largely unknown.⁵³  One interesting thing is that the directions of NGFRAP1 and CPXM1 in this LSC4 scoring system were different from those in the LSC17 equation for AML, implicating the different nature of LSCs between MDS and AML. Further studies of the underlying pathophysiology and mechanisms are worth exploring in the future. Overall, although these genes are regarded as LSC genes, their roles in MDS or other hematologic malignancies are still not fully understood, and further studies are needed to explore the functional pathways and pathogenesis of these genes in MDS.

There are limitations of this study. Firstly, we used array-based approaches rather than next-generation sequencing (NGS)–based methods to quantify gene expression. The advantages of NGS include higher sensitivity for genes with low expression and the broader dynamic range of expression levels in comparison with microarray. However, several reports have revealed that there is a strong concordance between microarray and RNA-sequencing data, and they have similar performance for the prediction of clinical end points.^54,55  Therefore, gene expression profiles generated by array-based methods are still valuable resources for biomarker discovery. Secondly, unlike mutations, which are clear-cut, standardizing the quantification of gene expression is a daunting task when it comes to applicability of this scoring system in clinical practice. To apply the scoring system in a clinical setting like BCR-ABL1 levels in chronic myelogenous leukemia, there are still many problems to solve. Using novel technologies such as the NanoString nCounter system or targeted RNA sequencing, the LSC signature–based scoring system can be applied in clinical practice more easily. In conclusion, this is the first study integrating the LSC gene signature into risk stratification of MDS patients. This study provides a novel, powerful, and biologically significant scoring system, which serves as an independent prognostic factor for OS and LFS in MDS patients. This integrated prognostic system refines the prognostic prediction models and might guide the therapeutic decision and possible LSC-targeted therapy in the future. Lastly, although the reliability and reproducibility of this prognostic model had been validated by 2 independent cohorts, prospective studies with more standardization and more precise quantification methods, such as NGS or real-time polymerase chain reaction, are warranted.

The authors would like to acknowledge the service provided by Department of Laboratory Medicine, Department of Medical Research, and Division of Hematology, Department of Internal Medicine, NTUH.

This study was supported by grants from the Ministry of Health and Welfare, Taiwan (project number MOHW107-TDU-B-211-114009) and the Ministry of Science and Technology, Taiwan (project numbers MOST 107-2314-B-002-013 and MOST 108-2314-B-002-011).

Contribution: Y.-H.W. and C.-C.L. were responsible for data collection and management, statistical analysis and interpretation, literature research, and manuscript writing; C.-Y.Y. was responsible for data management and statistical analysis; C.-L.H. assisted in statistical analysis; H.-A.H. and C.-H.T. were responsible for data collection and management; and W.-C.C. and H.-F.T. planned, designed, and coordinated the study over the entire period and wrote the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Wen-Chien Chou, Department of Laboratory Medicine, National Taiwan University Hospital, No. 7, Chung-Shan S. Rd, Taipei City 10002, Taiwan; e-mail: wchou@ntu.edu.tw; and Hwei-Fang Tien, Department of Internal Medicine, National Taiwan University Hospital, No. 7, Chung-Shan S. Rd, Taipei City 10002, Taiwan; e-mail: hftien@ntu.edu.tw.