A Cori cycle-related gene signature predicts prognosis, immune microenvironment, and drug response in breast cancer

Introduction

Breast cancer is the most prevalent and burdensome malignancy among women worldwide. Recent regional reports indicate that in China, breast cancer has surpassed lung cancer for the first time as the most commonly diagnosed cancer in women, accounting for approximately 19.9% of all new cancer cases (1). Although significant advances have been made in screening, diagnosis, and comprehensive treatment—including surgery, radiotherapy, endocrine therapy, and the widespread application of targeted agents—treatment outcomes and long-term survival remain highly variable due to the disease’s substantial molecular heterogeneity and clinical diversity (2). Existing prognostic models, such as those based on tumor, node, metastasis (TNM) staging, histological grade, and receptor status, have notable limitations. For instance, TNM staging primarily focuses on anatomical factors and often fails to account for biological variability, leading to suboptimal risk stratification in heterogeneous subtypes like triple-negative breast cancer (TNBC), where recurrence rates remain high despite similar staging (3). Similarly, multigene assays like Oncotype DX and MammaPrint, while useful for hormone receptor-positive cases, exhibit reduced accuracy in human epidermal growth factor receptor 2 (HER2)-enriched or TNBC subtypes and do not integrate metabolic or immune features, potentially overlooking therapy resistance mechanisms (4). Therefore, the discovery of novel molecular biomarkers and the development of more accurate prognostic evaluation systems are of great importance for advancing personalized therapy and improving patient survival.

Metabolic reprogramming is a well-established hallmark of cancer, playing a pivotal role in tumorigenesis, progression, and therapy resistance (5). Among metabolic pathways, the Cori cycle, which facilitates lactate recycling in the body, serves as a crucial mechanism by converting lactate produced in peripheral tissues back to glucose in the liver, thereby maintaining systemic energy homeostasis (6). Recent studies have shown that this cycle is aberrantly activated in multiple malignancies, not only supporting rapid tumor proliferation under hypoxic conditions but also modulating immune cell functions and the metabolic state of the tumor microenvironment (TME), thereby influencing disease progression (7). Cori cycle-related genes (CCRGs) encode key metabolic enzymes and transporters, and their dysregulation is closely associated with glycolytic flux (8), lactate metabolism (9), and redox balance (10), suggesting their potential involvement in breast cancer malignancy and treatment response. However, the overall prognostic value, immunomodulatory mechanisms, and clinical translational potential of CCRGs in breast cancer have not been systematically evaluated.

This study aims to integrate multi-omics data with machine learning approaches to construct a prognostic risk model based on CCRGs in breast cancer. We comprehensively assess its predictive performance and biological foundations, and further explore its implications for immunotherapy response, drug sensitivity, and pan-cancer relevance. Our findings may provide new directions for metabolically targeted research and personalized treatment in breast cancer. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2013/rc).

Methods

Data acquisition and preprocessing

Microarray expression profiles of breast cancer samples were downloaded from the Gene Expression Omnibus (GEO) under accession numbers GSE20685 and GSE21653. Bulk RNA-seq data from The Cancer Genome Atlas (TCGA)-breast invasive carcinoma (BRCA) cohort were obtained via the University of California, Santa Cruz (UCSC) Xena platform. For microarray data, probes were annotated to gene symbols based on the GPL570 platform. Probes matching multiple genes were excluded. When multiple probes corresponded to a single gene, the average expression value was used. A total of 326 and 240 breast cancer samples were included from GSE20685 and GSE21653, respectively. For TCGA-BRCA cohort, samples with survival time >0 and complete survival information were retained. Ensembl IDs were converted to gene symbols, resulting in 482 breast cancer samples from TCGA. Common genes across TCGA and GEO datasets were retained for integration, and batch effects were corrected using the ComBat method from the sva package. Efficacy of batch correction was confirmed via principal component analysis (PCA) visualization. CCRGs were defined and screened based on the following criteria: we obtained a curated list of 17 genes from the WP_CORI_CYCLE gene set in the Molecular Signatures Database (MSigDB), which encompasses genes encoding enzymes and transporters involved in the Cori cycle, including both the glycolytic phase and the gluconeogenic phase. These genes were selected for their documented roles in lactate metabolism, glycolytic flux, and systemic energy homeostasis under normal and pathological conditions (Table S1). This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Machine learning-based model construction

To develop a robust prognostic signature, we employed an integrated machine learning framework, Mime1. The 17 CCRGs from the “WP_CORI_CYCLE” set were used as input. The following steps were executed: (I) the TCGA-BRCA cohort was used as the training set, while all three cohorts (TCGA-BRCA, GSE20685, GSE21653) were used for validation; (II) univariate Cox regression was applied for preliminary screening of candidate genes; (III) the ML.Dev.Prog.Sig function was executed in “all“ mode, integrating multiple algorithms. The random forest node size was set to 5, and a fixed seed (seed =123) was used to ensure reproducibility. Model performance was comprehensively evaluated. The concordance index (C-index) was calculated across all cohorts to assess overall predictive accuracy. For the optimal model, the median risk score (RS) was used as the cutoff for stratifying patients into high- and low-risk groups. Kaplan-Meier analysis with confidence intervals was performed to visualize survival differences. Time-dependent receiver operating characteristic (ROC) curves were plotted at 1, 3, and 5 years, and area under the curve (AUC) values were calculated using Kaplan-Meier estimation methods to quantify the model’s discriminative ability at different time points.

Immune infiltration analysis

Multiple deconvolution algorithms, including Cellular Infiltration By Estimating Relative Subsets Of Transcriptomes (CIBERSORT), Molecular Context Profiling (MCP)-counter, Tumor Immune Estimation Resource (TIMER), quantitative Immune sequencing (quanTIseq), immunophenoscore (IPS), and Expression Signature for Tumor Infiltrating Stromal and Immune cells Method And Toolbox for Evaluation (ESTIMATE), were applied using the IOBR package to estimate immune cell infiltration levels from TCGA-BRCA gene expression data. Spearman correlation analysis was conducted between immune infiltration estimates, the RS, and the expression of signature genes. Results were visualized in a composite heatmap, with separate panels for each algorithm and significance indicators. Correlations between signature genes/RS and immune checkpoint gene expression were also analyzed and visualized. The predictive performance of the RS was further validated in the IMvigor210 cohort of patients receiving anti-programmed death-ligand 1 (PD-L1) immunotherapy. RS for IMvigor210 samples were calculated using the SuperPC model trained on TCGA data. Survival differences between high- and low-risk groups were assessed using Kaplan-Meier analysis and log-rank tests.

Drug sensitivity analysis

Drug sensitivity for TCGA-BRCA samples was predicted using the pRRophetic package, which employs ridge regression models trained on Cancer Genome Project (CGP) data. Sensitivity predictions were generated for 45 anticancer drugs. Pearson correlation analysis was performed between predicted drug sensitivity values and the calculated RS. Results for clinically relevant breast cancer drugs were visualized in a multi-panel scatter plot showing correlation coefficients and significance values.

Enrichment analysis

