The public transcriptome expression matrices and detailed clinical data of BC were collected from The Cancer Genome Atlas (TCGA) database¹ (113 normal breast samples and 1,109 BC samples), the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) database² (1,904 BC samples), and GSE96058 in the Gene Expression Omnibus (GEO) database³ (3,409 BC samples), respectively. Besides, 1,082 BC patients of TCGA were selected as the training set, while 1,903 BC samples of METABRIC and 3,409 BC patients were chosen as external validation sets after excluding patients without overall survival information. In total, 284 lactate metabolism-related genes were acquired from the Molecular Signature Database v7.5.1 (MSigDB)⁴. Furthermore, 205 overlapping LMGs were screened out for further analyses after intersecting the above 284 LMGs with the total genes in TCGA-BRCA, METABRIC, and GSE96058 datasets (Supplementary Figure S1).

We obtained the information on the somatic mutation and copy number variation (CNV) of BC from TCGA-BRCA in the UCSC Xena database⁵. The somatic mutation frequency of 205 LMGs was analyzed by the R package “maftools” (20) and visualized in an oncoplot waterfall plot. Additionally, the CNV frequency of LMGs was calculated and visualized in a bi-directional column chart. The differentially expressed genes (DEGs) of LMGs were identified after comparing the normal breast samples and BC samples in TCGA-BRCA dataset with the threshold set to |log2FC| >1 and false-discovery rate (FDR) <0.05 using the “edgR” R package (21). A heatmap and a volcano plot of these significant DEGs were shown subsequently.

To elucidate the underlying prognostic significance of 205 LMGs in BC, overall survival (OS)-associated LMGs were identified through univariate Cox hazard regression analysis with p < 0.05 in the TCGA-BRCA (n = 1,082) and METABRIC (n = 1,903) cohorts, respectively (Supplementary Tables S1, S2). Moreover, the overlapping OS-associated LMGs were extracted to further research. Meanwhile, the expression levels and the location on chromosomes of those eligible LMGs were illustrated by the “RCircos” R package (22), and the correlation characteristics among these LMGs were demonstrated in a correlation matrix plot.

Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis was performed in the training cohort to construct the statistical prognostic signature utilizing the candidate OS-related LMGs. Next, 5 optimal LMGs were dug out to establish the prognostic model for BC patients, while the expression levels and prognostic significance of each signature-contained LMG were depicted, respectively. In according with the predictive prognostic signature, Lactate Metabolism Index (LMI) could be calculated for each BC patient by employing the following formula:

To make data and plots more intuitionistic, a linear transformation was carried out to adjust the LMI in each dataset using the following formula:

Afterwards the BC patients in each cohort could be separated into the high- and low-LMI groups by the cutoff of the median LMI value. Principal component analysis (PCA) was carried out to evaluate the classification accuracy of the signature. To discover the feasibility of the signature, a K-M analysis of OS was executed between the high- and low-LMI groups in three datasets separately.

To decipher the availability of LMI applied in actual clinical issues further, boxplots with the Kruskal test were exhibited to compare the distribution of adjusted LMI value in various degrees of diverse clinicopathologic parameters available in the three datasets. Besides, heatmaps were shown to unveil the relevance between each signature-included LMGs’ expression level and several clinical indicators, comprising of LMI, T stage, N stage, AJCC stage, PAM50 subtypes, and survival status in the training set in addition to LMI, tumor size, positive nodes, PAM50 subtypes, AJCC stage, and survival status in two validation datasets.

Subsequently, it was delineated whether the LMI was an independent prognostic predictor in BC by univariate and multivariate Cox regression analyses. Based on the results above, a lactate metabolism-correlated clinical nomogram that integrated LMI, age, and AJCC stage in TCGA-BRCA was established through the “rms” and “regplot” R packages (23). To assess the satisfactory predictive discrimination of nomogram, the calibration curve (24) and the decision curve analysis (DCA) plot were portrayed for BC patients.

“GSVA” R package (25) was employed to clarify the different biological functions and signaling pathways between high- and low-LMI groups in the training set. “c2.cp.kegg.v7.5.1.symbols.gmt” [KEGG] was retrieved from MSigDB as the reference molecular signature database and the pathways with adjusted p-values <0.05 were considered significant. Ultimately, the most significant pathways were displayed in a heatmap.

To verify the potential implications for immunotherapy based on LMI, the expression levels of several immunologic checkpoints, comprising PD-1, PD-L1, CTLA4, CD96, VSIR, and TIGIT, were compared between the high- and low-LMI groups using the Wilcox test.

To dissect the TME landscape between the two LMI subgroups, the ESTIMATE algorithm (26) was implemented to calculate the estimate scores, immune scores, and stromal scores for further predicting tumor purity and analyzing the TME. Moreover, the CIBERSORT deconvolution algorithm (27) was utilized to estimate the abundance of 22 tumor immune-infiltrating cell types in the training set.

