A glycolysis-related signature to improve the current treatment and prognostic evaluation for breast cancer

Introduction

As a frequently diagnosed female malignancy (Siegel et al., 2021; Shang & Xu, 2022), breast cancer (BRCA) accounts for 15% (Heer et al., 2020) of all cancer-correlated deaths (Nik-Zainal et al., 2016). Invasive and metastatic characteristics (Kim et al., 2018) of BRCA are the leading causes of death to women aged 40 to 79 years old and males aged 60 to 79 years old (Siegel et al., 2023). The latest National Cancer Statistics Report showed that there were 297,790 and 59,910 new cases and deaths of BRCA in 2023, respectively (Siegel et al., 2023), exhibiting an obvious increasing trend as compared to 268,600 new cases and 41,760 death cases in 2019 (Siegel, Miller & Jemal, 2019). The regulating molecular mechanisms of the tumorigenesis and progression of BRCA are decisive to the heterogeneity of the cancer, which also poses great challenge to clinical treatment and disease prognosis for BRAC patients (Lüönd, Tiede & Christofori, 2021). High-throughput sequencing promotes the identification of molecular markers and development of effective treatment and intervention strategies for BRAC patients. This also suggests that developing reliable prognostic models can improve the overall survival (OS) of BRCA patients.

Metabolic reprogramming is an emerging hallmark of cancers. “Aerobic glycolysis” or the “Warburg effect” refers to the phenomenon that cancer cells favor glycolysis over mitochondrial oxidative phosphorylation in the breakdown of glucose even there is plentiful oxygen. A widely accepted explanation for this metabolic shift in cancer cells is that enhanced glycolysis is utilized by dividing cells as a flexible assembly line to produce glycolytic intermediates for multiple biosynthetic pathways to support rapid cell growth. In hepatocellular carcinoma, Li et al. (2017) demonstrated that genistein inhibits aerobic glycolysis through inactivating hypoxia inducible factor-1α (HIF-1α) to downregulate glucose transporter (GLUT1) and HexokinaseII. This suggested that developing glycolysis-related signatures may be a novel direction for cancer treatment and prognostic evaluation. Zhang, Zhang & Yu (2019) identified four glycolysis-related genes (DDIT4, AGRN, HMMR and AKR1A1) closely associated with the survival outcome of lung cancer patients. The glycolytic signature model established by Wang et al. (2019) can accurately predict the survival outcome of patients with endometrial cancer. This study contributed to the previous research on the role of glycolytic genes in BRCA by exploring the correlation between glycolysis and BRCA prognosis based on public databases. A prognostic model was developed with glycolysis-related genes and the correlation between risk score and immune infiltration landscape was analyzed. Furthermore, the impact of the key risk factors on the biological functions of BRCA cells was comprehensively explored. In conclusion, our study provided a novel direction for predicting clinical outcomes in BRCA patients.

Materials and Methods

Results

Glycolysis was closely related to the progression of BRCA

The hallmark_glycolysis score was calculated for each tumor and para-cancer sample from TCGA, and it was found that tumor samples has higher glycolysis score (Fig. 1A). Based on the optimal survival cutoff of glycolysis score, the samples in the TCGA-BRCA and GSE20685 dataset were categorized into high- and low-glycolysis groups. The KM survival analysis revealed worse prognosis (Fig. 1B), disease free survival (DFS, Fig. 1C) and progress free survival (PFS, Fig. 1D) and poor prognosis (Fig. 1E) of patients in high-glycolysis groups. Moreover, carcinogenic signaling pathway score for each TCGA-BRCA patient was calculated. Analysis on the correlation between the glycolysis score and carcinogenic signaling pathway score showed that the glycolysis score was significantly positively correlated with multiple carcinogenic signaling pathways, especially the PI3 kinase Akt signaling (PI3K), epidermal growth factor receptor (EGFR), hypoxia, vascular endothelial growth factor (VEGF) pathway, and tumor necrosis factor-α (TNFα) (Fig. 1F) and was significantly negatively correlated with the MCP-counter and TIMER immune scores of multiple immune cells such as the Neutrophils, NK cells and myeloid dendritic cells (Fig. 1G), macrophage, CD8_T cells and CD4_T cells (Fig. 1H). This indicated that overactive glycolysis could inhibit the functions and number of immune cells.

Construction of a glycolysis-related risk model

