A multiple breast cancer stem cell model to predict recurrence of T1–3, N0 breast cancer

Background Local or distant relapse is the key event for the overall survival of early-stage breast cancer after initial surgery. A small subset of breast cancer cells, which share similar properties with normal stem cells, has been proven to resist to clinical therapy contributing to recurrence. Methods In this study, we aimed to develop a prognostic model to predict recurrence based on the prevalence of breast cancer stem cells (BCSCs) in breast cancer. Immunohistochemistry and dual-immunohistochemistry were performed to quantify the stem cells of the breast cancer patients. The performance of Cox proportional hazard regression model was assessed using the holdout methods, where the dataset was randomly split into two exclusive sets (70% training and 30% testing sets). Additionally, we performed bootstrapping to overcome a possible biased error estimate and obtain confidence intervals (CI). Results Four groups of BCSCs (ALDH1A3, CD44+/CD24−, integrin alpha 6 (ITGA6), and protein C receptor (PROCR)) were identified as associated with relapse-free survival (RFS). The correlated biomarkers were integrated as a prognostic panel to calculate a relapse risk score (RRS) and to classify the patients into different risk groups (high-risk or low-risk). According to RRS, 67.81 and 32.19% of patients were categorized into low-risk and high-risk groups respectively. The relapse rate at 5 years in the low-risk group (2.67, 95% CI: 0.72–4.63%) by Kaplan-Meier method was significantly lower than that of the high-risk group (19.30, 95% CI: 12.34–26.27%) (p <  0.001). In the multiple Cox model, the RRS was proven to be a powerful classifier independent of age at diagnosis or tumour size (p <  0.001). In addition, we found that high RRS score ER-positive patients do not benefit from hormonal therapy treatment (RFS, p = 0.860). Conclusion The RRS model can be applied to predict the relapse risk in early stage breast cancer. As such, high RRS score ER-positive patients do not benefit from hormonal therapy treatment. Electronic supplementary material The online version of this article (10.1186/s12885-019-5941-5) contains supplementary material, which is available to authorized users.


Background
More than 50% of patients with breast cancer are classified into the early-stage (T 1-3 N 0 M 0 ) group [1]. Despite systemic adjuvant therapy dramatically increasing the clinical outcome of patients with early breast cancer, relapse still occurs in more than 20% of patients after surgery within 10 years [2]. Relapse, including recurrence both at local or distant sites, is the main cause for patient deaths, and thus remains an unmet challenge for a curative treatment of breast cancer. It is pivotal to identify patients at risk of relapse at early stages in hopes of improving clinical outcomes, especially within the subgroup of node-negative females, defined as a relatively indolent disease based on pathologic features. Recently, several multigene assays have been developed for earlystage breast cancer patients [3]. Multigene assays are able to provide more prognostic information than traditional parameters in several tumour types [4][5][6][7][8][9][10][11], and several of them have been adopted by the oncology guidelines for treatment. One example is 21-gene expression profiling, which has been widely accepted in clinical practice [12].
As reported, breast cancer is a tumour with high heterogeneity. Although recent advancements have further divided this heterogeneous disease into distinct subgroups by gene expression profiling (GEP) assays, among other methods, several intriguing findings revealed that a small subset of cells isolated from different subgroups of breast cancers exhibit remarkable similar biological behaviours. These subset of cells were defined as cancer stem cells (CSCs) and reported to be responsible for the heterogeneity. Accumulating evidence has proved that CSCs retain the critical characteristics of normal stem cells, such as ability self-renewal and the capacity of proliferation, which contribute significantly to therapeutic resistance and breast cancer relapse [13][14][15][16][17]. In addition, several articles indicated that some CSCs might be derived from normal stem cells, which suggested that normal mammary stem cells might share similar identifying markers [18][19][20]. Mammary stem cell markers or combined markers have been certified in different stages of stem cells in breast cancer, including ALDH, CD44, CD24, ITGA6/EpCAM, and PROCR. [21][22][23][24][25][26]. Some of these markers and combined markers (i.e., CD44 + /CD24 low ALDH + and ITGA6 + ) are considered to correlate with poor prognosis in breast cancer [21,27,28], because they also identified a BCSC subpopulation [14,21,26,29]. In addition, it has been suggested that ITGA6 + /EpCAM + mammary luminal progenitor cells were possible transformation targets in basallike breast cancers, which have close associations with poor prognosis. In addition, it was reported that ITGA6 may define the mesenchymal population and is necessary for CSC function [30][31][32]. PROCR was reported to be highly expressed in myoepithelial cells of the mammary gland. In a recent study, Wang D et al. identified PROCR as a marker of multipotent mammary stem cells. They found that PROCR-positive mammary cells exhibited epithelial-to-mesenchymal transition (EMT) characteristics, and had high tumorigenesis ability in vivo, which suggested that PROCR-positive mammary cells might be one of the progenitor populations for breast CSCs (BCSCs) [24]. Furthermore, PROCR also promotes tumour metastasis in cancer cell lines [33,34].
To explore the prognostic role of mammary stem cell (MSCs) and BCSC markers, we have studied the ALDH family (including ALDH1A1, ALDH1A3, ALDH3A1, ALDH4A1, ALDH6A1, and ALDH7A1), PROCR, and ITGA6/EpCAM. In a medium cohort of patients in previous studies, these findings revealed that ALDH1A3, PROCR, ITGA6 + , ITGA6 + /EpCAM − and ITGA6 − / EpCAM + were correlated with reduced RFS or overall survival of these breast cancer patients [35][36][37]. In this study, we defined these markers and CD44 + /CD24 low as BCSC-associated markers and employed these biomarkers to label stem cells among patients with early stage breast cancer. ALDH1A3, CD44 + /CD24 − , ITGA6, and PROCR were shown to be closely associated with RFS. Then, they were integrated into the prognostic panel to calculate an RRS. Patients were then divided into two distinct risk groups, which effectively shows promise in predicting prognosis and treatment. In addition, several EMT transition associated markers, proliferation factors and other clinicopathological parameters were also included in our study to improve the efficiency of our model.