Likewise, to investigate the correlation between LMI signature and cytokines in TME, several essential cytokines were picked out to make a comparison between high- and low LMI subgroups in expression levels, including IL-1B, IL-2, IL-6, IL-10, IL-18, TNF, IFNG, GZMA, and GZMB.

Human breast cancer cell lines were purchased from the American Type Culture Collection. All cell lines were cultured following standard guidelines. All cell lines were maintained without antibiotics in an atmosphere of 5% CO₂ and 99% relative humidity at 37°C. Cell lines were passaged for fewer than 6 months and were authenticated by short tandem repeat analysis. No mycoplasma infection was found for all cell lines.

Total RNA of cells was obtained with RNA-Quick Purification Kit (ES-RN001, Shanghai Yishan Biotechnology Co., Shanghai, China). The quantitative real-time PCR (qRT-PCR) plate was employed from NEST (402301, Wuxi NEST Biotechnology Co., Jiangsu, China). The primer sequences are shown in Supplementary Table S3. RNA levels were determined by qRT-PCR in triplicate on a Bio-Rad CFX96 using the SYBR Green method (RR420A, Takara, Mountain View, CA, USA). The RNA levels were normalized against β-actin RNA using the comparative Ct method.

All statistical analyses were completed via R software (Version 4.0.2, http://www.R-project.org). The discrepancy in the expression level of signature-encompassed LMGs in normal breast and BC samples, ESTIMATE algorithm-calculated scores, checkpoints, and cytokines in the low- and high-LMI groups were detected by the Wilcox test. The comparison of each Kaplan–Meier (KM) curve that occurred in this study was accomplished by the log-rank test. Also, Kruskal tests were performed to discover the differences in adjusted LMI values in various clinical parameters. Univariate and multivariate Cox regression analyses were employed to screen out the OS-related LMGs and the independent prognostic indicators of OS for BC. The correlation matrix plot was portrayed under Spearman’s correlation test. Statistical significance was confirmed as p-value <0.05, and all p-values were bilateral.

Initially, we assessed the global alterations of 205 LMGs in somatic mutation and copy number variation (CNV). As shown in Figure 1A, the top 10 genes with the highest somatic mutation rates were included in the heatmap, with the highest mutation frequency distributed in TP53. For frequency of CNV, the result showed that there were common CNV mutations among LMGs, and the top 20 genes in amplified and deleted CNV status separately were displayed in Figure 1B. Additionally, to investigate the differential expression of these LMGs, we compared their mRNA expression levels between 1082 tumor samples and 113 normal breast samples with the threshold of |log2FC| > 1 and FDR <0.05, and the results were demonstrated through a heatmap (Figure 1C) and a volcano plot (Figure 1D). Furthermore, to identify prognostic LMGs in breast cancer for the following research, the univariate Cox regression analyses were conducted to screen out OS-related genes in both TCGA and METABRIC datasets (Figure 1E). In total, 25 and 79 significant OS-related genes were obtained respectively and 9 overlapping genes (RPS14, MECP2, OCRL, RRM2B, GOT2, AIFM1, SDHA, SLC19A1, and CYC1) were included for further analysis after taking intersecting of the results above (Figure 1F). Moreover, the location of chromosomes and expression level of the 9 genes were illustrated by a Circos plot (Figure 1G). Finally, a correlation network plot was used to unravel the correlation features among 9 eligible FMGs (Figure 1H).

By performing the LASSO Cox regression analysis with 9 candidate genes in breast cancer patients of TCGA-BRCA training dataset, 5 pivotal genes were unearthed to establish the prognostic signature, Lactate Metabolism Index, namely LMI (Figures 2A, B), including RPS14, SLC19A1, CYC1, RRM2B, OCRL. Besides, the investigation of expression levels and survival capability of every signature-contained gene was further conducted by boxplots (Figure 2C) of mRNA expression levels and KM survival curves of OS (Figure 2D). From the results, we found that mRNA expressions of SLC19A1, CYC1, RRM2B, and OCRL were clearly elevated in breast cancer, while RPS14 expression was considerably decreased. Additionally, the RNA expression levels of SLC19A1 and RPS14 were validated in human breast cancer cell lines (Supplementary Figure S2), the result of which demonstrated that SLC19A1 was significantly promoted in breast cancer cell lines including MDA-MB-231, T47D, and SK-BR-3 while RPS14 declined significantly in BC cell lines comparing with breast epithelial cell line MCF10A. And for the survival analyses of OS, high expressions of SLC19A1, CYC1, RRM2B, and OCRL and downregulation of RPS14 were significantly related to more unfavorable OS in breast cancer, which further confirmed the validity of selected genes. Eventually, the prognostic signature was established as follows: LMI = Expression of SLC19A1 * 0.010382 − Expression of RPS14 * 0.000045 + Expression of CYC1 * 0.000408 + Expression of RRM2B * 0.000005 + Expression of OCRL * 0.006258.

To further verify the prognostic value of LMI in breast cancer, patients of the TCGA-BRCA training set and two validation sets (METABRIC and GSE96058) were then independently segregated into high-LMI and low-LMI subgroups according to the median value of LMI in each dataset (Figure 3A). As expected, we found that the number of deaths increased with LMI increasing in both the training and validation cohorts (Figure 3B). Meanwhile, PCA was applied to demonstrate the distribution patterns of the two subgroups in two-dimensional graphs (Figure 3C). Finally, the KM curves for the OS of breast cancer patients showed that patients with high LMI had considerably poorer survival probability across all cohorts (Figure 3D), which proved the predictive accuracy of LMI in the prognoses of breast cancer.

Further evaluation of the relationship between LMI and various clinicopathological factors in breast cancer patients was also conducted. In the TCGA-BRCA training cohort (Figure 4A), prominent discrepancies were found between LMI and clinical characteristics, including survival status, T, N, AJCC stage, and PAM50 subtypes. Similarly, conspicuous differences were also observed in validation cohorts, containing survival status, tumor size, positive nodes, stage, and PAM50 subtypes in the METABRIC set (Figure 4B) and survival status, tumor size, positive nodes, and PAM50 subtypes in the GSE96058 set (Figure 4C). Notably, as the LMI levels were higher, the degree of clinical features tended to become more severe, which was displayed in both the training and validation sets and reconfirmed the predictive value of LMI in breast cancer. Meanwhile, heatmaps were used to demonstrate the correlation analyses between LMI-contained genes and clinicopathological features in TCGA-BRCA (Figure 4D), METABRIC (Figure 4E), and GSE96058 (Figure 4F) sets.

Concerning whether LMI signature could be an independent prognostic predictor for breast cancer patients, univariate and multivariate Cox regression analyses were performed in the TCGA dataset. And the results illuminated that age, T stage, N stage, AJCC stage and LMI were significantly linked to OS in the univariate Cox analysis (Figure 5A), while only age, AJCC stage and LMI were still independent prognostic factors in the multivariate Cox analysis (Figure 5B). Next, a clinicopathologic nomogram, which contained LMI, age and AJCC stage, was developed to predict individual OS of 2-, 3- and 5-years (Figure 5C). Moreover, the calibration plot was portrayed to confirm the predictive consistency of the nomogram, which displayed great fitness (Figure 5D). And decision curve analysis (DCA) showed that the nomogram attained a greater clinical net benefit than any single clinical feature (Figure 5E). Collectively, the prognostic nomogram based on LMI signature has excellent prediction performance for OS of breast cancer patients.

To clarify the different biological functions and signaling pathways between high- and low-LMI groups in the training set, “GSVA” enrichment analysis was exploited to detect that there are substantial differences in the metabolic pathways and immune-related pathways in the two LMI subgroups of BC (Figure 6A). For instance, ether lipid metabolism, cytokine cytokine-receptor interaction, and arachidonic acid metabolism were primarily enriched in the high-LMI group, while selenoamino acid metabolism, pyrimidine metabolism, B-cell receptor signaling pathway, T-cell receptor signaling pathway, and natural killer cell-mediated cytotoxicity were considerably enriched in the low-LMI group.

As the aforementioned result unveiled the immune-related pathways with significant differences between the high- and low-LMI groups, considerable attention was paid to the TME landscape and characteristics of the two LMI subgroups. ESTIMATE algorithm was utilized to calculate the stromal score, immune score, and ESTIMATE score to evaluate and quantify the TME, the results of which unraveled that the low-LMI group of BC obtained the higher ESTIMATE score, stromal score, and immune score together with the diminished tumor purity than the high-LMI group (Figure 6B). Furthermore, the comparison of expression levels of the candidate immune checkpoints within the two LMI groups manifested that PD-1, CTLA4, CD96, VSIR, TIGIT except PD-L1 were significantly augmented in the low-LMI group (Figure 6C), which hinted that BC patients with lower LMI might acquire a more enhanced response to immunotherapy targeting the checkpoints above. Aiming to estimate the distribution of immune infiltrating cells in the TME of BC in various LMI groups, the CIBERSORT algorithm was performed to disclose that CD8+ T cells, gamma delta T cells, resting NK cells, resting dendritic cells, and neutrophils were notably enriched in TME of the low-LMI group while macrophages M0 and M2 were markedly strengthened in the high-LMI group (Figure 6D). Collectively, the results above shed light on that BC patients of the lower LMI could attain an immune-activated TME, while BC patients of the higher LMI might gain an immunosuppressive TME that contributed to the immune escape of tumor cells and a worse prognosis.

Additionally, it is well established that cytokines play a crucial role in the immune TME. Therefore, the boxplots in TCGA-BRCA dataset show that the expression levels of IL-2, IL-6, IL-18, TNF, IFNG, GZMA, and GZMB were promoted in the low-LMI group (Figure 7A). Likewise, in the METABRIC dataset, IL-1B, IL-6, IFNG, and GZMA were enhanced in the low-LMI group (Figure 7B). Meanwhile, IL-1B, IL-2, IL-6, IL-10, IL-18, IFNG, GZMA, and GZMB were elevated in BC patients of the low-LMI group (Figure 7C).