To elucidate biological functions and pathways associated with differentially expressed genes (DEGs) between high- and low-risk subgroups, enrichment analysis was performed. DEGs were identified using the limma package on TCGA-BRCA RNA-seq count data, with samples grouped by median RS. Genes with a P value <0.05 and |log₂fold change| >1 were considered significantly differentially expressed. A total of 931 DEGs were identified, including 600 upregulated and 331 downregulated genes. Gene Ontology (GO) enrichment analysis for biological processes (BPs), molecular functions (MFs), and cellular components (CCs), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were conducted using clusterProfiler. Terms with an adjusted P value <0.05 [Benjamini-Hochberg (BH) method] were considered significantly enriched. The top significant terms were visualized using dot plots. Gene Set Enrichment Analysis (GSEA) was performed using the Hallmark gene sets (h.all.v2025.1.Hs.symbols.gmt) to identify pathways showing coordinated differential expression. GSEA was run with 10,000 permutations and a minimum gene set size of 10. Results were visualized to show the normalized enrichment score (NES) and running enrichment score for selected gene sets.

Somatic mutation analysis

Mutation data for high- and low-risk groups were extracted using the subsetMaf function from the maftools package. The mafCompare function was used to systematically compare mutational profiles between groups, with a minimum mutation count of 5 to filter low-frequency events. Odds ratios (ORs) and P values (Fisher’s exact test) were calculated to identify significantly differentially mutated genes. The top 20 genes with the most significant differences were visualized using coBarplot, with normalization to compare mutation proportions between groups. The correlation between tumor mutational burden (TMB) and RS was evaluated: TMB was calculated per sample using the tmb function (mutations per megabase), and Spearman correlation analysis was applied. Differences in TMB between high- and low-risk groups were assessed using the Wilcoxon rank-sum test.

Nomogram construction and evaluation

A nomogram integrating the RS with key clinical parameters was developed and validated using the TCGA-BRCA cohort. Clinical variables included age, pathologic T/N stage, overall American Joint Committee on Cancer (AJCC) stage, Prediction Analysis of Microarray 50 (PAM50) subtype, receptor status, and radiotherapy history. Univariate Cox regression was performed to identify variables significantly associated with overall survival (OS), followed by multivariate Cox regression to identify independent prognostic factors. The final nomogram was constructed using the rms package based on the multivariate Cox proportional hazards model, incorporating the RS and significant clinical predictors. The nomogram was visualized with points assigned to each variable value and total points corresponding to 1-, 3-, and 5-year survival probabilities. Internal validation was performed via bootstrap resampling (1,000 repetitions) to generate calibration plots comparing predicted and observed survival rates at 1, 3, and 5 years. Decision curve analysis (DCA) was conducted using the rmda package to evaluate the clinical utility of the nomogram compared to individual clinical factors by calculating net benefit across threshold probabilities. Time-dependent ROC analysis was performed using the survivalROC package to assess the discriminative ability of the nomogram at 1, 3, and 5 years.

Pan-cancer analysis

A systematic pan-cancer analysis of the WP_CORI_CYCLE gene set was performed using the TCGAplot R package and multi-omics data from 33 cancer types in TCGA. Analyses included: differential expression of the gene set between tumor and adjacent normal tissues across cancer types (visualized with paired boxplots); association between gene set expression and TMB and microsatellite instability (MSI) (visualized with radar plots); and prognostic analysis using forest plots to estimate hazard ratios (HRs) for the gene set across cancer types. All visualizations were generated using default parameters of the TCGAplot functions.

Single-cell expression and functional analysis