Materials and methods
Breast cancer patient dataset Clinical information from 1036 patients with breast invasive ductal carcinoma (BIDC) diagnosed from 2006 to 2011 was collected from West China Hospital. After selection, 407 patients were enrolled into our study. All the patients were adult females and were treated with mastectomy or lumpectomy to negative margins and with axillary lymph node dissection. Axillary nodes of patients were observed to be without metastasis under microscope. Patients with local invasion and distant metastasis identified initially were ineligible. Patients with neoadjuvant chemotherapy were removed from our study group to avoid its impact on the characteristics of tumour cells in paraffin embedded tissues. Patients enrolled in the study were considered to be early-stage BIDC and defined as entire datasets. The end-point of follow-up was occurrence of local recurrence or distant metastasis. Detailed information of this dataset is listed in Additional file 4: Table S1.
Detailed information and specificity of these antibodies were shown in Additional file 5: Table S2, Additional file 1: Figure S1, respectively.

Statistical analysis and model construction
The associations between relapse-free survival (RFS) and the expression panel were analysed by the Cox proportional hazard regression model [41]. To investigate the effectiveness of the BCSC-associated biomarker panel for clinical outcome prediction, we assigned each patient a risk score according to a linear combination of the expression level of BCSC-associated markers. The RRS for sample i using the information from the significant biomarkers was calculated as follows: RRS ¼ P 4 j¼1 WjÃSj: In the above formula, Sj is the IHC score for biomarker j, and Wj is the weight of the IHC score of biomarker j. Weights were obtained by the coefficients derived from the univariate Cox proportional hazard regression [42]. The RRS was calculated out by the receiver operating characteristic curve (ROC, non-parametric test), which identifies the cut-off value based on the maximum sums of specificity and sensitivity in the ROC curve. Meanwhile, to investigate the association between the relapse and other clinicopathological variables, univariate Cox proportional hazard regression analysis was adopted using clinicopathological factors (including age, tumour size, histological grade, ER status, PR status and HER2 status), proliferation factors (Ki67), and EMT related factors (including Twist and Slug) in the dataset. The cutoff values of ER, PR, HER2 and Ki67 were 1, 1%, 1+/2+, and 14%, respectively, according to the standards of clinical practice. For twist and slug, the final score was 0 to 12 as the cut-off value for the analyses to obtain significant results. Furthermore, multivariate Cox proportional hazard regression analysis was applied to investigate whether the predictive value of the panel was independent of other clinical variables.
The model was established using the and holdout methods, an approach to out-of-sample evaluation, where the dataset was randomly split into two exclusive sets (70% training and 30% testing sets) [43]. The model was then trained on the training group and tested on the testing group 10 times. Additionally, bootstrapping was used to overcome a possible biased error estimate and obtain confidence intervals (CI). We reported the 95% CI of the coefficients, hazard ratio, and relapse rate for each model. Statistical analyses were performed using GraphPad Prism version 6 and R 3.4.0. To enroll more effective biomarkers and clinicopathological factors into further modelling, a p-value less than 0.1 was defined as statistically significant in the univariate Cox Proportional Analysis. Then, potential significant factors were enrolled into the multivariate Cox Proportional Analysis, with the p-value less than 0.05 considered to be statistically significant. The detail was shown in Additional file 3: Figure S3.