The DEGs were identified between high-glycolysis and low-glycolysis groups in the TCGA-BRCA and GSE20685 datasets (Figs. 2A, 2B), and there were 299 overlapping genes (Fig. 2C). Pathway enrichment analysis revealed that these overlapping genes were mainly enriched in cell cycle, progesterone-mediated oocyte maturation, and oocyte meiosis (Fig. 2D), indicating that glycolysis promoted the proliferation of cancer cells. Next, univariate Cox regression and LASSO Cox regression analysis were performed to determine the key prognostic genes from these 299 genes (Fig. 3A). Multivariate Cox regression analysis eventually selected five prognostic genes for BRCA and their regression coefficient was calculated (Fig. 3B). The formula of the model was RiskScore $=\left(-0.102\ast \mathrm{F}\mathrm{R}\mathrm{E}\mathrm{M}1\right)+\left(-0.051\ast \mathrm{N}\mathrm{P}\mathrm{Y}1\mathrm{R}\right)$ + $\left(-0.075\ast \mathrm{C}\mathrm{C}\mathrm{L}19\right)+0.074\ast \mathrm{P}\mathrm{S}\mathrm{C}\mathrm{A}+\left(-0.061\ast \mathrm{C}\mathrm{L}\mathrm{I}\mathrm{C}6\right)$. Each patient was assigned with a RiskScore and the optimal survival cutoff of the RiskScore could effectively distinguish between high- and low-risk patients. KM survival analysis revealed that the overall survival (OS) of high-risk patients were significantly worse (p < 0.0001, Fig. 3C). ROC analysis demonstrated that the AUC value of 1-, 3- and 5-year survival was 0.75, 0.71 and 0.65, respectively (Fig. 3D), showing a strong classification performance of the RiskScore. Comparison on the expression levels of the five model genes in different risk groups showed that only the high-risk patients had high-expressed PSCA (Fig. 3E), which was consistent with the results from multivariate Cox regression analysis as it showed that PSCA was a risk factor (Hazard ratio > 1 and p < 0.02). In addition, the OS of high-risk patients in the validation set (GSE20685) was significantly worse as compared to the low-risk patients (Fig. 3F), with 1-, 3- and 5-year AUC value of 0.75, 0.75 and 0.69, respectively (Fig. 3G). PSCA was high-expressed in high-risk patients as a risk factor (Fig. 3H). Based on the expression of the five model genes, principal component analysis (PCA) showed that the five genes could effectively divide the patients in both TCGA and GSE20685 cohorts into two risk groups (Figs. 3I, 3J).

Screening independent factors and development of a nomogram

The RiskScore increased with lower clinical grade in the TCGA and GSE20685 cohorts (Figs. 4A, 4B). Univariate Cox regression analysis showed that the age, RiskScore, T. stage, N. stage, M. stage, and stage were significant prognostic factors for BRCA (p < 0.001, Fig. 4C). Multivariate Cox regression analysis further showed that the age, M. stage, stage, and RiskScore were significantly independent prognostic factors (p < 0.001, Fig. 4D). A nomogram model for predicting 1-, 3- and 5-year survival was developed (Fig. 4E). The calibration curves of the nomogram at 1-, 3- and 5-year points was highly consistent with the standard curves (Fig. 4F), indicating an accurate prediction of the nomogram. Also, the decision curve displayed a significantly higher benefit of the nomogram than the extreme curve (Fig. 4G), and the nomogram had strong prediction performance than other clinicopathological features.

Immune infiltration differences among different risk groups

In the TCGA cohort, the low-risk group showed higher stromal score, immune score and ESRIMATE score (Fig. 5A). MCP-counter analysis demonstrated that this low-risk group was associated with the higher infiltration levels of immune cells, such as the CD8 T cells, NK cells, T cells, B lineage cells, endothelial cells, and myeloid dendritic cells (Fig. 5B). The results of CIBERSORT calculation revealed that the low-risk group had significantly more CD8+ T cells, M1 macrophages, NK cells activated, and mast cells (Fig. 5C). Moreover, the immune infiltration in the validation set (GSE20685) was also quantified. The results demonstrated that the low-risk samples were associated with higher stromal score and immune score (Fig. 5D) and had higher infiltration levels of lineage cells, myeloid dendritic cells, T cells, CD8 T cells, NK cells (Figs. 5E, 5F). These finding suggested that the abundance of immune cells were significantly reduced in the high-risk group and their normal functions were also inhibited in TME.