Single-cell expression levels of CCRGs were obtained from the stomics database (https://db.cngb.org/stomics/analysis/) using datasets STDS0000027, STDS0000038, and STDS0000049. Functional annotations of these genes were analyzed using four breast cancer single-cell RNA-seq datasets (EXP0052, EXP0053, EXP0054, EXP0055; GEO accessions: GSE77308, GSE75688, GSE75367, GSE86978) from the CancerSEA database. Correlations between gene expression and 14 functional states (e.g., angiogenesis, apoptosis, cell cycle) were visualized via heatmap.

Statistical analysis

All statistical analyses were conducted using R (v4.5.0). Survival differences were assessed via Kaplan-Meier analysis with log-rank test. Model performance was evaluated using time-dependent ROC curves, AUC, and the C-index. Group comparisons employed Wilcoxon rank-sum or Kruskal-Wallis tests. Correlation analyses used Spearman or Pearson methods. A two-sided P value <0.05 denoted statistical significance.

Results

CCRGs are associated with breast cancer prognosis

We first integrated three breast cancer cohorts after removing batch effects for subsequent model construction and evaluation. PCA demonstrated that the batch effects were effectively eliminated (Figure 1A). We then analyzed the somatic mutation landscape of CCRGs in the TCGA-BRCA cohort. Only phosphofructokinase, platelet (PFKP), aldolase, fructose-bisphosphate A (ALDOA), and phosphoglycerate kinase 1 (PGK1) exhibited mutation frequencies above 1%, predominantly missense mutations (Figure 1B). Differential expression analysis revealed that solute carrier family 2 member 4 (SLC2A4) and glutamic-pyruvic transaminase (GPT) were significantly downregulated, while solute carrier family 2 member 1 (SLC2A1) was significantly upregulated in breast cancer tissues (Figure 1C). Univariate Cox analysis identified PGK1 [HR =3.02, 95% confidence interval (CI): 1.86–4.89, P<0.001], PFKP (HR =1.38, 95% CI: 1.03–1.84, P=0.03), hexokinase 1 (HK1) (HR =1.70, 95% CI: 1.14–2.53, P=0.009), and phosphoglycerate mutase 1 (PGAM1) (HR =2.35, 95% CI: 1.11–4.98, P=0.03) as risk factors for OS in the TCGA-BRCA cohort, as well as for disease-specific survival (DSS). Additionally, HK1 (HR =1.55, 95% CI: 1.01–2.38, P=0.043) was also a risk factor for progression-free interval (PFI) in this cohort (Figure 1D).

Figure 1 Characteristics of CCRGs in breast cancer. (A) PCA plot of the three integrated breast cancer cohorts after batch effect removal. (B) Somatic mutation landscape of CCRGs in the TCGA-BRCA cohort. (C) Differential expression analysis of CCRGs in the TCGA-BRCA cohort. (D) Univariate Cox analysis of the association between CCRGs and patient OS, DSS, DFI, and PFI in the TCGA-BRCA cohort. BRCA, breast invasive carcinoma; CCRG, Cori cycle-related gene; DFI, disease-free interval; DSS, disease-specific survival; OS, overall survival; PCA, principal component analysis; PFI, progression-free interval; TCGA, The Cancer Genome Atlas.

Construction of a CCRG-based risk signature for breast cancer

To investigate the potential of CCRGs in breast cancer risk prediction, we employed multiple machine learning models. Three CCRGs significantly associated with OS—HK1, PGK1, and PGAM1—were selected for model construction using the integrated cohorts. Several models, including survival-SVM, SuperPC, and stepCox + SuperPC, demonstrated consistent performance (Figure 2A). The SuperPC model was chosen for subsequent analysis. Patients across all cohorts were assigned RS based on this model and stratified into high- and low-risk groups using the median score. Kaplan-Meier analysis revealed significantly worse OS in the high-risk group (Figure 2B), with HRs of 3.75 (P<0.001) in TCGA-BRCA, 1.84 (P=0.006) in GSE20685, and 1.91 (P=0.005) in GSE21653. The model achieved AUC values for predicting 1-year OS of 0.861, 0.513, and 0.710 in the three cohorts, respectively; for 3-year OS, AUCs were 0.749, 0.560, and 0.606; and for 5-year OS, AUCs were 0.750, 0.574, and 0.607 (Figure 2C).

Figure 2 Construction of a prognostic risk signature based on CCRGs in breast cancer. (A) Distribution of C-index values for various machine learning models based on CCRGs. (B) Kaplan-Meier survival curves comparing OS between high- and low-risk groups defined by the SuperPC model. (C) ROC curves evaluating the accuracy of the SuperPC model in predicting 1-, 3-, and 5-year OS. AUC, area under the curve; BRCA, breast invasive carcinoma; CCRG, Cori cycle-related gene; CI, confidence interval; GBM, gradient boosting machine; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; OS, overall survival; PC, principal component; ROC, receiver operating characteristic; RSF, random survival forest; SVM, support vector machine; TCGA, The Cancer Genome Atlas.

To demonstrate the added value of our three-gene signature, we compared its prognostic performance with individual genes and existing glycolysis-related signatures. Univariate Cox models for single genes in the TCGA-BRCA cohort yielded lower AUC values (Figure S1): HK1 (1-year: 0.472, 3-year: 0.577, 5-year: 0.655), PGK1 (1-year: 0.868, 3-year: 0.756, 5-year: 0.736), and PGAM1 (1-year: 0.591, 3-year: 0.666, 5-year: 0.633). The integrated signature outperformed these (1-year: 0.861, 3-year: 0.749, 5-year: 0.750), indicating synergistic predictive power. Compared to published signatures, our model showed superior 1-year AUC (0.861) versus a four-gene glycolysis signature (AUC: 0.740–0.806) (11), a six-gene signature (AUC: 0.719) (12), and an 11-gene signature (1-year: 0.719, 3-year: 0.762, 5-year: 0.742) (13), with comparable 3- and 5-year performance. The C-index for our signature (0.78) was also higher than the six-gene model’s (0.764) (12), underscoring its robustness in a multi-cohort validation.

Association of the CCRG-based SuperPC signature with clinicopathological features

To elucidate the relationship between the CCRG-derived risk signature and clinicopathological characteristics, we performed multiple subgroup analyses. A scatter plot showed a significant positive correlation between patient age and RS (r=0.085, P=0.049, Figure 3A). Significant differences in RS were observed across subgroups of pathological stage, T stage, PAM50 subtype, estrogen receptor (ER) status, progesterone receptor (PR) status, and HER2 status. Higher RS were associated with advanced pathological and T stages, the HER2-enriched PAM50 subtype, ER-negative and PR-negative status, and HER2-positive status (Figure 3B-3I).

Figure 3 Correlation of the CCRG-derived risk signature with clinicopathological features. (A) Scatter plot showing the correlation between patient age and risk score. Comparison of risk scores among subgroups based on (B) pathological stage, (C) T stage, (D) N stage, (E) radiotherapy status, (F) PAM50 subtype, (G) ER status, (H) PR status, and (I) HER2 status. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. CRG, Cori cycle-related genes; ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; N, node; PAM50, Prediction Analysis of Microarray 50; PR, progesterone receptor; T, tumor.

The CCRG-derived risk signature predicts immunotherapy response

To investigate the relationship between CCRGs and the tumor immune microenvironment (TIME) in breast cancer, we employed multiple algorithms for immune cell infiltration estimation. As shown in Figure 4A, the RS was significantly negatively correlated with the infiltration of CD4⁺ T cells, natural killer cells, and M2 macrophages, but positively correlated with M1 macrophages and Treg cells (Table S2). For instance, using CIBERSORT estimates in the TCGA-BRCA cohort, RS showed negative correlations with CD4⁺ T cells (r=−0.185, P<0.001), natural killer cells (r=−0.15, P=0.008), and M2 macrophages (r=−0.16, P<0.001), but positive correlations with M1 macrophages (r=0.19, P<0.001) and Treg cells (r=0.14, P<0.001). Furthermore, RS was negatively correlated with stromal score and IPS. The trends for PGK1 and PGAM1 were similar to RS, whereas HK1 often showed opposite correlations (Figure 4A). We also assessed correlations between CCRGs/RS and immune checkpoint/activity-related genes (Figure 4B, Table S3). RS, PGK1, and PGAM1 were predominantly negatively correlated with TBX2 but positively correlated with most other genes. Specifically, RS was positively correlated with PD-L1 (CD274; r=0.10, P=0.03), CTLA4 (r=0.14, P=0.003), and IDO1 (r=0.20, P<0.001), but negatively with TBX2 (r=−0.22, P<0.001). PGK1 and PGAM1 followed similar patterns (e.g., PGK1 with PD-L1: r=0.11, P=0.01; PGAM1 with CTLA4: r=0.15, P<0.001), while HK1 showed inverse trends (e.g., with CTLA4: r=−0.20, P<0.001). Conversely, HK1 was largely negatively correlated with most genes but positively correlated with TBX2. Additionally, in the IMvigor210 cohort treated with immunotherapy, high-risk patients showed significantly worse survival (P=0.045, Figure 4C).

Figure 4 Association of CCRGs and the derived signature with the immune landscape and therapy response in breast cancer. (A) Heatmap showing correlations between CCRGs/RS and parameters of the breast cancer immune landscape. (B) Heatmap showing correlations between CCRGs/RS and immune checkpoint/activity-related genes in breast cancer. (C) Kaplan-Meier survival analysis for high- and low-risk groups in the IMvigor210 immunotherapy cohort. *, P<0.05; **, P<0.01; ***, P<0.001. AZ, atezolizumab; CCRG, Cori cycle-related gene; CP, cytotoxic potential; EC, endothelial cell; ICG, immune checkpoint gene; IPS, immunoscore; MCP, monocyte chemoattractant protein; MHC, major histocompatibility complex; NK, natural killer; RS, risk score; SC, stem cell.

The CCRG-derived risk signature is associated with drug sensitivity

To further explore the relationship between the CCRG-derived risk signature and drug sensitivity, we evaluated its association with nine commonly used therapeutic agents. Correlation analysis (Figure 5) revealed that RS was significantly negatively correlated with sensitivity to docetaxel (r=−0.2, P=1.2e−05), paclitaxel (r=−0.16, P<0.001), and cisplatin (r=−0.12, P=0.01), but positively correlated with sensitivity to lapatinib (r=−0.11, P=0.02), temsirolimus (r=−0.24, P=1.3e−07), and metformin (r=−0.14, P=0.003).

Figure 5 Correlation analysis between the CCRG-derived risk score and sensitivity to common breast cancer therapeutic drugs. CCRG, Cori cycle-related gene.

BPs and pathways associated with the CCRG-derived signature

Enrichment analyses were performed to identify BPs and pathways linked to the CCRG-derived signature. Differential expression analysis identified genes significantly dysregulated between high- and low-risk groups (Figure 6A). GO enrichment analysis indicated these genes were involved in processes such as regulation of membrane potential, chromatid segregation, amino acid transport, and cell cycle regulation (Figure 6B). KEGG pathway analysis associated them with cell cycle, PPAR signaling pathway, neuroactive ligand-receptor interaction, IL-17 signaling pathway, and PI3K-AKT signaling pathway (Figure 6C). GSEA revealed that hallmark gene sets including glycolysis, unfolded protein response, estrogen response, epithelial-mesenchymal transition (EMT), and mitotic spindle were enriched in the high-risk group, while interferon gamma response, myogenesis, and myelocytomatosis (MYC) targets were enriched in the low-risk group (Figure 6D).

Figure 6 Biological processes and pathways associated with the CCRG-derived risk signature. (A) Volcano plot of DEGs between high- and low-risk groups in the TCGA-BRCA cohort. (B) GO and (C) KEGG enrichment analyses of DEGs between high- and low-risk groups. (D) Hallmark gene sets associated with the risk score. BP, biological process; BRCA, breast invasive carcinoma; CC, cellular component; CCRG, Cori cycle-related gene; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; NES, normalized enrichment score; NS, not significant; TCGA, The Cancer Genome Atlas.

Somatic mutation landscape across risk subgroups

Analysis of somatic mutations revealed several genes with significantly different mutation frequencies between the high- and low-risk groups, including tumor protein P53 (TP53), reelin (RELN), cadherin 1 (CDH1), and phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA) (Figure 7A). Most genes, except CDH1 and PIK3CA, had higher mutation frequencies in the high-risk group. Spearman correlation analysis showed a significant positive correlation between RS and TMB (r=0.389, P=7.6e−18, Figure 7B). Consistently, the high-risk group exhibited a significantly higher TMB than the low-risk group (Figure 7C), indicating an association between the CCRG-derived risk signature and somatic mutations.