Characteristics of patients and IHC results
The mean age of the patients was 49.3 ± 9.9 years. The youngest patient was 23 years old while eldest one was 78 years old. Among the 407 patients, the median follow-up was 66 months, and relapse was observed in 42 (10.3%) patients during five years after diagnosis, consistent with results published in the literature. The characteristics of clinicopathological, proliferation, and EMT related factors of the 407 patients are depicted in Table 1 and Additional file 4: Table S1. IHC staining was performed on slides of paraffin embedded blocks of those 407 BIDC samples. Results are shown in Fig. 1. We also performed IHC in tissues of patients with reductional mammoplasty. The prevalence of BSCCs biomarkers in reductional mammoplasty samples were shown in Additional file 2: Figure S2.

Construction and validation of the RRS model
A univariate analysis was performed to test whether the expression level of each BCSC-associated marker was related to differences of patient RFS. Among all the BCSC related biomarkers, four biomarkers (ALDH1A3, CD44 + / CD24 − , ITGA6 + , and PROCR) were confirmed to be statistically correlated with patient RFS ( Table 2). The RRS formula according to the expression coefficient of those 4 BCSC-associated biomarkers for survival is listed as follows: RRS = 0.30× (score of ALDH) + 0.34× (score of CD44 + /CD24 − ) + 0.24× (score of ITGA6) + 0.56× (score of PROCR). Therefore, patients were classified into highrisk and low-risk group individually using the optimal RRS (RRS corresponding to the maximum sum of specificity and sensitivity in the ROC curve) as the cut-off value. With the aid of the method described in the Materials and Methods, the cut-off value was calculated to be 2.05.
Then, Kaplan-Meier analysis showed that the proportion of patients in the low-risk group who were free of relapse at 5 years (97.68, 95% CI: 97.37-98.00%) was significantly higher than that in the high-risk group (81.33, 95% CI: 80.50-82.16%) (p < 0.001) in the training group. In another exclusive group (the testing group), the proportion of patients in the low-risk group who were free of relapse at 5 years (96.82, 95% CI: 95.88-97.76%) was also higher than that in the high-risk group (82.13, 95% CI: 79.93-84.33%) (p < 0.001). Distributions of risk score, relapse status and BCSC-associated biomarker expression of patients in the training group and testing group is displayed in Table 3 and Fig. 2.
Among all the clinicopathological factors (including age at diagnosis, tumour size, histological grade, ER status, PR status and HER2 status), proliferation factors (Ki67), EMT related factors (including Twist and Slug), age at diagnosis and tumour size were considered potential significant factors in the univariate survival analysis. These factors were then fully enrolled to the multivariate Cox model with RRS. In a multiple Cox model, RRS demonstrated significant predictive power that was independent of tumour size and age at diagnosis in both the training group (p < 0.001) and testing group (p = 0.014) ( Table 4).