Immunotherapy and drug sensitivity analysis for different risk groups

The low-risk group with lower IPS were more likely to benefit from receiving immunotherapy treatment (Fig. 6A). Meanwhile, comparison on the antigen presentation of each risk group showed that the expression levels of related genes (such as the HLA-DPA1, HLA-DPB1 and HLA-DQA1) were higher in the low-risk group (Fig. 6B), while those of immune checkpoint inhibitor-related genes (such as the PDCD1, BTN2A2, VTCN1, ADORA2A, CTLA4 and TIGIT) were significantly higher in low-risk patients (Fig. 6C). Above these results indicated greater immunotherapy benefit for low-risk patients. The drug sensitivity analysis revealed that PHA-665752, DMOG, AZ628, AP-24534, Crizotinib, CP466722 and HG-6-64-1 were positively correlated with the low-risk group (Fig. 6D), and that pyrimethamine was negatively correlated with the low-risk group and had lower IC50 value in the high-risk patients (Figs. 6D, 6E). This suggested that pyrimethamine treatment could benefit high-risk patients more.

PSCA was a key risk factors affecting the proliferation, invasion and migration of BRCA cells

Our model revealed that high-expressed PSCA was related to a higher RiskScore, and the function of PSCA was further verified using tumor cells. In the Western blot assay, the protein expression in the siPSCA#1 and siPSCA#2 groups was significantly downregulated (p < 0.05) as compared to the normal control group (Fig. 7A), indicating that the two independent small interfering RNAs had a high silencing efficiency. The transwell assay showed that the cell number was significantly reduced after silencing PSCA (Fig. 7B), suggesting a suppressed cell invasion activity. Additionally, in the colony formation assay, the number of cell colony numbers in the PSCA silencing groups were significantly less than that in the siNC groups (Fig. 7C), which indicated an inhibited cell proliferation ability. Meanwhile, the wound-healing assay showed that the cell migration ability was suppressed after PSCA silencing (Fig. 7D). These results validated that the PSCA was a key risk factors affecting the cell proliferation, invasion and migration in BRCA, and that high-expressed PSCA promoted BRCA progression.

Discussion

Recent findings in the studies of energy metabolism extend our traditional understanding of cancer as primarily a genetic disease to a metabolic disease (Upadhyay et al., 2013). The distinct metabolic characteristic of cancer cells, alternatively known as the “Warburg effect”, means that cancer cells tend to favor glycolysis over oxidative phosphorylation even when oxygen is abundant (Vaupel, Schmidberger & Mayer, 2019). Emerging evidence indicates that the majority of cancer cells exhibit the “Warburg” metabolic phenotype, leading to a faster production of ATP compared to oxidative phosphorylation (Graziano et al., 2017; Yan et al., 2017; Guo et al., 2021). This study selected and analyzed genes related to glycolysis applying a series of bioinformatics analyses. Subsequently, a risk assessment model was developed utilizing multivariate and LASSO Cox regression analyses to select potential biomarkers for patients suffering from BRCA. The prognostic prediction of the model was highly accurate, highlighting its potential significance in the clinical diagnosis and management of BRCA patients.

In recent years, several studies have developed risk models for cancers using glycolysis-related genes. Based on TCGA and GEO databases and univariate Cox and LASSO regression analyses, Guo et al. (2021) constructed a risk model with five glycolysis-related genes to predict the prognosis for patients with prostate cancer. Liu, Wang & Li (2022) also screened genes related to glycolysis and immunity from public databases and developed a risk score model consisting of six glycolysis-immunity-related genes (ALDOC, VEGFA, PADI3, IGSF11, MIPOL1, and HGR) for oral squamous cell carcinoma. In addition, another previous study also screened six genes as independent prognostic genes for colon adenocarcinoma by performing enrichment analysis on the glycolysis-associated gene set but this model was not suitable for patients with rectal adenocarcinoma (Liu et al., 2022). These findings suggested a strong potential of using glycolysis-related genes as potential targets for cancer therapy.