Figure 7 Evaluation of the association between the risk signature and somatic mutations. (A) Genes with significantly different mutation frequencies between high- and low-risk groups. (B) Scatter plot of the correlation between RS and TMB. (C) Comparison of TMB between high- and low-risk groups. ***, P<0.001. RS, risk score; TMB, tumor mutational burden.

Construction of an RS-based nomogram for breast cancer

To further analyze the prognostic value of RS in conjunction with clinicopathological features, we performed Cox regression analyses. Univariate Cox analysis (Table 1) identified RS (HR =9.9, P<0.0001), age (HR =1.1, P<0.001), pathological N stage (HR =1.5, P=0.03), and pathological stage (HR =2.3, P=0.001) as risk factors, while radiotherapy (HR =0.45, P=0.02) was a protective factor. Multivariate Cox analysis confirmed RS (HR =14.80, P<0.001) and age (HR =1.04, P=0.002) as independent risk factors, and radiotherapy (HR =0.38, P=0.02) as an independent protective factor. Consequently, a nomogram was constructed based on RS, age, and radiotherapy status (Figure 8A). ROC curves showed that the nomogram achieved AUC values of 0.868, 0.773, and 0.763 for predicting 1-, 3-, and 5-year OS, respectively (Figure 8B), outperforming RS alone. Calibration curves demonstrated good agreement between the nomogram-predicted and observed 1-, 3-, and 5-year OS probabilities (Figure 8C). DCA indicated that the nomogram provided a higher standardized net benefit compared to individual clinicopathological features (Figure 8D).

Figure 8 Construction and validation of a breast cancer nomogram. (A) Nomogram built based on RS, age, and radiotherapy status selected from multivariate Cox analysis. (B) ROC curves evaluating the accuracy of the nomogram in predicting 1-, 3-, and 5-year OS. (C) Calibration curves for 1-, 3-, and 5-year OS, comparing predicted probabilities from the nomogram against actual observations. (D) Decision curve analysis comparing the standardized net benefit of the nomogram against other clinicopathological features. AUC, area under the curve; N, node; OS, overall survival; ROC, receiver operating characteristic; RS, risk score; T, tumor.

Pan-cancer analysis of CCRGs

To further evaluate the impact of CCRGs at a pan-cancer level, we calculated gene set variation analysis (GSVA) scores for the CCRG gene set and performed Cox and correlation analyses. CCRG scores were significantly higher in most cancer tissues compared to adjacent normal tissues, except in head and neck squamous cell carcinoma (HNSC) (Figure 9A). Cox analysis revealed that the CCRG score was significantly associated with prognosis in multiple cancers (Figure 9B), including BRCA (HR =1.99, P=0.01), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) (HR =2.797, P=0.007), HNSC (HR =2.664, P<0.001), kidney chromophobe (KICH) (HR =18.721, P=0.04), liver hepatocellular carcinoma (LIHC) (HR =2.275, P=0.02), lung adenocarcinoma (LUAD) (HR =2.486, P=0.001), mesothelioma (MESO) (HR =2.26, P=0.03), pancreatic adenocarcinoma (PAAD) (HR =2.961, P=0.002), and sarcoma (SARC) (HR =2.063, P=0.02). Furthermore, the CCRG score showed significant correlations with TMB and MSI across various cancers (Figure 9C,9D). These results suggest that CCRG activity is upregulated in cancers, plays a pro-tumorigenic role, and adversely affects prognosis.

Figure 9 Pan-cancer analysis of CCRGs. (A) Comparison of CCRG scores between tumor and normal tissues across cancer types. (B) Forest plot of HRs from Cox analysis of the association between CCRG score and prognosis in various cancers. (C) Pearson correlation between CCRG score and TMB across cancers. (D) Pearson correlation between CCRG score and MSI across cancers. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. CI, confidence interval; CCRG, Cori cycle-related gene; HR, hazard ratio; MSI, microsatellite instability; TMB, tumor mutational burden.

Single-cell analysis of CCRGs

To clarify the expression of RS-related CCRGs (HK1, PGK1, PGAM1) at the single-cell level, we analyzed several breast cancer spatial transcriptomics datasets (STDS0000027, STDS0000038, STDS0000049). A common observation was the relatively higher expression of PGK1 and PGAM1 compared to HK1, with widespread expression across most cell types (Figure 10A). Furthermore, we investigated the functional associations of these genes using four breast cancer single-cell transcriptomics datasets (CancerSEA ID: EXP0052, EXP0053, EXP0054, EXP0055). The analysis indicated that HK1 was negatively correlated with processes like DNA repair, DNA damage, cell cycle, and invasion. PGK1 was positively correlated with EMT, metastasis, stemness, and invasion. PGAM1 was positively correlated with proliferation, DNA repair, and DNA damage, but negatively correlated with angiogenesis (Figure 10B).

Figure 10 Expression and functional analysis of CCRGs at the single-cell level. (A) Expression of HK1, PGK1, and PGAM1 in breast cancer spatial transcriptomics data from the STOMICS database. (B) Functional correlation analysis of the genes based on the CancerSEA database. CCRG, Cori cycle-related gene; EMT, epithelial-mesenchymal transition.

Discussion