Assessment of the RRS model in the entire dataset Assessment of the RRS model in univariate survival analysis (Kaplan-Meier method)
To validate our findings, the RRS model was assessed in the entire dataset (n = 407). By using the same cut-off value of training groups, patients in the entire dataset were classified into the high-risk group (n = 131) and low-risk group (n = 276) (Fig. 3a). Patients with high risk scores demonstrated significantly reduced RFS when compared to those with low risk scores (log-rank test p < 0.001) (Fig. 3b). The relapse rate at 5 years was 19.30% (95% CI: 12.34-26.27%) and 2.67% (95% CI: 0.72-4.63%) in the high-risk group and low-risk group, respectively. Distributions of risk score, relapse status and BCSC-associated biomarker expression of each patient in the entire datasets were then analysed (Fig. 3c).

Assessment of the RRS model in multivariate survival analysis (cox proportional analysis)
In the entire dataset, the correlation between RFS and clinicopathological factors (including age, tumour size, histological grade, ER status, PR status and HER2 status), proliferation factors (Ki67), EMT related factors (including Twist and Slug) was analysed by Kaplan-Meier method. Reduced RFS was only demonstrated in patients with smaller tumour size (log-rank p = 0.032) and younger age (log-rank p = 0.016) ( Table 1). Then, multivariate survival analyses were adopted to explore the association between relapse and age as well as tumour size. As a result, younger age, larger tumour size and RRS were implied to be significant predictors of relapse (Table 5).

Hormone therapy benefit in different groups
Among the 407 patients, there were 282 ER-positive and 125 ER-negative patients. We found that our panel worked in both of these two subgroups (Fig. 4a, b). In the ER-positive group, all patients were treated with chemotherapy, whereas only 89.72% (n = 253) of these patients received hormone therapy. Our results demonstrated no difference for the RFS between those hormone-treated patients and non-treated patients in the high-risk score group (p = 0.860 Fig. 4d). However, in the low-risk score group, patients in the treated group showed remarkably longer RFS than those in the nontreated group (p = 0.038, Fig. 4c), which indicated that patients with a high-risk score may not benefit from the traditional hormone therapy.