This study identified PSCA as a risk gene (hazard ratio > 1) for BRCA, while FREM1, NPY1R, CCL19 and CLIC6 were the protective genes for the cancer (hazard ratio < 1). PSCA encodes cell membrane glycoprotein involved in signal transduction, cell adhesion and proliferation (Wu et al., 2020). This antigen has been reported as a critical marker in several cancers, such as prostate, bladder and pancreatic cancers (Saeki et al., 2010), and its high expression is related to a worse survival (Kang et al., 2016). A stable knockdown of PSCA can delay the growth, proliferation, metastasis and invasion of prostate cancer cells (Liu et al., 2017). In addition, the expression of PSCA is restricted to the tumor cells of BRCA and is upregulated in invasive BRCA but there is no clear association between PSCA-expression and treatment outcome, suggesting that PSCA cannot directly serve as a prognostic marker for patients with BRCA (Link et al., 2017). Consistently, this study also found that the PSCA was a key BRCA-promoting factor by enhancing the proliferation, invasion and migration ability of BRCA cancer cells, and that PSCA might be a potential treatment target. FREM1 encodes extracellular matrix (ECM) protein and plays a crucial part in the adhesion process of basement membrane and dermis (Petrou, Makrygiannis & Chalepakis, 2008). Previous study analyzed FREM1 expression profile with immunohistochemistry (IHC) and observed that a lower expression of FREM1 in BRCA is related to lower immune cell infiltration and worse overall survival (OS) (Li et al., 2020). Meanwhile, using various machine learning algorithms, Xu et al. (2024) screened FREM1 as one of the core genes for smog disease and found that the high FREM1 subgroup had significantly enriched glycosaminoglycan degradation and glycosphingolipid biosynthesis-ganglio series. NPY1R, an estrogen‑responsive gene, is high-expressed in BRCA patients and mediates the inhibitory action to tumor cells. High expression of NPY1R as an independent prognostic factor is normally associated with higher endocrine sensitivity and better survival (Bhat et al., 2022). Based on ex vivo and in vivo studies, Choong et al. (2021) indicated that influencing islet peptide tyrosine tyrosine (PYY (1–36)) through NPY1R dependence improves glucose homeostasis in patients with type 2 diabetes. Chemokine CCL19 plays a pivotal role in tumor-suppressive functions (Bhat et al., 2021), for example, CCL19 promotes naïve T-cell survival via its anti-apoptotic activity (Link et al., 2007), enhances the endocytosis capacity of dendritic cells (DCs) (Wu et al., 2023), and induces the production of TNFα, IL-12 and IL-1β, indirectly leading to increased proliferation of T cells and natural killer (NK) cells (Laufer & Legler, 2018). Hayashi et al. (2021) observed that overexpression of adipose-specific CCL19 can drive insulin resistance and weight gain in mice. In addition, CLIC6 is a channel protein involved in cell proliferation and homeostasis maintenance at early tumor stage, and its expression level is downregulated as malignant development progresses (Shu et al., 2024). However, the significance of CLIC6 in the glycolysis in cancer patients remained unknown, and our findings suggested that this gene was implicated in the glycolysis and progression of BRCA.

Overall, we constructed a RiskScore model for the prognostic prediction of BRCA based on the gycolytic metabolism characteristics. These key candidate genes in the RiskScore were associated with the adhesion, proliferation, homeostasis maintenance and immune cell protection process. High-risk BRCA patients exhibited lower immune infiltration and immunotherapy evasion, indicating that the immune microenvironment homeostasis was disrupted and recruited immune cells were reduced. For high-risk patients, pyrimethamine may be a more effective drug. Our study provided novel molecular targets for BRCA treatment but still had some limitations. First, the sample size of the dataset used in this study was relatively small. Therefore, it is important to collect more samples of BRCA patients, including those from different regions and ethnicities. In addition, further animal model experiments should be carried out to verify the role of the key genes in BRCA progression.

Conclusions

The current study developed a reliable five-gene glycolysis-related prognostic model as an independent prognostic factor for BRCA. High-risk patients were closely correlated with lower immune infiltration, higher immunotherapy evasion and worse prognosis. In addition, a nomogram model integrating the RiskScore, age and stage can accurately predict both long- and short-term survival for BRCA patients with different RiskScores. PSCA was a key risk factor that affected the cell proliferation, invasion and migration of BRCA cells. It is hoped that the current model can become a reliable clinical tool to help clinicians make treatment decisions for BRCA patients with different risks.

Supplemental Information