Our study developed a robust metabolic gene-based prognostic signature for breast cancer derived from the CCRG set and comprehensively validated its clinical utility, molecular relevance, and predictive power across multiple dimensions. The risk model, constructed using a machine learning-integrated approach, effectively stratified patients into high- and low-risk groups with distinct survival outcomes across independent cohorts, underscoring its reproducibility and generalizability.

Notably, the three-gene signature comprising HK1, PGK1, and PGAM1 reflects key aspects of glycolytic flux and gluconeogenesis—processes central to the Cori cycle. Elevated expression of these genes was significantly associated with poor prognosis, suggesting that dysregulated lactate recycling and glucose metabolism may fuel aggressive tumor behavior. This is consistent with previous studies highlighting the role of metabolic reprogramming in promoting proliferation, invasion, and treatment resistance in breast cancer. In particular, HK1 has been implicated in multiple cancer types as a promoter of tumor growth and metastasis. For instance, HK1 is overexpressed in colorectal cancer and is an independent predictor of poor survival (14). Similarly, in gastric cancer, high HK1 expression correlates with lymphatic metastasis and advanced TNM stage (15). In ovarian cancer, HK1 upregulation is associated with poor prognosis and promotes cell proliferation, invasion, and glycolysis (16). These findings align with our results, supporting the notion that HK1 is not only a metabolic enzyme but also a key oncogenic driver.

The uniqueness of our signature lies in its targeted focus on Cori cycle genes, distinguishing it from broader glycolysis-related models that include 4–11 genes (e.g., PGK1-inclusive signatures) but overlook the specific lactate-glucose recycling axis (11-13). By integrating machine learning (SuperPC), our parsimonious model achieves enhanced prognostic accuracy (e.g., higher 1-year AUC and C-index) compared to individual genes or these larger signatures, as demonstrated in our comparative analysis. This added value extends to its associations with immunosuppressive microenvironments and altered drug sensitivities—features that provide novel insights into metabolic-immune crosstalk and personalized therapy, beyond the survival prediction offered by prior models (11-13).

Specifically, PGAM1 has emerged as a critical regulator of glycolysis and tumor progression across multiple cancers. In breast cancer, PGAM1 promotes immunosuppressive M2 macrophage polarization through C-C motif chemokine ligand 2-mediated recruitment and JAK-STAT pathway activation, thereby remodeling the TME toward an immunosuppressive state (17). This aligns with our findings that the RS negatively correlates with stromal score and antitumor immune cells while positively correlating with immunosuppressive elements. Furthermore, PGAM1 inhibition has been shown to enhance ferroptosis and synergize with anti-programmed cell death 1 (PD-1) immunotherapy in hepatocellular carcinoma (18), suggesting a broader role in immune evasion and response to immunotherapy. PGAM1 is subject to complex regulatory mechanisms, including succinylation, phosphorylation, and protein-protein interactions, which modulate its enzymatic activity and non-metabolic functions. For instance, pyruvate kinase M2 can phosphorylate PGAM1 at H11 to enhance glycolysis shunts in cancer (19), while aspirin modulates PGAM1 succinylation via NF-κB/HAT1 signaling to restrict glycolysis in liver cancer (20). These regulatory layers highlight PGAM1 as a dynamic node in metabolic adaptation and tumor progression.

PGK1, as a crucial glycolytic enzyme, has been extensively reported to facilitate tumor progression across various cancers. In breast cancer, PGK1 is upregulated and promotes glycolysis, proliferation, and metastasis (21). Hypoxia-induced PGK1 expression enhances glycolysis and supports tumor survival under low-oxygen conditions (22). Furthermore, PGK1 is regulated by post-translational modifications such as succinylation and acetylation, which modulate its activity and stability, thereby influencing metabolic reprogramming and tumor growth (23,24). In TNBC, 3-oxoacid CoA-transferase 1-mediated succinylation of PGK1 at K146 stabilizes the protein and promotes aerobic glycolysis and immune escape by upregulating PD-L1 expression (25). Additionally, in luminal A breast cancer cells, acetylation of PGK1 at K323 enhances its enzymatic activity, glycolysis, and metastatic potential (24). These findings underscore the multifaceted role of PGK1 in driving breast cancer progression and immune evasion.

A particularly compelling finding was the strong association between the RS and immune infiltration patterns. The negative correlation with stromal score, IPS, and antitumor immune cells (e.g., CD4⁺ T and natural killer cells), along with a positive correlation with immunosuppressive elements (Tregs, M1 macrophages), implies that CCRG-driven metabolic remodeling may foster an immunosuppressive microenvironment. Quantitatively, this is evidenced by reduced infiltration of effector cells and elevated expression of checkpoints like PD-L1 in high-risk cases, aligning with reports of PGK1/PGAM1 promoting immune evasion via M2 polarization and PD-L1 upregulation (17,25,26). While our analyses reveal correlations, existing literature supports potential causal mechanisms linking high CCRG expression to immunosuppression. One hypothesis posits that tumor cells with elevated glycolysis (driven by HK1, PGK1, and PGAM1) actively shape the immunosuppressive environment through secretion of lactate, a Cori cycle byproduct abundant in the Warburg effect. Lactate acidifies the TME, directly inhibiting effector functions of CD4⁺ T cells and natural killer cells (e.g., via suppression of IFN-γ and granzyme B) while promoting polarization of macrophages toward an immunosuppressive M2 phenotype and expansion of Tregs, as demonstrated in breast cancer models (27,28). Alternatively, shared upstream drivers such as hypoxia-inducible factor-1α (HIF-1α) may orchestrate both phenomena: HIF-1α transcriptionally activates CCRGs to enhance glycolytic flux and lactate production, while simultaneously upregulating immunosuppressive checkpoints like PD-L1 and skewing immune cell metabolism (e.g., Treg switch from glycolysis to oxidative phosphorylation for enhanced suppression), fostering immune evasion (27,29). These hypotheses warrant experimental validation, such as through lactate-neutralizing agents or HIF-1α inhibitors in preclinical models.

Regarding the potential relationship between the Cori cycle and the breast cancer immune microenvironment, the cycle’s activation in tumors exacerbates lactate accumulation, which not only fuels cancer cell metabolism but also creates an acidic TME that impairs immune surveillance. For example, Cori cycle-driven lactate export from tumor cells can recruit myeloid-derived suppressor cells and inhibit dendritic cell maturation, further dampening adaptive immunity (28). In breast cancer, this interplay may explain subtype-specific immune coldness, such as in TNBC, where high glycolytic activity correlates with reduced T-cell infiltration (27). Shared regulators like HIF-1α reinforce this link by co-activating CCRGs and immunosuppressive pathways, suggesting a bidirectional metabolic-immune axis that our signature captures (27,29). This could partly explain the limited efficacy of immune checkpoint inhibitors (ICIs) in breast cancer and suggests that targeting metabolic pathways might reverse immune evasion. Interestingly, HK1 has also been shown to contribute to immune evasion in other contexts. For example, in hepatocellular carcinoma, HK1-containing extracellular vesicles derived from hepatic stellate cells promote tumor progression by enhancing glycolysis and modulating the TME (30). Similarly, PGK1 has been implicated in immune regulation. In HNSC, high PGK1 expression correlates with immune cell infiltration and poor prognosis (31). In glioblastoma, PGK1 succinylation regulated by the HIF1α/ATF3/P4HA1 axis promotes immune suppression and tumor growth (32). These studies highlight the role of PGK1 in shaping the TIME. Moreover, PGAM1 has been directly implicated in immune suppression. In colon cancer, PGAM1 inhibition reduces TAM infiltration and enhances CD8⁺ T cell recruitment, synergizing with anti-PD-1 therapy (26). Similarly, in melanoma, PGAM1 knockdown reduces immune checkpoint markers (e.g., PD-L1) and modulates EMT and apoptosis-related proteins (33), further underscoring its role in immune evasion and tumor progression.