Discussion
An increasing number of females are diagnosed with node negative invasive breast carcinomas. Even though most of patients with early-stage breast cancer have a favourable outcome, the 5-year rate of local relapse or distant metastasis in our dataset is still up to 10.3%. As metastatic diseases are challenging to cure, accurate evaluation for prognosis and more efficacious treatments are needed. In our present study, we developed and validated a novel prognostic model based on 4 BCSC-associated biomarkers to improve our accuracy of predicting disease recurrence in patients with early stage BIDC (T 1-3 N 0 M 0 ). The four biomarkers incorporated into our predictive model have been shown to be involved in stem cell ability in vivo and in vitro, including self-renewal ability and tumorigenic capacity, which could contribute greatly to metastasis of BIDC in vitro and in vivo, or in tumour tissues [21][22][23][24][25][44][45][46]. The holdout methods were adopted to establish our RRS model, which assisted us to obtain a stable model to calculate RRS in our study. Our model was further validated in the entire dataset. The AUC value of ROC curve is 0.781 which indicated that the RRS is a good classifier for relapse among patients with early stage breast cancer. The difference in the risk of relapse between patients with low risk scores and those with highrisk scores was large and statistically significant. There are 276 (67.81%) patients who were classified in the lowrisk group, while only 32.19% of patients were included    in the high-risk group, and their rate of relapse at 5 years was 19.30 and 2.67%, respectively. Therefore, the application of the RRS predictor provides a good estimate of the risk of local or distant recurrence in individual patients.
We also enrolled other biomarkers in the univariate survival analysis in the training set, such as age, tumour size, histological grade, Ki67, and EMT related biomarkers. All those parameters have been reported to play critical roles in accelerating the presence of distant metastasis or local relapse [47,48]. Despite the fact that EMT has been reported to produce cells with stem celllike properties [49], we found that no parameter showed significantly different RFS in different subgroups of EMT related biomarkers. In this study, smaller tumour size was validated as an independent factor protecting patients from relapse. When the RRS was combined  with data pertaining to tumour size to predict the risk of relapse, the relapse score remained statistically significant in a multivariate analysis. Due to poor compliance of our patients, in the ERpositive subgroups, only 89.72% of patients received endocrine therapy systematically. The results indicated that only patients with low risk responded well to endocrine therapy, while those with high risk showed no difference between the treated group and untreated group. A previous study revealed that mesenchymal-like BCSCs in hormone-sensitive luminal breast cancers were one of the reasons for hormone-resistant [50]. Similar to above finding, there was evidence suggesting that BCSCs should be partially responsible for the endocrine-resistant capacity of breast cancer cells. This is due to the fact that CSCs could only respond to treatment by virtue of paracrine signalling pathway from adjacent differentiated ER-positive tumour cells [51][52][53][54], which were probably responsible for the endocrine-resistance in the high-risk group.
The RRS not only offers an approach to predict therapeutic sensitivity but also provides a new perspective to eliminate BCSCs in early stage breast cancer. As been reported, BCSCs were not as sensitive to hormone therapy and conventional chemotherapy as non-BCSC tumours. Thus, targeting BCSCs clinically might enhance the therapeutic sensitivity among patients with high risk scores. The most promising CSC treatment strategies that target Notch, Hedgehog, Wnt and many other BCSC self-renewal pathways provide a number of opportunities for new clinical trials. 20 In addition, the strategy of "destemming" CSCs, including inducing Fig. 4 Kaplan-Meier analysis for RFS using RRS in the subgroups stratified by ER status and endocrine therapy. a Kaplan-Meier curves for earlystage BIDC patients with ER-positive status. b Kaplan-Meier curves for early-stage BIDC patients with ER-negative status. c Kaplan-Meier curves for ER-positive patients with high risk scores stratified by endocrinotherapy. d Kaplan-Meier curves for ER-positive patients with low risk scores stratified by endocrinotherapy CSC differentiation or inhibiting self-renewal capacity were also recommended [55]. Combination of BCSCtargeted therapy and traditional therapy may provide our patients with high-risk scores more effective therapeutic strategies. However, the study of CSCs remains an enigma, and further exploration is needed.
In terms of limitations, this study was a retrospective analysis that selected patients who had not received neoadjuvant chemotherapy after resection in early stage breast cancer, which may lead to a selection bias of patients with a relative lower risk of recurrence. However, all our patients included in this study were T 1-3 N 0 M 0 by the TNM staging system, and the majority of them did not receive neoadjuvant chemotherapy, according to the NCCN guideline [12]. The total study size is modest in absolute numbers, and some subgroup analyses may be underpowered; however, this is one of the largest cohorts of well-characterized early stage breast cancer that employed a BCSC biomarker panel as a prognosis model. The shortcomings of this panel should not be ignored. First of all, though IHC staining is the most common method for semi-quantified the protein expression level in carcinomous tissues, the subjectivity of evaluation of this method couldn't be avoided. Secondly, the selection of antibodies should be cautiously considered, as their quality will affect the result of IHC staining directly. Performing immunofluorescence staining and q-RT PCR may help us obtain a relative exact result; however, these two methods also have their disadvantages in assessing BCSCs.

Conclusion
Though previous studies have combined different BCSCs biomarkers for assessing prognosis in different types of breast cancer, such as three-negative, HER2positive and metastatic breast cancer [56][57][58][59], no BCSC-associated biomarkers have been combined to form a model for evaluating the relapse risk of earlystage breast cancer. We propose that BCSCs could be used as a panel in prognostic or predictive tests of early-stage breast cancer. Here, we conducted a prospectively designed validation study of a multi-biomarker panel in a cohort of patients with early-stage BIDC. In addition, this panel is promising for prediction of early-stage BIDC recurrence, the efficacy of which warrants further validation in a large-scale cohort. In addition, it reminds us that further consideration is needed to explore new therapeutic managements for high-risk patients with therapeutic resistance. In addition, it is of practical significance that the panel only involves the use of routine slides of the tumour tissues and five antibodies, which is not as time-consuming and expensive as other gene profiles.