Three bulk RNA sequencing datasets (GSE1456, GSE7390, GSE86166) and a single-cell sequencing dataset (GSE161529) were extracted from the Gene Expression Omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo/). Corresponding clinical information was obtained from the easyGEO website (https://tau.cmmt.ubc.ca/eVITTA/easyGEO/). GSE1456 (GPL96, 159 samples) [14] and GSE7390 (GPL96, 198 samples) [15] were combined and used as the training set, while GSE86166 (GPL15048, 366 samples) [16] was used to validate the risk model. Among which, 126 of 159 patients have received tamoxifen and its combination in GSE1456 cohort. Patients in GSE7390 and GSE86166 cohorts have not received neoadjuvant therapy. Data preprocessing was conducted as described previously [17, 18, 19], and R package “limma” (version 3.56.2) was used to remove batch effects. GRGs were obtained from the Reactome database (https://reactome.org/).

Dataset GSE161529 was downloaded from GEO and analyzed using the R package “Seurat” (version 5.0.3). Each gene was expressed in at least three cells. The number of expressed genes in cells was between 300 and 6000. Filters were applied with mitochondrial genes $>$20%, ribosomal genes 3%, and hemoglobin genes $>$0.1%. The gene expression matrix was normalized, and the top 2000 most variable genes were identified and used as input for subsequent principal component analysis (PCA). Data was clustered by the FindClusters function (resolution = 0.5). Finally, the cell clusters were annotated based on markers from CellMarker2.0 (https://bio.tools/cellmarker_2.0) [20] . Differential marker genes for cell types were identified by the “FindAllMarkers” function (logFC = 0.25, min.pct = 0.25, p 0.05). The R package “AUCell” (version 1.26.0) [21] was used to calculate GAS expression enrichment scores in single cells.

The R package “ConsensusClusterPlus” (version 1.66.0) was used to determine members in the dataset and the number of possible clusters. At each iteration, 80% of the samples are subsampled and each subsample is divided into k groups (maximum k = 6) by the k-means algorithm based on Euclidean distance. This process was repeated 1000 times. Subsequently, the cumulative distribution function (CDF) curve was used to identify the optimal number of clusters.

Univariate Cox analysis was used to identify the GAS and to calculate regression coefficients. Genes with p 0.05 were used to develop risk scores. Overall survival (OS) and relapse free survival (RFS) were evaluated using the R package “survminer” (version 0.4.9). The risk score was calculated based on the formula: risk score = coef (gene1) $\times$ expr (gene1) + coef (gene2) $\times$ expr (gene2) +…+ coef (geneN) $\times$ expr (geneN). Patients were divided into high- and low-risk groups according to the median risk score. The R package “timeROC” (version 0.4) was used to perform time-dependent receiver operating characteristic (ROC) curve analysis. This model was validated using the GSE86166 cohort.

Single-sample Gene Set Enrichment Analysis (ssGSEA) was used to evaluate the abundance of immune cell infiltrate in each group. Differences in the expression of immune-inhibitory and immune-stimulatory proteins between HGAS and low GAS (LGAS) was calculated and visualized using the R package “ggplot2” (version 3.5.0).

WGCNA was performed as described in a previous study [19]. The optimal soft power threshold ($\beta$ = 6) was then calculated to construct a scale-free network. The module characteristic genes and clinical traits were correlated, together with the gene significance (GS) of genes in the main module. Hub genes were identified based on GS $>$0.5 and module membership (MM) $>$0.8.

Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of hub genes were performed using the R package “clusterProfiler” (version 4.10.1). Analysis of protein-protein interaction (PPI) was performed using the GENEMANIA (https://genemania.org/) database. The “c2.cp.kegg.v2023.1.Hs.symbols.gmt” gene set was obtained from the Molecular Signature Database (https://www.gsea-msigdb.org/gsea/msigdb) and inter-group differences were analyzed. Gene set enrichment analysis was performed using the R package “GSVA” (version 1.50.1) [22].

The normal breast epithelial cell line MCF-10A and three human BC cell lines (MDA-MB-231, MCF-7, BT474) were obtained from Procell (Wuhan, Hubei, China). All cell lines were validated by STR profiling and tested negative for mycoplasma. The MCF-10A, MDA-MB-231 and MCF-7 cell lines were cultured as described previously by our group [18]. BT474 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Procell, Wuhan, China) supplemented with 10% fetal bovine serum (FBS; Procell, China), and 1% penicillin and streptomycin (Solarbio, Shanghai, China). All cells were kept at 37 °C in a humidified atmosphere containing 5% CO₂.

Total RNA was extracted from cell lines using TransZol reagent (TransGen, Beijing, China) according to the manufacturer’s instructions. The primer sequences were: CCNB2, CAACCCACCAAAACAACA (forward), AGAGCAAGGCATCAGAAA (reverse); GAPDH, GCACCGTCAAGGCTGAGAAC (forward), TGGTGAAGACGCCAGTGGA (reverse). Among which, the relative expression levels of mRNA were calculated by the 2^{-Δ⁢Δ⁢Ct} method, and GAPDH was used as an internal control.

Tissue samples from two triple-negative BC (TNBC) patients who did not undergo neoadjuvant therapy were collected from the First Affiliated Hospital of Nanchang University. Immunohistochemistry was performed using the mouse monoclonal anti-CCNB2 (1:500 dilution, Abcam, Cambridge, UK) antibody as described previously [23]. All experiments were performed in compliance with the relevant regulations. Ethics approval was obtained from the Medical Research Ethics Committee of the First Affiliated Hospital of Nanchang University [Grant no. (2023) CDYFYYLK (04-036)]. Patients provided written informed consent.

Drug response data was obtained from the Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org/). The R package “oncoPredict” (version 0.2) [24] was used to analyze drug sensitivity in HGAS and LGAS groups, which were compared using the Wilcoxon rank sum test.

All statistical analyses were performed using R language (version 4.3.1). Statistical differences between two groups were evaluated using the Wilcoxon test. p 0.05 was considered statistically significant, and all statistical tests were two-sided.

After removing batch effects, a total of 85 glycolytic genes were initially identified from a combined dataset (Supplementary Fig. 1A,B, Fig. 1A), of which 21 were then selected using univariate Cox analysis (Supplementary Fig. 1C, Supplementary Table 1). Consensus clustering analysis revealed that BC patients could be divided into two clusters based on the best parameter (k = 2) (Fig. 1B,C, Supplementary Fig. 1D). Heatmap analysis showed that 15 of the 21 glycolytic genes were upregulated in cluster 1, which was thus defined as HGAS. In addition, 4 of the 21 glycolytic genes were downregulated in cluster 2, which was defined as LGAS (Fig. 1D). BC patients in the HGAS group had significantly worse RFS and OS compared to the LGAS group (Fig. 1E). In addition, BC patients in the HGAS group showed higher tumor grade and a higher death and recurrence rate (Table 1).

Fig. 1.

Identification of HGAS in BC. (A) Venn diagram. (B) Consensus score matrix of all samples when k = 2. (C) The principal component analysis (PCA) map. (D) Heatmap of mRNA expression levels of glycolysis genes shows distinct clustering of samples into HGAS and LGAS, the red font represents upregulation and blue font represents downregulation. (E) Kaplan–Meier analysis showing OS and RFS of the two clusters of BC patients. *p 0.05, **p 0.01, ****p 0.0001; ns, not significant. HGAS, high glycolytic activity signature (GAS); LGAS, low GAS; BC, breast cancer; OS, overall survival; RFS, relapse free survival.

We next constructed a prognostic risk model using the coefficient (Supplementary Table 1) by univariate Cox analysis. This was used to assess GAS as a predictor of OS. The model showed that BC patients in the high-risk group had significantly worse OS than those in the low-risk group (Fig. 2A,B). The predictive performance of the risk score for OS was evaluated by time-dependent ROC curves. The area under the curve (AUC) was 0.78 at 1-year, 0.69 at 3-years, and 0.66 at 5-years (Fig. 2C). The nomogram revealed that the risk score contributed more to the prognosis than clinical features (Fig. 2D,E). Similar results were obtained using the validation dataset (Fig. 2F–J). Collectively, these results indicate that the prognostic risk model developed here was an excellent factor in predicting OS.

Fig. 2.

Construction and validation of a Glycolytic Activity Signature (GAS) risk model. Distributions of OS status, OS and risk score in the training set (A) and test set (F). Kaplan-Meier analysis of OS in the training set (B) and test set (G). Receiver operating characteristic (ROC) curve for 3-year, 5-year and 10-year survival in the training set (C) and test set (H). Nomogram developed by integrating the signature risk-score with the clinicopathologic features in the training set (D) and test set (I). Calibration curves of nomogram for predicting overall survival at 3-year, 5-year and 10-years in the training set (E) and test set (J). ID, identity document; OS, overall survival; AUC, area under the curve.

The immune microenvironment is known to significantly affect therapeutic effects and prognosis of BC. We therefore compared the immune response between HGAS and LGAS groups using several different immune analyses. The ESTIMATE analysis showed that BC in the HGAS group had a lower stromal score and a higher immune score compared to BC from the LGAS group (Supplementary Fig. 2). Single-sample gene set enrichment analysis (ssGSEA) results also showed the HGAS group had a relatively high immune status compared to the LGAS group (Fig. 3A), consistent with the ESTIMATE results. Eighteen of 22 immune cell types showed relatively higher infiltration in the HGAS group compared to LGAS group, including regulatory T cells, myeloid derived suppressor cells (MDSCs) and so on. Four of 22 immune cell types showed relatively lower infiltration in the HGAS group, including plasmacytoid dendritic cells, macrophages, eosinophils and mast cells.

Fig. 3.

Immune response comparison between the HGAS and LGAS BC groups. (A) Immune cell infiltration. (B) Immune modulator molecules. (C) Differences in scores for infiltrating immune cells. (D) Differences in scores for immune-related functions. (E) Correlation between glycolysis-related genes (GRGs) and immune cells in HGAS group. (F) Pathway enrichment analysis in the two groups. *p 0.05, **p 0.01, ***p 0.001 and ****p 0.0001. ns, no significance; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The HGAS and LGAS BC groups were also compared in terms of their immunomodulator expression level (antigen presentation, cell adhesion, co-inhibitors, co-stimulators, ligands, receptors, etc.). Most immunomodulators were expressed at higher levels in the HGAS group, with the exception of selectin P (SELP), nectin cell adhesion molecule 2 (PVRL2), C-X-C motif chemokine ligand 12 (CXCL12) and endothelin receptor type B (EDNRB) (Fig. 3B). In addition, ssGSEA enrichment scores were significantly greater in the HGAS group, with the exception of mast cells and Th2 cells (Fig. 3C). Similarly, the HGAS group showed enrichment for inflammatory conditions, checkpoints, cytolytic exertion, T cell co-inhibition, T cell co-stimulation, para-inflammation, and type 1-related IFN reactions. HGAS BC also showed higher immune function scores for major histocompatibility complex (MHC) class I, Antigen-presenting cells (APC) co-inhibition, APC co-stimulation, chemokine receptor (CCR), and human leukocyte antigen (HLA) (Fig. 3D). Thus, the correlations between GRGs and infiltrated immune cell types in HGAS group were analyzed (Fig. 3E). GRGs showed higher positively correlation with most immune cells, myeloid-derived suppressor cell (MDSC) in particularly (Fig. 3E).

Pathway enrichment analysis was used to reveal functional differences between the two groups (Fig. 3F). Almost all pathways were upregulated in HGAS BC compared to LGAS BC, including metabolism pathways such as the pentose phosphate pathway and methionine metabolism. Of note, the cell cycle pathway was significantly increased in the HGAS group (Fig. 3F), suggesting that it could be induced by high glycolytic activity in tumor cells.

To further explore the relationship between high glycolytic activity and the cell cycle, we first analyzed GAS-related genes by constructing co-expression modules using WGCNA. All candidate genes were clustered and divided into 16 modules according to the best soft threshold ($\beta$ = 4) with a cutoff R² value of 0.92 (Fig. 4A,B, Supplementary Fig. 3A–C). Tumor grade was used as a key trait to identify GAS-related genes (Fig. 4C). The blue module was found to be most closely related to GAS (r = –0.64, p = 4 $\times$ 10^-40) and tumor grade (r = 0.57, p = 8 $\times$ 10^-31), and was therefore identified as the best module for further analysis. Subsequently, 22 and 19 genes were identified based on an MM $>$0.8 and GS $>$0.5 (Fig. 4D–G), respectively. Ultimately, 17 intersected genes were identified as hub genes for additional analysis (Fig. 4H).

Fig. 4.

Construction of co-expression modules using WGCNA in GEO data. (A) Heatmap depicting the topological overlap matrix among genes based on co-expression modules. (B) Clustering dendrograms of genes with dissimilarity based on topological overlap, together with assigned module colors. (C) Relationships of modules with GAS and grade traits. (D,F) Correlations between blue modules and GAS and grade traits. (E,G) Candidate genes in the blue module related to both traits. (H) Hub genes. WGCNA, weighted correlation network analysis; GEO, Gene Expression Omnibus.

GO enrichment and pathway enrichment analysis were performed to explore the potential function of GAS-related genes. GO analysis indicated that various biological processes were enriched, especially cell cycle-related processes such as regulation of the metaphase/anaphase transition and the mitotic phase transition (Fig. 5A). For cellular component ontology, many genes were mapped into the chromosomal region spindle. For molecular function, the genes were primarily mapped into cyclin-dependent protein serine/threonine kinase regulator activity and ATP hydrolysis activity (Fig. 5A). Importantly, pathway enrichment analysis revealed significant enrichment for the cell cycle and cellular senescence (Fig. 5B). Functional analysis of GAS-related genes therefore highlighted the importance of the cell cycle, suggesting potential roles for glycolytic activity and the cell cycle in the regulation of BC progression. GAS-regulated genes were then mapped to the cell cycle (Fig. 5C). Four GAS-related genes were found to be associated with the cell cycle, with CCNB2 participating in the most pathways (Fig. 5B,C). Moreover, PPI results indicated that CCNB2 interacts with many of the key regulatory proteins in the cell cycle (Fig. 5D). Collectively, these findings suggest that CCNB2 is a core factor linking the cell cycle with glycolytic activity in BC.

Fig. 5.

Functional analysis and key gene of the cell cycle. (A) Gene Ontology (GO) analysis of hub genes. (B) Pathway enrichment analysis of hub genes. (C) Pathways of cell cycle-related hub genes. (D) Protein-protein interactions. BP, Biological process; CC, cellular component; MF, Molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

We next explored the functions of CCNB2 in BC. Elevated CCNB2 expression was observed in BC tissue compared to normal breast tissue in the The Cancer Genome Atlas (TCGA)-BRCA and GSE16228 cohorts (Fig. 6A,B). This high expression level was then validated in three BC cell lines (Fig. 6C). CCNB2 protein was also highly expressed in BC patients (Fig. 6D). In addition, the results showed that CCNB2 was a potential prognostic factor, with high expression being associated with significantly worse OS and RFS in BC patients (Fig. 6E,F).

Fig. 6.

Cyclin B2 (CCNB2) expression analysis and prognostic significance. CCNB2 expression in the TCGA cohort (A), in GSE16228 (B), and in BC lines (C); Immunohistochemistry analysis of CCNB2 expression in breast normal and tumor tissue based on the Human Protein Atlas (HPA). Images available from v23.0.proteinatlas.org (D); Kaplan–Meier analysis of OS (E) and RFS (F) in BC patients according to CCNB2 expression. *p 0.05, ****p 0.0001. BRCA, Breast invasive carcinoma; TCGA, The Cancer Genome Atlas.

The relationship between high glycolytic activity and CCNB2 expression was further investigated using scRNA-seq data. First, ineligible cells were filtered out to yield 46046 core cells (normal, 9376; Her2, 10934; luminalA/B, 10232; TNBC, 15504) for subsequent analysis of different tumor cell subtypes and normal samples (Supplementary Fig. 4). Fifteen cell clusters were finally identified (Fig. 7A, Supplementary Fig. 5B–D). High CCNB2 expression was observed in proliferating cells in the TNBC subtype (Fig. 7B). Importantly, proliferating cells were the primary cell type in TNBC samples (Fig. 7C). Increased expression of CCNB2 was also observed in TNBC patients (Fig. 7D). A high AUC value was also observed for the 21 HGAS-related genes in proliferating cells in the TNBC subtype (Fig. 7E). These results highlight the close link between elevated glycolytic activity and CCNB2 expression.

Fig. 7.

Single-cell clustering analysis. (A) Results of cell annotation. (B) Uniform Manifold Approximation and Projection (UMAP) plots showing the expression of CCNB2 in different cell clusters. (C) Cell types in different samples. (D) Immunohistochemistry for CCNB2 expression in normal and tumor tissue from triple-negative BC (TNBC) patients. Scale bar: 80 µm. (E) AUC value for glycolytic activity in different cell types. CAFs, cancer-associated fibroblasts; HER2, Human epidermalgrowth factor receptor-2; ER, estrogen receptor.

Since CCNB2 can potentially interact with a variety of cell cycle-related proteins, we further investigated the sensitivity of cell cycle-related clinical drugs in relation to GAS and CCNB2. In all, BC patients in HGAS group were more sensitive to 13 cell cycle-related clinical drugs due to the lower half-maximal inhibitory concentration (IC₅₀) score, including dinaciclib, ribociclib, camptothecin, vinblastine, cytarabine, docetaxel, fluorouracil, paclitaxel, irinotecan, gemcitabine, topotecan, teniposide and vinorelbine (Fig. 8A). Among them, BC patients in the high CCNB2 expression group were more sensitive to 9 cell cycle-related clinical drugs due to the lower IC₅₀ score, including dinacicilib, ribocicilib, vinblastine, docetaxel, fluorouracil, paclitaxel, gemcitabine, teniposide and vinorelbine (Fig. 8B). The drug score (IC₅₀) of each of the 9 cell cycle-related clinical drugs showed a significant negative correlation with CCNB2 expression (Fig. 8C). Collectively, these findings suggest that these 9 medications may be more effective in patients with higher HGAS levels and CCNB2 expression.

Fig. 8.

The IC₅₀ score of cell cycle drugs. Drug score (IC₅₀) in HGAS and LGAS BC groups (A), and in CCNB2 low- and high-level groups (B); (C) Correlation between CCNB2 expression and drug score (IC₅₀) to cell cycle drugs. The lower the drug score (IC₅₀, the number in Y-axis) is, the higher sensitivity of the patients to the drug under the condition of high-level GAS or CCNB2. *p 0.05, **p 0.01, ***p 0.001 and ****p 0.0001. IC₅₀, half-maximal inhibitory concentration. ns, no significance.