To further elucidate the relationship between our CCRG signature and breast cancer treatment schemes, as well as the role of these genes in tumor-immune cell interactions, we note that high-RS correlate with resistance to chemotherapeutics like taxanes and cisplatin but sensitivity to targeted agents such as lapatinib and temsirolimus, suggesting metabolic vulnerabilities exploitable in precision medicine. For instance, targeting PGK1 disrupts glycolysis and sensitizes breast cancer cells to ferroptosis inducers (e.g., via interaction with glutathione peroxidase 4), potentially synergizing with immunotherapy by reversing immune escape (34). PGAM1 suppression remodels the TME in TNBC by reducing M2 macrophage infiltration and enhancing T-cell recruitment, improving responses to ICIs (35). HK1, through its interaction with VDAC, promotes tumor cell apoptosis and may augment chemotherapy efficacy when inhibited (36). Moreover, CCRGs mediate tumor-immune crosstalk: PGK1 in macrophages promotes aerobic glycolysis in tumor cells, fostering an immunosuppressive niche (37), while breast cancer cells reprogram natural killer cells to aid metastasis via metabolic signals (23). These insights position our signature as a biomarker for selecting patients for metabolic-targeted therapies (e.g., PGK1 inhibitors) combined with immunotherapy, warranting clinical trials to validate these interactions.

To enhance the clinical translational potential of our signature, it can be integrated with existing clinical indicators through advanced scoring systems like nomograms, as demonstrated in our study (Figure 8A), which outperformed individual factors in AUC. This combined approach could build a more accurate prognostic system by weighting metabolic RS alongside clinicopathological variables via multivariate Cox models or machine learning ensembles, enabling refined risk stratification. For postoperative adjuvant treatment decisions, high-risk patients (e.g., elevated CCRG scores with advanced stage) might benefit from intensified regimens, such as adding metabolic inhibitors (e.g., PGK1/PGAM1-targeted drugs) to standard chemotherapy or endocrine therapy, while low-risk cases could avoid overtreatment. In TNBC or HER2⁺ subtypes, the signature’s link to immunosuppression suggests prioritizing immunotherapy combinations (e.g., anti-PD-L1 with glycolysis inhibitors). Translational steps include developing quantitative polymerase chain reaction (qPCR)-based assays for CCRG expression in routine biopsies, prospective cohort validation, and incorporation into clinical trials [e.g., National Clinical Trial (NCT) trials evaluating metabolic-immune therapies], potentially improving precision oncology and patient outcomes (38).

Furthermore, the risk model correlated with sensitivity to conventional chemotherapies and targeted agents. The negative association with taxanes and platinum agents indicates potential resistance mechanisms active in high-risk cases, possibly mediated through metabolic adaptation. Conversely, positive correlations with lapatinib and temsirolimus suggest vulnerabilities that could be therapeutically exploited, supporting the concept of metabolic subtype-specific treatment strategies. In line with this, studies have shown that HK1 inhibition sensitizes cancer cells to metabolic inhibitors and chemotherapy. For instance, HK1 loss enhances the efficacy of metformin in ovarian cancer (39), and in multiple myeloma, HK1-negative/HK2-positive cells are sensitive to HK2-targeted therapy combined with oxidative phosphorylation inhibitors (40). PGK1 inhibition has also been shown to enhance radiosensitivity in esophageal squamous cell carcinoma by increasing reactive oxygen species (ROS) production and inhibiting the Akt/mTOR pathway (41). Moreover, targeting PGK1 with small molecules such as masitinib suppresses hepatocellular carcinoma growth under hypoxic conditions (22). Notably, PGAM1 has been linked to chemoresistance in various cancers. In ovarian cancer, PGAM1 promotes paclitaxel resistance through enhanced pyruvic acid production (42). Inhibition of PGAM1 sensitizes cancer cells to metabolic and chemotherapeutic agents, highlighting its potential as a therapeutic target to overcome drug resistance (18,26).

Enrichment analyses revealed that high-risk tumors are characterized by activation of cell cycle, PI3K-AKT, and glycolysis pathways, along with suppression of immune and interferon response pathways—a molecular landscape typical of aggressive, poorly immunogenic carcinomas. The significant correlation with TMB and specific mutational patterns (e.g., TP53, PIK3CA) further supports the interplay between metabolic dysregulation and genomic instability. Notably, HK1 is regulated by oncogenic signaling pathways such as AKT, which directly phosphorylates HK1 to enhance glycolysis and tumor growth (43). Additionally, KRAS4A directly binds and regulates HK1, highlighting a mechanism by which oncogenic mutations can directly modulate metabolic activity (44). PGK1 is also regulated by oncogenic pathways. In esophageal cancer, hypoxia-induced PGK1 expression promotes progression via the MYH9/GSK3β/β-catenin pathway (45). In colorectal cancer, PRMT1-mediated arginine methylation of PGK1 at R206 enhances its phosphorylation and promotes glycolysis and tumorigenesis (46). These studies illustrate the complex regulation of PGK1 by oncogenic signals and post-translational modifications. Meanwhile, PGAM1’s role in these processes is further supported by studies showing its involvement in Wnt/β-catenin signaling in breast cancer (47) and in regulating ASS1 expression via the cAMP/AMPK/CEBPB pathway (48), linking glycolysis to arginine metabolism and tumor growth.

Pan-cancer analysis revealed that CCRG dysregulation is a ubiquitous phenomenon across cancer types, showing consistent overexpression in tumors and association with poor prognosis, TMB, and MSI. This implies a conserved oncogenic role for Cori cycle genes beyond breast cancer, warranting further investigation into their utility as pan-cancer metabolic biomarkers. Specifically, HK1 is overexpressed in a variety of cancers and is linked to adverse outcomes, making it a broad-spectrum metabolic marker (14-16). PGK1 is also overexpressed in multiple cancers, including breast cancer (21), gastric cancer (49), and hepatocellular carcinoma (22), and is associated with poor prognosis and therapeutic resistance.

At single-cell resolution, we observed distinct expression patterns and functional associations of HK1, PGK1, and PGAM1. While PGK1 and PGAM1 were linked to pro-tumorigenic processes such as invasion and proliferation, HK1 showed a more complex, context-dependent role. These findings emphasize the importance of understanding metabolic gene function at the cellular and microenvironmental level. For example, HK1’s role can vary depending on cellular context: it promotes glycolysis and growth in some cancers (50), yet in others, its loss may paradoxically enhance malignancy but increase sensitivity to glycolytic inhibition (51). PGK1 also exhibits context-dependent functions. In liver cancer, mitochondrial import of PGK1 mediated by mcPGK1 promotes metabolic reprogramming and self-renewal of tumor-initiating cells (52). In osteosarcoma, PGK1 inhibition by icariside II suppresses EMT and metastasis (53). These studies highlight the diverse roles of PGK1 in different cancer types and cellular contexts.

