Prostate cancer (PCa) is one of the most frequently diagnosed solid malignancies in men and has become the fifth leading cause of male cancer death worldwide (1, 2). Due to the widespread use of prostate-specific antigen (PSA) serum test and the improvement of overall longevity, the incidence of this disease is increasing. Although radical prostatectomy (RP) leads to a favorable rate of cancer control, approximately 25% of all patients experience biochemical recurrence (BCR), which is determined by rising of serum PSA levels within 10 years of RP (3, 4). BCR indicates increased risk for metastasis and castration resistant PCa (5, 6). Exploring gene expressions that are closely correlated with BCR is of great importance. In this regard, more informative markers for assessing increased risk of BCR are highly needed.

Reprogramming of energy metabolism, especially abnormal activation of glycolysis (also known as the Warburg effect) in the presence of oxygen, has been recognized as one of the central hallmarks of cancer (7, 8). Cancer cells exhibit a higher level of glucose consumption and consequent lactate production (9). Glycolysis facilitates conversion of nutrient uptake into biomass and thereby sustains the rapid cancer cell growth (10). Moreover, activation of glycolysis coupled to increased production of lactic acid, promotes multiple cancer-promoting processes, including tumor invasion and metastasis, angiogenesis, suppression of the local anticancer immune response, as well as hypoxia resistance (11, 12). It has been reported that PCa cells develop high rates of glucose consumption in the metastatic stage (13, 14). PCa patients with highly glycolytic tumors have a worse prognosis (15). However, the association between the glycolysis-related gene signature and biochemical recurrence of PCa is largely unknown. We hypothesized that a signature based on glycolysis-related genes might be a concise and practical tool for predicting the BCR of PCa patients after RP. Here, we develop a glycolysis-based five-gene signature for BCR after RP using large-scale gene expression data obtained from The Cancer Genome Atlas (TCGA) and validate it using two Gene Expression Omnibus (GEO) dataset (GSE54460 and GSE70769).

All analyses were processed using the R software (v. 3.6.3).

The overall study procedures were summarized in Figure 1 . Univariate Cox proportional hazard model, robust likelihood-based survival model and multivariate Cox proportional hazards model were used to establish a five-gene signature that can predict the BCR free survival.

Considering that one of the major turning points in PCa development is the progression to BCR, it is highly desirable to effectively assess PCa patients with increased risk of BCR (20, 21). Accumulating researches have been made to identify biomarkers to improve the prediction of patients with BCR. It has been demonstrated that certain clinical parameters, including pathologic stage, Gleason score, lymphonode metastasis, and lymphovascular invasion, are associated with BCR (22–24). Moreover, several gene signatures have been established to predict BCR after prostatectomy. For instance, it has been shown that the cell cycle progression (CCP) score, an RNA expression signature based on the levels of 31 CCP genes, can predict BCR free survival (25). Jiang et al. extracted 696 differentially expressed genes from the TCGA dataset and developed a 15-gene signature (SigMuc1NW) using Elasticnet for prediction of BCR (26). In addition, signatures based on non-coding RNA, such as long non-coding RNAs (lncRNAs) and microRNAs (miRNAs), also have been established to predict BCR. Using lncRNAs that were differentially expressed between tumor and normal prostate tissues, as well as between high and low Gleason score tumor tissues, Shao et al. constructed a seven-lncRNA signature that can predict BCR (27). Moreover, five miRNAs (miR-30c-5p/31-5p/141-3p/148a-3p/miR-221-3p) were identified as independent prognostic biomarkers for BCR (28). However, the predictive value of glycolysis-related gene signature in BCR remains largely unknown. Given the important roles of elevated glycolysis in cancer development and progression, we would like to explore whether it is possible to establish a robust glycolysis-based gene signature to predict the BCR in PCa patients.

As compared with a single biomarker, integrating multiple biomarkers into an aggregated model with bioinformatics analysis would substantially improve the predictive performance (29, 30). Here, we performed a multistep analysis to identify a glycolysis-based gene signature which could predict BCR free survival in patients with PCa. Considering that differentially expressed genes between normal tissues and malignant specimens may not be associated with BCR at all, a univariable Cox analysis was carried out for the primary screening. Robust likelihood-based survival modeling, which selects predictive factors based on the partial likelihood of the Cox model, is commonly used in construction of predictive signature for cancers (31, 32). We used robust likelihood-based survival analysis and multiply Cox regression model to establish a five-gene signature with prediction value. Cut-off scores are often set arbitrarily and vary between reports. ROC curve analysis can be used as an alternative method in the selection and validation of cut-off scores for determining clinically relevant threshold (33, 34). We used ROC analysis to identify the optimal cut-off point for dividing patients into low risk and high risk groups. We found that the BCR free survival time of high risk group was significantly shortened. Importantly, results from the validation cohorts confirmed the robustness of the glycolysis-based five-gene signature, suggesting the excellent performance and consistency of this model throughout the training cohort and two validation cohorts. These data indicate that the five-gene signature exhibits a robust prediction value for BCR free survival in patients with PCa.