Several limitations should be acknowledged. First, the study relied on retrospective public data; prospective validation is required before clinical translation. Second, although bioinformatic analyses suggest mechanistic links, experimental studies are needed to confirm causal relationships between CCRG expression and malignant phenotypes. Finally, the precise biological role of each gene within the Cori cycle and its interaction with immune cells remains to be fully elucidated.

Conclusions

In conclusion, we have developed and validated a CCRG signature that robustly predicts prognosis, immune contexture, and drug response in breast cancer. This signature provides novel insights into metabolic immunosuppression and offers an actionable tool for risk stratification and tailored therapy. Future work should focus on functional validation and integration into clinical trial design, particularly exploring combination therapies targeting HK1, PGK1, and related metabolic pathways.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2013/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2013/prf

Funding: This study was supported by Medical and Health Science and Technology Project of Zhejiang Province (Nos. 2023KY505, 2023KY046 and 2025HY0046) and Traditional Chinese Medicine Science and technology Project of Zhejiang Province (No. 2023ZL273).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-2013/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Zheng R, Zhang S, Zeng H, et al. Cancer incidence and mortality in China, 2016. J Natl Cancer Cent 2022;2:1-9. [Crossref] [PubMed]
2. Pan JW, Zabidi MMA, Ng PS, et al. The molecular landscape of Asian breast cancers reveals clinically relevant population-specific differences. Nat Commun 2020;11:6433. [Crossref] [PubMed]
3. Wang X, Zhao L, Song X, et al. Genomic and transcriptomic analyses identify distinctive features of triple-negative inflammatory breast cancer. NPJ Precis Oncol 2024;8:265. [Crossref] [PubMed]
4. De Placido P, Di Rienzo R, Pietroluongo E, et al. Insights on the association of anthropometric and metabolic variables with tumor features and genomic risk in luminal early breast cancer: Results of a multicentric prospective study. Eur J Cancer 2025;221:115409. [Crossref] [PubMed]
5. Aden D, Sureka N, Zaheer S, et al. Metabolic Reprogramming in Cancer: Implications for Immunosuppressive Microenvironment. Immunology 2025;174:30-72. [Crossref] [PubMed]
6. Bartoloni B, Mannelli M, Gamberi T, et al. The Multiple Roles of Lactate in the Skeletal Muscle. Cells 2024;13:1177. [Crossref] [PubMed]
7. Passarella S, Schurr A. l-Lactate Transport and Metabolism in Mitochondria of Hep G2 Cells-The Cori Cycle Revisited. Front Oncol 2018;8:120. [Crossref] [PubMed]
8. Yang M, Cai W, Lin Z, et al. Intermittent Hypoxia Promotes TAM-Induced Glycolysis in Laryngeal Cancer Cells via Regulation of HK1 Expression through Activation of ZBTB10. Int J Mol Sci 2023;24:14808. [Crossref] [PubMed]
9. Chen D, Liu P, Lu X, et al. Pan-cancer analysis implicates novel insights of lactate metabolism into immunotherapy response prediction and survival prognostication. J Exp Clin Cancer Res 2024;43:125. [Crossref] [PubMed]
10. Kulik U, Moesta C, Spanel R, et al. Dysfunctional Cori and Krebs cycle and inhibition of lactate transporters constitute a mechanism of primary nonfunction of fatty liver allografts. Transl Res 2024;264:33-65. [Crossref] [PubMed]
11. Zhang X, Wang J, Zhuang J, et al. A Novel Glycolysis-Related Four-mRNA Signature for Predicting the Survival of Patients With Breast Cancer. Front Genet 2021;12:606937. [Crossref] [PubMed]
12. He M, Hu C, Deng J, et al. Identification of a novel glycolysis-related signature to predict the prognosis of patients with breast cancer. World J Surg Oncol 2021;19:294. [Crossref] [PubMed]
13. Zhang D, Zheng Y, Yang S, et al. Identification of a Novel Glycolysis-Related Gene Signature for Predicting Breast Cancer Survival. Front Oncol 2020;10:596087. [Crossref] [PubMed]
14. He X, Lin X, Cai M, et al. Overexpression of Hexokinase 1 as a poor prognosticator in human colorectal cancer. Tumour Biol 2016;37:3887-95. [Crossref] [PubMed]
15. Gao Y, Xu D, Yu G, et al. Overexpression of metabolic markers HK1 and PKM2 contributes to lymphatic metastasis and adverse prognosis in Chinese gastric cancer. Int J Clin Exp Pathol 2015;8:9264-71.
16. Li Y, Tian H, Luo H, et al. Prognostic Significance and Related Mechanisms of Hexokinase 1 in Ovarian Cancer. Onco Targets Ther 2020;13:11583-94. [Crossref] [PubMed]
17. Zhang D, Wang M, Ma S, et al. Phosphoglycerate mutase 1 promotes breast cancer progression through inducing immunosuppressive M2 macrophages. Cancer Gene Ther 2024;31:1018-33. [Crossref] [PubMed]
18. Zheng Y, Wang Y, Lu Z, et al. PGAM1 Inhibition Promotes HCC Ferroptosis and Synergizes with Anti-PD-1 Immunotherapy. Adv Sci (Weinh) 2023;10:e2301928. [Crossref] [PubMed]
19. Wang Y, Shu H, Qu Y, et al. PKM2 functions as a histidine kinase to phosphorylate PGAM1 and increase glycolysis shunts in cancer. EMBO J 2024;43:2368-96. [Crossref] [PubMed]
20. Wang YF, Zhao LN, Geng Y, et al. Aspirin modulates succinylation of PGAM1K99 to restrict the glycolysis through NF-κB/HAT1/PGAM1 signaling in liver cancer. Acta Pharmacol Sin 2023;44:211-20. [Crossref] [PubMed]
21. Qiu A, Wen X, Zou Q, et al. Phosphoglycerate Kinase 1: An Effective Therapeutic Target in Cancer. Front Biosci (Landmark Ed) 2024;29:92. [Crossref] [PubMed]
22. Zhang Q, Zhang Y, Fu C, et al. CSTF2 Supports Hypoxia Tolerance in Hepatocellular Carcinoma by Enabling m6A Modification Evasion of PGK1 to Enhance Glycolysis. Cancer Res 2025;85:515-34. [Crossref] [PubMed]
23. Guo Z, Zhang Y, Wang H, et al. Hypoxia-induced downregulation of PGK1 crotonylation promotes tumorigenesis by coordinating glycolysis and the TCA cycle. Nat Commun 2024;15:6915. [Crossref] [PubMed]
24. Gao X, Pan T, Gao Y, et al. Acetylation of PGK1 at lysine 323 promotes glycolysis, cell proliferation, and metastasis in luminal A breast cancer cells. BMC Cancer 2024;24:1054. [Crossref] [PubMed]
25. Zhang H, Ling M, Zhang Y, et al. OXCT1 promotes triple negative breast cancer immune escape via modulating succinylation modification of PGK1. Commun Biol 2025;8:1033. [Crossref] [PubMed]
26. Wang C, Zhang M, Li S, et al. A phosphoglycerate mutase 1 allosteric inhibitor restrains TAM-mediated colon cancer progression. Acta Pharm Sin B 2024;14:4819-31. [Crossref] [PubMed]
27. Wang C, Xue L, Zhu W, et al. Lactate from glycolysis regulates inflammatory macrophage polarization in breast cancer. Cancer Immunol Immunother 2023;72:1917-32. [Crossref] [PubMed]
28. Zhao P, Wang S, Jiang J, et al. Targeting lactate metabolism and immune interaction in breast tumor via protease-triggered delivery. J Control Release 2023;358:706-17. [Crossref] [PubMed]
29. Kao TW, Bai GH, Wang TL, et al. Novel cancer treatment paradigm targeting hypoxia-induced factor in conjunction with current therapies to overcome resistance. J Exp Clin Cancer Res 2023;42:171. [Crossref] [PubMed]
30. Chen QT, Zhang ZY, Huang QL, et al. HK1 from hepatic stellate cell-derived extracellular vesicles promotes progression of hepatocellular carcinoma. Nat Metab 2022;4:1306-21. [Crossref] [PubMed]
31. Wang P, Wang YY, Xu YL, et al. Phosphoglycerate-kinase-1 Is a Potential Prognostic Biomarker in HNSCC and Correlates With Immune Cell Infiltration. Cancer Genomics Proteomics 2023;20:723-34. [Crossref] [PubMed]
32. Yang S, Zhan Q, Su D, et al. HIF1α/ATF3 partake in PGK1 K191/K192 succinylation by modulating P4HA1/succinate signaling in glioblastoma. Neuro Oncol 2024;26:1405-20. [Crossref] [PubMed]
33. Niu W, Yang Y, Teng Y, et al. Pan-Cancer Analysis of PGAM1 and Its Experimental Validation in Uveal Melanoma Progression. J Cancer 2024;15:2074-94. [Crossref] [PubMed]
34. He Y, Luo Y, Huang L, et al. Novel inhibitors targeting the PGK1 metabolic enzyme in glycolysis exhibit effective antitumor activity against kidney renal clear cell carcinoma in vitro and in vivo. Eur J Med Chem 2024;267:116209. [Crossref] [PubMed]
35. Zhang D, Wang M, Wang W, et al. PGAM1 suppression remodels the tumor microenvironment in triple-negative breast cancer and synergizes with anti-PD-1 immunotherapy. J Leukoc Biol 2024;116:579-88. [Crossref] [PubMed]
36. Ren YL, Li ZF, Chen K, et al. Chronic Dietary Exposure to Environmental Levels of Glyphosate Increases the Risk of Reproductive Dysfunction in Male Mice. Environ Sci Technol 2025;59:15705-19. [Crossref] [PubMed]
37. Zhang Y, Yu G, Chu H, et al. Macrophage-Associated PGK1 Phosphorylation Promotes Aerobic Glycolysis and Tumorigenesis. Mol Cell 2018;71:201-215.e7. [Crossref] [PubMed]
38. Gradishar WJ, Moran MS, Abraham J, et al. NCCN Guidelines® Insights: Breast Cancer, Version 4.2021. J Natl Compr Canc Netw 2021;19:484-93. [Crossref] [PubMed]
39. Šimčíková D, Gardáš D, Hložková K, et al. Loss of hexokinase 1 sensitizes ovarian cancer to high-dose metformin. Cancer Metab 2021;9:41. [Crossref] [PubMed]
40. Bischof H, Cisarova K, Burgstaller S, et al. Targeting hexokinase 2 to induce breast cancer cell senescence. Br J Pharmacol 2025; Epub ahead of print. [Crossref]
41. Chen J, Luo H, Wu X, et al. Inhibition of Phosphoglycerate Kinase 1 Enhances Radiosensitivity of Esophageal Squamous Cell Carcinoma to X-rays and Carbon Ion Irradiation. Front Biosci (Landmark Ed) 2025;30:36430. [Crossref] [PubMed]
42. Feng Y, Zhang X, Zhang S, et al. PGAM1 Promotes Glycolytic Metabolism and Paclitaxel Resistance via Pyruvic Acid Production in Ovarian Cancer Cells. Front Biosci (Landmark Ed) 2022;27:262. [Crossref] [PubMed]
43. Yu Y, Wang S, Wang Y, et al. AKT1 Promotes Tumorigenesis and Metastasis by Directly Phosphorylating Hexokinases. J Cell Biochem 2024;125:e30613. [Crossref] [PubMed]
44. Amendola CR, Mahaffey JP, Parker SJ, et al. KRAS4A directly regulates hexokinase 1. Nature 2019;576:482-6. [Crossref] [PubMed]
45. Xu JC, Wu LF, Chen TY, et al. Hypoxia-induced PGK1 expression promotes esophageal squamous cell carcinoma progression via stimulating MYH9-mediated GSK3β/β-catenin signalling. Clin Transl Med 2025;15:e70376. [Crossref] [PubMed]
46. Liu H, Chen X, Wang P, et al. PRMT1-mediated PGK1 arginine methylation promotes colorectal cancer glycolysis and tumorigenesis. Cell Death Dis 2024;15:170. [Crossref] [PubMed]
47. Wang Y, Liu W, Lai X, et al. PGAM1: a potential therapeutic target mediating Wnt/β-catenin signaling drives breast cancer progression. Discov Oncol 2025;16:161. [Crossref] [PubMed]
48. Liu M, Li R, Wang M, et al. PGAM1 regulation of ASS1 contributes to the progression of breast cancer through the cAMP/AMPK/CEBPB pathway. Mol Oncol 2022;16:2843-60. [Crossref] [PubMed]
49. Gu C, Xia Y, Lu C, et al. TRIM50 inhibits glycolysis and the malignant progression of gastric cancer by ubiquitinating PGK1. Int J Biol Sci 2024;20:3656-74. [Crossref] [PubMed]
50. Ni Y, Zhuang Z. DDX24 promotes tumor progression by mediating hexokinase-1 induced glycolysis in gastric cancer. Cell Signal 2024;114:110995. [Crossref] [PubMed]
51. Tseng PL, Chen CW, Hu KH, et al. The decrease of glycolytic enzyme hexokinase 1 accelerates tumor malignancy via deregulating energy metabolism but sensitizes cancer cells to 2-deoxyglucose inhibition. Oncotarget 2018;9:18949-69. [Crossref] [PubMed]
52. Chen Z, He Q, Lu T, et al. mcPGK1-dependent mitochondrial import of PGK1 promotes metabolic reprogramming and self-renewal of liver TICs. Nat Commun 2023;14:1121. [Crossref] [PubMed]
53. Hu J, Chen J, Zhao C, et al. Icariside II inhibits Epithelial-Mesenchymal transition in metastatic osteosarcoma by antagonizing the miR-194/215 cluster via PGK1. Biochem Pharmacol 2025;236:116838. [Crossref] [PubMed]