We identified five glycolysis-related mRNAs (GYS2, STMN1, PPFIA4, KDELR3 and ABCB6) which were associated with BCR after RP. Among the 5 genes, the expression levels of the genes with negative coefficients (GYS2 and KDELR3) were increased in BCR free patients. Moreover, STMN1, PPFIA4, and ABCB6 have positive coefficients, and their expression levels were upregulated in patients with BCR. It is of note that bulk tumor mass is composed of diverse cells, including malignant, stromal and immune cells (35). Moreover, the datasets in current study are based on the conventional bulk-level molecular profiling. The molecular profiles of these approaches represent an average readout from all cell types within the tissue. The averaging over the individual cells leads to information loss (36). Therefore, the change of the expression levels of the five genes may partially be due to the alterations in stromal amount. Recent developments in single-cell RNA sequencing (scRNA-seq) have enabled the transcriptomes of single cells to be sequenced in a high throughput manner (37). ScRNA-seq provides a comprehensive and precise analysis of the cancer cell genome (38). Further studies based on the new powerful approach in large cohort of PCa patients will provide new insights into the characteristics of the disease and facilitate the exploration of new markers and therapeutic targets.

Some of the 5 genes have been implicated to be involved in cancers, including PCa. GYS2 encodes a protein that catalyzes the rate-limiting step in the synthesis of glycogen (39). It has been reported that GYS2 acted as a tumor suppressor in hepatocellular carcinoma (HCC) (40). The research found that GYS2 inhibited the proliferation of HCC cells via a negative feedback loop with p53 (40). A series studies have demonstrated the oncogenic role of STMN1 in various kinds of cancers (41–43). STMN1 was overexpressed in PCa and its expression was associated with the malignant behavior of the disease (44). MiR-34a, a tumor suppressor miRNA, inhibited the progressive phenotypes of PCa cells via directly regulating STMN1 (45). PPFIA4 belongs to the liprin-alpha gene family and inhibition of PPFIA4 reduced pancreatic cancer cell proliferation and invasion (46). Moreover, suppression of PPFIA4 promoted chemosensitivity of small lung cancer (SCL) cells under hypoxia (47). KDELR3 is a gene which encodes a member of the KDEL endoplasmic reticulum protein retention receptor family. The research by Marie et al. showed that silencing of KDELR3 reduced lung colonization of melanoma cells in experimental metastasis assays via regulating the metastasis suppressor, KAI1 (48). In addition, overexpression of ABCB6, a member of the ATP-binding cassette (ABC) transporter superfamily, could enhance the accumulation of protoporphyrin IX and improve the efficacy of 5-aminolevulinic acid-based photodynamic therapy in glioma (49). Karatas et al. reported that the expression of ABCB6 was elevated in PCa tissues as compared with that in normal prostate tissues. Consistently, they also found that ABCB6 was overexpressed in recurrent PCa in comparison with non-recurrent PCa (50). However, the unrecognized roles of GYS2, PPFIA4 and KDELR3 in PCa are worth further investigating to identify the biological functions and underlying mechanisms of theirs in the development and progression of the disease. Further delineation of molecules from the signature will provide new insights into the etiology of the disease and might uncover potential therapeutic targets.

In conclusion, we conducted an integrated study to develop a glycolysis-based five-gene signature for the prediction of the BCR free survival of PCa patients after RP. Future prospective clinical trials are warranted to evaluate the clinical utilization of this signature.

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/, The Cancer Genome Atlas prostate adenocarcinoma dataset.

ZX, LX, and LL: data acquisition, analysis, and manuscript writing. HL, JJ, MP, and YH: interpretation of data and statistical analyses. HX, HG, and YL: study concept, design, and supervision. All authors contributed to the article and approved the submitted version.

The funding for this project was provided by the National Natural Science Foundation of China (No. 82073050 and 81772486); Guangdong Basic and Applied Basic Research Foundation (No. 2019A1515012046 and 2020A1515010049); Guangdong Medical Scientific Research Foundation (No. A2020343); and Young Teachers Cultivate Projects of Sun Yat-sen University (No. 20ykpy60).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.