The effect of a novel glycolysis-related gene signature on progression, prognosis and immune microenvironment of renal cell carcinoma
BMC Cancer volume 20, Article number: 1207 (2020)
Glycolysis is a central metabolic pathway for tumor cells. However, the potential roles of glycolysis-related genes in renal cell carcinoma (RCC) have not been investigated.
Seven glycolysis-related gene sets were selected from MSigDB and were analyzed through GSEA. Using TCGA database, the glycolysis-related gene signature was constructed. Prognostic analyses were based on the Kaplan–Meier method. The cBioPortal database was employed to perform the mutation analyses. The CIBERSORT algorithm and TIMER database were used to determine the immunological effect of glycolytic gene signature. The expressions in protein level of eight glycolytic risk genes were determined by HPA database. Finally, qPCR, MTT and Transwell invasion assays were conducted to validate the roles of core glycolytic risk genes (CD44, PLOD1 and PLOD2) in RCC.
Four glycolysis-related gene sets were significantly enriched in RCC samples. The glycolytic risk signature was constructed (including CD44, PLOD2, KIF20A, IDUA, PLOD1, HMMR, DEPDC1 and ANKZF1) and identified as an independent RCC prognostic factor (HR = 1.204). Moreover, genetic alterations of glycolytic risk genes were uncommon in RCC (10.5%) and glycolytic risk signature can partially affect immune microenvironment of RCC. Six glycolytic risk genes (except for IDUA and HMMR) were over-expression in A498 and 786-O renal cancer cells through qPCR test. MTT and Transwell assays revealed that silencing of CD44, PLOD1 and PLOD2 suppressed the proliferation and invasion of renal cancer cells.
The glycolysis-related risk signature is closely associated with RCC prognosis, progression and immune microenvironment. CD44, PLOD1 and PLOD2 may serve as RCC oncogenes.
Renal cell carcinoma (RCC) is the third most common genitourinary cancer, accounting for approximatively 2.2% of new cancer cases worldwide . It has a high mortality rate, which results in almost 200,000 deaths/year, accounting for 1.8% of all cancer-related deaths . Radical nephrectomy is the main treatment for early stage RCC patients, with a 72% 5-year overall survival rate (OSR) . Unfortunately, approximatively 20–30% of patients have metastatic symptoms at the time of diagnosis, with only a 12% 5-year OSR . In addition, 20–30% of patients diagnosed with T1–2 stages suffered from metastasis within 1 to 2 years after surgery . Although molecular-targeted drugs improve survival time of metastatic patients, the median survival time is still less than 3 years . Therefore, it is still urgent and challenging to explore RCC molecular mechanisms and find potential therapeutic targets and prognostic biomarkers.
With the last decade progress in oncology, the reprogramming of energy metabolism has been considered as a new hallmark of cancers . Particularly, the tumor cells’ carbohydrate metabolism pattern was found to be extremely different from that of normal cells. Even in the presence of abundant oxygen, cancer cells still largely switch their metabolism to glycolysis, a metabolic process called Warburg effect . Glycolysis is an enzymatic reaction, that is collectively regulated by 3 limiting enzymes (including HK,PFK and PKM2) and some equilibrium enzymes of fructose-bisphosphatase (including ALDO, ENO, GAPDH) . Glycolysis-related genes, which encode these regulatory enzymes, were shown to be closely related to malignant progression and prognosis of many cancers, and therefore, are considered as potential therapeutic targets . Hexokinase (HK) is overexpressed in rectum , breast  and gastric cancers  and is associated with poor prognosis. Phosphofructokinase 1(PFK1), a key enzyme in glycolysis flux control , is also up-regulated in breast cancer and can promote malignant proliferation and invasion through p-AKT  activation. However, the potential effect of glycolytic genes on RCC prognosis and progression has not been investigated. Besides, the role of glycolysis in RCC is also not fully elucidated.
With the continuous development of bioinformatics and high-throughput technology, several biomarkers have been found to be closely associated with RCC prognosis and malignant progression . However, the predictive ability of single-gene biomarkers is limited in RCC prognosis analysis and multiple-genes signature maybe a better strategy . So far, only four studies have established the glycolytic risk signature in lung adenocarcinoma , endometrial cancer , hepatocellular carcinoma  and bladder cancer . Therefore, to explore the potential roles of glycolytic metabolism and glycolysis-related genes in RCC, a series of bioinformatic analyses were performed through TCGA database (609 samples). We found that glycolysis is significantly enriched in RCC and constructed a glycolysis-related risk signature. The novel risk signature was closely associated with RCC prognosis, progression and immune microenvironment. Moreover, its hub genes (CD44, PLOD1 and PLOD2) were proven to possess pro-cancer abilities through MTT and Transwell invasion assays. Our research provides a new method for RCC prognostic analysis and a new insight into the mechanisms of glycolytic genes in RCC progression.
The gene expression data and corresponding clinical information were obtained from the Cancer Genome Atlas (TCGA) public database. A total of 609 samples were analyzed, including 537 RCC and 72 normal or paracancerous samples. The clinical features of RCC patients were shown in Table 1. The data type of gene expression was transcriptome profiling and Gene Expression Quantification, and the type of clinical data was BCR-XML. The selected sample types were adenomas and adenocarcinomas.
Glycolysis-related gene sets and GSEA
Well-annotated gene sets, representing the universe of the biological processes, are critical for meaningful and insightful interpretation of large-scale genomic data . The Molecular Signatures Database (MSigDB) is a collection of annotated gene sets for GSEA software use (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) . Combination with the selecting strategies of previous glycolysis-related studies [22,23,24], we used ‘glycolysis’ and ‘glycolytic’ as the searching terms in the MSigDB database and furtherly screened out seven priori glycolysis-related gene sets. These gene sets included Hallmark Glycolysis, Go Glycolytic Process, KEGG Glycolysis Gluconeogenesis, Reactome Glycolysis, BioCarta Feeder Pathway and Module 306. The descriptions of these gene sets were shown in Supplementary Table 1.
GSEA is a computational method that determines whether an a priori defined set of genes shows statistically significant and concordant differences between two biological phenotypes . Enrichment statistics was based on weighted method. Phenotype labels were set as tumor versus normal samples. Permutation Type was phenotype. The number of permutations was set as 1000 and the gene sets were considered to significantly enriched in RCC samples, when the normalized enrichment score (NES) ≥1, nominal (NOM) p-value ≤0.05 and false discovery rate (FDR) q-value ≤0.25 were simultaneously satisfied. The detailed parameter settings of GSEA were presented in Supplementary Table 2.
Identification and functional analyses of glycolysis-related DEGs
TCGA data extraction and arrangement were performed by Perl (Practical Extraction and Report Language) version 5.28. The differential expression of glycolysis-related genes between tumor and normal samples was analyzed via the ‘Limma’ package in R software Ver3.6.2, in order to select differentially expressed genes (DEGs). When p-value < 0.05 and Log2FC absolute value ≥1, the gene was regarded as differentially expressive in tumor samples. Protein-protein interaction (PPI) networks of glycolysis-related genes were constructed using the String database  and Cytoscape software Ver3.4.0 .
To explore the biological function of glycolysis-related DEGs in RCC, we performed Gene Ontology (Go) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses through David database . The figures were drawn by the ggplot2 package (R software). The bioinformatics analysis of Go consisted in three categories, namely molecular function (MF), biological process (BP) and cellular components (CC). The pathways and biological functions with top ten enriched gene counts were selected for graphing.
Construction of the risk signature
Cox univariate analysis was used to preliminarily screen out the glycolysis-related genes which could affect RCC prognosis. Subsequently, the risk signature based on glycolysis-related genes was established through multivariate regression analysis. Since the expression levels of glycolysis-related genes, in the risk signature, were all high, Log2 transformation was used to normalize the expression profiles of these genes, which can prevent the generation of diminutive coefficients. The risk plot and the corresponding heatmap were drawn using the ‘pheatmap’ package in R software.
According to the glycolytic risk signature, the risk score of each RCC samples was further calculated. Subsequently, RCC samples were divided into high- and low-risk groups, according to the median of risk score. The prognostic difference between high- and low-risk groups was evaluated via the ‘survival’ package in R software. Moreover, cox univariate and multivariate analyses were successively performed to identify whether risk score was an RCC independent prognostic factor. To assess the applicative scope of glycolytic risk signature in RCC prognosis analysis, the prognostic difference was detected among RCC patients between different risk groups under the same clinical subgroup.
Receiver operating characteristic curve (ROC) was employed to assess predicting accuracy of glycolytic risk signature in RCC prognosis analysis. Decision curve analysis (DCA) was used to evaluate the clinical net benefit (NB) brought by the novel signature .
Mutation analyses of glycolysis-related genes
To better comprehend the genomics profiles of glycolysis-related genes, we performed mutation analyses using the cBioPortal online database . The OncoPrint tab summarizes genomic alterations in all queried genes across a sample set . The cancer types summary tab reveals mutation types and frequency information of queried genes.
Immune analyses of glycolysis-related genes
To determine the effect of glycolysis-related risk signature on immune microenvironment of RCC, immune infiltration analyses were performed through the CIBERSORT algorithm and TIMER database. The immune abundances of 22 leukocyte subtypes in each RCC samples were obtained by using the CIBERSORT algorithm. The landscape of immune distributions of each RCC samples were presented by ‘barplot’ package. The differential infiltration levels of 22 immune cells between high- and low-risk groups were visualized via the ‘vioplot’ package in R software.
The public database, Tumor Immune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/), was employed to assess the relationships between the somatic copy number alteration (SCNA) of eight risk signature genes and the infiltration levels of six immune cells (including B cells, CD4 + T cells, CD8 + T cells, neurphils, macrophases and dendritic cells).
Protein levels of glycolytic risk genes in the HPA database
The Human Protein Atlas (HPA) is a database with the aim to map all the human proteins in cells, tissues and organs using integration of various omics technologies (https://www.proteinatlas.org/) . The HPA database consists of six separate parts. Among that, the Tissue and the Pathology Atlas provide information regarding the expression profiles of specified genes in normal and tumor tissues on protein levels. All images of tissues in HPA database are stained by immunohistochemistry.
Two renal cancer cell lines (786-O and A498) and one normal human renal tubular epithelial cell line (HK-2) were purchased from the Institute of Biochemistry and Cell Biology at the Chinese Academy of Sciences (Shanghai, China). Cancer cells were cultured in RPMI-1640 medium (GIBCO-BRL, Invitrogen, Carlsbad, CA, USA) with 10% fetal bovine serum (FBS, Hyclone, Logan, Utah, USA). HK-2 cells were cultured in keratinocyte medium (KM, ScienCell, San Diego, California, USA) plus 1% keratinocyte growth supplement (KGS, Scien-Cell, San Diego, California, USA). All mediums were treated with 100 U/ml penicillin and 100μg/ml streptomycinm, and all cells were incubated at 37 °C with 5% CO 2.
Quantitative real-time PCR (qRT-PCR)
Total RNA was extracted using TRIzol reagent (G-Clone, Beijing, China). The cDNA was synthesized with the PrimeScript RT reagent Kit with gDNA Eraser (Takara, Japan). qRT-PCR was performed using the TB Green Premix Ex Taq II (Takara, Japan) in ABI 7500 Real-Time PCR instrument. The primers were listed in Supplementary Table 3. GAPDH was used as an internal control. The relative mRNA levels were calculated based on 2−ΔΔCt method.
The expression in mRNA level of CD44, PLOD1 and PLOD2 were silenced by specific small interfering RNAs (siRNAs). Si-CD44, si-PLOD1, si-PLOD2 and si-NC (negative control) were synthesized by GenePharma Biotechnology (Shanghai, China). 786-O and A498 cells were transfected by these siRNAs through the Lipofectamine 2000 (Invitrogen, Thermo Fisher, Waltham, MA, USA) Further experiments were performed 48 h after transfection.
Transfected cells were seeded into 96-well plate (5 × 103 cells per well) and incubated for 24, 48 and 72 h. Every detecting point (24 h/each time), 10 μL MTT reagent (Life science, NY, USA) was added to each well and incubated for 4 h. After removing the medium and washing each well by PBS, 100 μL DMSO was added into each well and the absorbance at 490 nm was measured by a microplate reader (ThermoFisher, Waltham, MA, USA).
Transwell invasion assay
The polycarbonate membrane was coated with the matrigel (200 ng/mL; BD Biosciences, Franklin Lakes, NJ, USA). Then transfected cancer cells were seeded with serum-free RPMI-1640 medium into the upper chambers (1 × 104 cells/well). The RPMI-1640 medium with 10% FBS was added into the bottom chambers. After 24 h incubation at 37 °C, the invasive cells adhering to the lower surface of the membrane were fixed by 4% paraformaldehyde and stained with 0.1% crystal violet. The invasive cells were photographed by a microscope.
All statistical analyses and graphing were performed by the R software (version 3.6.2) or GraphPad Prism (version 8.0.1). The Student’s t test was used to compare the expressive difference of glycolysis-related genes. The Kolmogorov–Smirnov test was used to evaluate the relationships between glycolytic risk score and RCC clinicopathological characteristics. Additionally, the survival analysis was based on the Kaplan–Meier method. A p-value < 0.05 was considered significant.
Glycolytic metabolism is significantly enriched in RCC samples
Five hundred thirty-seven RCC samples from TCGA database were included in our bioinformatic analyses (Table 1). Seven glycolysis-related gene sets were selected from the MSigDB for GSEA and the analytical results were shown in Fig. 1. Four gene sets were significantly enriched in RCC samples, namely Biocarta feeder pathway, Biocarta glycolysis pathway, Reactome glycolysis and Hallmark glycolysis (Table 2). Glycolysis-related gene sets (44.4%, 4/9) and glycolysis-related genes (77.0%, 261/339) were greatly enriched in RCC, suggesting that the glycolytic metabolism was prevalent in RCC, which is consistent with the general metabolic hallmark of cancers. After the GSEA preliminary screening, a total of 261 glycolysis-related genes from the 4 enriched gene sets were speculated to play important roles in RCC and were therefore selected for further analyses.
Glycolysis-related genes are differential expressed in RCC samples and regulate some cancer-related biofunction
To further screen out the hub genes that contribute to RCC malignant progression, we assessed the differential expression of glycolysis-related genes between RCC and normal samples. A total of 75 out of 261 glycolytic genes were significantly differentially expressed in RCC. The PPI network of 75 glycolysis-related DEGs was constructed via the String database and Cytoscape software (Fig. 2a). Among three kinds of glycolysis limiting enzymes, HK3 was up-regulated in RCC, while the expression of PKM family members were bidirectional regulated (PFKP and PFKFB4 were up-regulated; PFKFB1–3 were down-regulated). However, there was no expressive difference of PKM between RCC and normal samples.
GO enrichment analysis revealed that the subcellular functional localization of glycolytic genes was cytoplasm, and these genes may regulate glycolytic metabolism through the exosome pathway (Fig. 2b). Meanwhile, these glycolytic genes also engaged in cell proliferation and apoptosis. KEGG enrichment analysis showed that glycolysis-related genes not only involved in glycolysis, pentose phosphate and fructose metabolism pathways, but also were closely related to some cancer-related signal pathways, such as HIF-1, renal cell carcinoma and AMPK pathways (Fig. 2c).
Eight glycolytic DEGs constitute the glycolysis-related risk signature in RCC
A single gene generally cannot accurately predict cancer prognosis; however, a gene signature, consisting of multiple genes, could improve the accuracy of prognosis analysis. Through univariate regression analysis, 23 out of the 75 glycolysis-related DEGs were found to greatly affect RCC prognosis (Table 3). Subsequently, the glycolytic risk signature was constructed via multivariate regression analysis (Table 3). Risk score = 0.239 × CD44 + 0.203 × PLOD2 + 0.695 × KIF20A + 0.527 × IDU4–0.265 × PLOD1–0.998 × HMMR + 0.886 × DEPDC1 + 0.199 × ANKZF1 (the gene expression was normalized with log2 transformation). Eight risk signature genes were all over-expressed in RCC samples (Supplementary Figure 1) and their overexpression all led worse survival outcomes (Fig. 3c-j), which indicated that they may serve as oncogenes.
The risk score of each RCC patient was calculated using the risk formula and the patients were divided into high-and low-risk groups according to median of risk score (Fig. 3a). The proportion of death events in high-risk group was significantly higher than that in low-risk group (Fig. 3b). Moreover, patients in high-risk group were often associated with undesirable clinicopathological features (Fig. 3k-o). These results indicate that the risk signature based on glycolytic genes not only affect RCC prognosis but also is closely associated with RCC malignant progression. Besides, glycolytic risk signature also possessed a good predictive accuracy in prognosis analysis (AUC = 0.729) (Fig. 3p).
Glycolysis-related risk signature is valuable for RCC prognostic analyses
Comparing with the low-risk group, the high-risk group led to worse 3-year (62.9% vs 86.2%) and 5-year overall survival rates (42.4% vs 78.2%) (Fig. 4a). Moreover, the DCA curve indicated that taking glycolysis risk signature into the RCC prognostic analysis can bring significant clinical benefit when making clinical decision (Fig. 4b).
Through univariate and multivariate regression analyses, glycolytic risk score (HR = 1.204, p < 0.001) was identified as an independent RCC prognostic factor (Fig. 4cd). To further determine the applicable scope of glycolytic risk signature in RCC prognostic analysis, the prognostic difference between different risk groups under the same clinical subgroup, was compared (Fig. 4g-r). Glycolytic risk signature could sensitively distinguish the prognostic differences in patients with age ≥ 65, age < 65, Grade 2–3, M0–1 stages, N0 stage, clinical stage I-IV and T1–3 stages. Due to lack of samples, glycolytic risk signature did not detect the prognostic difference of RCC patients with Grade1/4, N1 and T4 stages (Supplementary Figure 2A-D). Therefore, there is ample evidence reiterates that glycolytic risk signature is a beneficial and universally applicable method for RCC prognostic analysis.
Genetic alteration of glycolytic risk signature is uncommon in RCC
To better understand the genomic characteristics of glycolytic genes in RCC, the cBioPortal database was used for mutation analysis. A summary of the genetic alterations in eight risk signature genes is shown in Fig. 5a. Although these glycolysis-related genes were significantly differentially expressed in RCC samples, their mutation frequency was not high. Based on TCGA database, glycolytic genes were only mutated in 37 out of the 354 RCC samples (10.5%).
HMMR and KIF20A possessed the highest mutation frequency (7%), while IDUA, PLOD1 and DEPDCA did not generate any mutations (Fig. 5a). Moreover, we also analyzed mutational type. Different glycolytic risk genes were provided with different mutational types (Fig. 5b-g). For example, deep deletion only occurred in the genetic alterations of ANKZF1 (Fig. 5g) and CD44 harbored only one kind of alteration pattern, namely mutation (Fig. 5b). Overall, amplification was the most common mutational form (Fig. 5h) and the rarest form was deep deletion (3%). All these results revealed that the glycolytic risk signature was relatively conserved in RCC and may not reflect genetic predisposition of RCC.
Glycolytic risk signature could affect immune microenvironment of RCC to some extent
We further analyzed the effect of eight glycolytic risk genes on RCC immune microenvironment. The immune abundances of 22 leukocyte subtypes in each RCC samples were exhibited in Supplementary Figure 3. The immune proportions of each sample were diverse. Moreover, the difference of immune abundance between high- and low-risk groups was determined (Fig. 6a). High risk groups brought increased infiltration levels of T cells follicular helper (p<0.001) and T cells regulatory (Tregs) (p<0.001) in RCC; Inversely, it could decrease the immune abundance of Macrophages M2 (p = 0.049), Dendritic cells resting (p = 0.022) and Mast cells resting (p<0.001) (Fig. 6a). However, the risk signature cannot bring any change on immune infiltration levels of CD4 and CD8 T cells. Combination with some immunological studies [33,34,35,36,37], the basic biofunction in immune response of these affected lymphocytes and the final effect of their alterations against antitumor immunological process were shown in Table 4. Alterations of immune abundance caused by high glycolytic risk did not always hinder antitumor immune process (Table 4).
Moreover, the relationships between SCNA of eight risk signature genes and the infiltration levels of six immune cells were explored through TIMER database (Fig. 6b-i). The arm-level gain in copy number of ANZKF1, CD44 and PLOD2 could lead an increment on infiltration levels of CD4 T cells (Fig. 6bci). However, as the most common mutational type (Fig. 5h), the amplification of these genes could not affect the infiltration levels of CD4 and CD8 T cells. Besides, the amplified alteration of other risk genes did not contribute to the increasing infiltration levels of macrophages, CD4 and CD8 T cells. In a word, the effect of glycolytic risk signature on immune microenvironment was limit and complex.
Glycolytic risk genes differentially express in protein levels
The histological expressions of glycolysis-related risk genes in normal and tumor tissues were exhibited in Fig. 7. There was significantly expressive difference of these genes in protein levels between normal and renal cancer tissues. Among that, the protein expression of PLOD1, PLOD2, KIF20A, IDUA and DEPDC1 were obviously up-regulated; However, discernable expressive difference was not found in CD44, HMMR and ANKZF1. The protein expression of CD44 and ANKZF1 were low in both normal and tumor tissues. And that of HMMR were even not detected.
CD44, PLOD1 and PLOD2 can promote the proliferation and invasion of renal cancer cells
The PPI network of eight glycolytic risk genes was constructed through the String database. As shown in Fig. 8a, CD44, PLOD1 and PLOD2 were in the central of the network suggesting these genes maybe the core genes among glycolysis-related risk signature. Therefore, we validated the biofunctions of CD44, PLOD1 and PLOD2 in renal cancer cells through further vitro experiments.
Except for IDUA and HMMR, other six risk genes were all up-regulated in renal cancer cells (Fig. 8b-d and supplementary Figure 4). To investigate the biofunctions of CD44, PLOD1 and PLOD2 in RCC, we synthesized siRNAs of these genes to decrease their expression. The qPCR tests confirmed that their expressions were significantly down-regulated in transfected renal cancer cells (Fig. 8e-g). MTT and Transwell invasion assays were employed to evaluate the proliferative and invasive abilities of cells. As expected, empty vector transfection (NC group) did not affect cell viability. However, silencing of CD44, PLOD1 and PLOD2 inhibited the proliferative abilities of A498 and 789-O cells (Fig. 8h-m). Moreover, knockdown of CD44, PLOD1 and PLOD2 can suppress the invasion of A498 and 789-O cells, compared with those transfected with si-NC (Fig. 9). As a result, we ascertained that CD44, PLOD1 and PLOD2 played a pro-cancer role in RCC through vitro experiments, reiterating that glycolytic risk signature was closely associated with RCC progression.
Since Otto Warburg first discovered in 1924 that tumor cells preferred to engage in aerobic glycolysis, this cellular process has gradually been confirmed as one of the cancer hallmarks, that is also known as the Warburg effect . In the present study, four glycolysis-related gene sets from the MSigDB were significantly enriched in RCC samples, suggesting that aerobic glycolysis is also prevalent in RCC. However, the reasons that tumor cells preferred to utilize glycolysis, an inefficient metabolic pattern, rather than oxidative phosphorylation for their proliferation, remain entirely unclear. Based on previous studies, we speculated there are three reasons for this tumoral metabolic preference.
First of all, proliferating cells have important metabolic requirements beyond ATP . Cell proliferation requires a large number of nucleotides, amino acids, and lipids. A molecule of glucose provides a maximum of 36 molecules of ATP, but only 2 molecules of NADPH and 6 molecules of carbon . During cell proliferation, massive glucose consumption is not a prerequisite for ATP production, but is essential for the synthesis of biological products that are needed for proliferation. For instance, during synthesizing the palmitate, an important component of cellular membranes, a molecule of glucose provides five times ATP as much as needed, but less than 20% carbon and NADPH as required . Moreover, the pentose phosphate pathway, a branch from glycolysis, is the main source of nucleotides and NADPH during proliferation process . Similarly, our KEGG analysis also indicated glycolytic genes significantly enriched in the pentose phosphate pathway. Meanwhile, NADPH supplied by the pentose phosphate pathway can protect cancer cells from oxidative stress . Thus, the glycolytic metabolism is inefficient for ATP production but efficient for cell proliferation. Secondly, aerobic glycolysis can induce tumor cells’ therapeutic resistance. Classic chemotherapy treatments generally target rapidly dividing cells; however, quiescent tumor cell populations evade therapeutic destruction. It has been demonstrated that proliferating tumor cells can induce the accumulation of lactic acid through prompting glycolysis metabolism, which leads to cellular acidosis and quiescence [42, 43]. When external treatments are terminated, quiescent cells can reenter proliferative state and complete their therapeutic evasion . Thirdly, glycolysis can maintain a high ratio of ATP/ADP, which leads to the continuous activation of cell proliferation [45, 46]. Although glycolysis is not a highly efficient ATP producer, it can regulate the activity of adenylate kinases to buffer the declining ATP production [38, 47]. In summary, the glycolytic metabolism can provide stimulating signals and a sufficient biomass for cell proliferation, and assist in tumor cells therapeutic resistance. The GSEA results also revealed that glycolysis is prevalent in RCC, which is consistent with previous studies .
After screening, we found that eight glycolytic risk genes were differentially expressed in RCC; however, their mutation frequency is not high (37/354, 10.5%). This suggests that the aberrant expression of these genes may be a result of post-transcriptional regulations or translation modifications. For example, over-expressed DEPDC1 is negatively regulated by miR-26b, which facilitates cell proliferation in triple negative breast cancer (TNBC) . Moreover, the loss of the von Hippel-Lindau (VHL) gene is the most prominent genetic alteration in RCC, that is associated with over 80% of RCC cases . Downregulated VHL elevated the accumulation of HIF-1 and led to the activation of hypoxia-inducible genes, such as PLOD1 and PLOD2, which is a key metabolic step for the transition into aerobic glycolysis [51,52,53]. Similarly, KEGG analysis also indicated that glycolytic genes are involved in the HIF-1 signaling pathway (Fig. 2d), suggesting that the screened genes may play a crucial role in glycolytic metabolic switch. Mutation analysis is valuable for the identification of tumorous susceptible genes, which can be used as favorable markers for early diagnosis and therapeutic target prediction . In the present study, the highest mutational frequency among glycolytic risk genes was only 7%, which demonstrates that the glycolytic risk signature may not serve as RCC susceptibility genes.
With the development of high-throughput genomics, many genes have been identified as potential biomarkers of tumor prognosis and progression, and as therapeutic targets. In the present study, we constructed a glycolysis-related risk signature and validated the biofunctions of its hub genes (CD44, PLOD1 and PLOD2) in RCC. CD44, a non-kinase transmembrane glycoprotein, is regard as a marker of cancer stem cell (CSC) and can regulate the properties of CSCs, including cellular plasticity, self-renewal, invasiveness and treatment resistance . Meanwhile, Gan Yu et al. demonstrated that miRNA-34a suppressed the proliferation and metastasis of renal cancer cells by inhibiting CD44 . In the current research, silencing CD44 also can weaken proliferative and invasive abilities of renal cancer cells, which coincided with previous research . PLOD1 and PLOD2 are the key enzymes of collagen synthesis and extracellular matrix formation, therefore can drive tumorigenesis and angiogenesis in many cancers . Although PLOD1 was proven to promote aggressiveness of bladder cancer cells , the potential roles of PLOD1 and PLOD2 in RCC were not elucidated. Through MTT and Transwell invasion assays, we confirmed the pro-cancer abilities of PLOD1 and PLOD2, which provided new basis for the molecular mechanism of RCC tumorigenesis.
The glycolysis-related gene signature was confirmed as an independent prognostic factor of RCC and can distinguish the prognostic difference of the patients under the same clinical subgroup (Fig. 4). Therefore, we believed that the novel risk signature was valuable for RCC clinical assessment. Firstly, glycolytic risk signature is an important supplement to RCC prognosis analysis. Single TNM staging system cannot accurately predict RCC prognosis. Correa et al. have found that the predictive ability of TNM model was significant variable over time and its predictive accuracy was not satisfactory (C-index = 0.60) . Nevertheless, as DCA curve displayed (Fig. 4b), taking glycolysis risk signature into the RCC prognostic analysis can add net benefit when making clinical decision. Secondly, dividing RCC patients into different risk groups drives individualized treatment strategies. When patients are classified into high-risk group, their follow-up time could be shortened and some adjuvant therapy could be attempted.
Tumor immune microenvironment commonly reflects the immune status of cancer, which could give support to formulate therapeutic strategy. CD8 T cells could differentiate into cytotoxic T cells, which exhibit cytotoxicity against tumor cells . CD4 T cells could assist antitumor immunity and deficient CD4 T cells reduce the response of cytotoxic T lymphocytes (CTLs) . Therefore, these two types of lymphocytes are the core ingredients for tumor cellular immunity. However, the glycolytic gene signature cannot affect infiltration levels of CD4 and CD8 T cells, which reveals that the impact of risk signature on immune microenvironment of RCC is limit. Besides, alterations of immune abundance caused an intricate effect on antitumor immune process (Table 4). For example, Mast cells are bidirectional to tumor progression . On one hand, Mast cells promote cancer growth, stimulate neoangiogenesis and remodel tissue through releasing potent proangiogenic factors such as vascular endothelial growth factor (VEGF) and basic fibroblast growth factor (bFGF) . On the other hand, they can suppress cell proliferation, inhibit immunologic stimulation and cell mobility by secreting chymase, tryptase, TNF-α, IL-1 and IL-6 . Hence, we difficultly determine the final effect of decreasing mast cells (in high risk group) on antitumor immunity.
It’s believed that our findings could bring some advance to RCC field. Firstly, potential therapeutic targets. The role of Warburg’s effect in RCC has been widely confirmed and blocking the VHL-HIF-Glycolysis axis has been considered as a potential therapeutic strategy for renal cancer [62, 63]. For example, 2-deoxy-D-glucose, a glycolysis inhibitor, can not only kill various pathological subtypes of renal cancer cells, including clear cell RCC (ccRCC), papillary RCC and the rare subtype chromophobe RCC, but also enhance susceptibility of ccRCC to pazopanib treatment . In the present study, we demonstrated that inhibitions of three core glycolytic genes suppressed the proliferation and invasion of renal cancer cells, indicating that these genes may serve as potential therapeutic targets of RCC. Secondly, assisting in prognostic analysis. Using TCGA database, we established a novel glycolysis-related prognostic model, which could assess the risk levels of RCC patients with the same TNM stage and contribute to make prognostic assessment more accurately. Thirdly, prompt for genetic regulatory mode. VHL mutations occur in more than 80% RCC cases . Invalidation of VHL factor stimulates the expression of glycolytic genes and increases the glycolytic flux through stabilizing HIF . However, in mutation analysis, we found that only 10.5% of RCC samples harbored with the mutations of glycolytic genes, suggesting that the ectopic expression of glycolytic genes in RCC may be accomplished by epigenetic modification rather than mutation. For example, miR-122-5p can promote proliferation and invasion of renal cancer cells by targeting PKM2 (a glycolytic limiting enzyme) . LncARSR-miR-34a-5p-HK1 axis (HK is another glycolytic limiting enzyme) facilitates the progression of rectal cancer by promoting glycolysis metabolism . In summary, our discoveries provide new insights into the roles of glycolytic genes in RCC from the perspectives of prognosis, molecular mechanism and regulatory mode.
Naturally, there are some limitations in our study. Firstly, the biofunctions of glycolysis-related genes in RCC have not been verified through vivo models and clinical samples. Secondly, glycolytic prognostic model should be tested in actual clinical cases. Thirdly, we did not assess the impact of fasting blood glucose (FBG) and anti-diabetes drugs on final results. According to some metabolomic discoveries and oncological researches, we found that some hypoglycemic drugs and patients’ FBG certainly regulated the tumor progression, which may affect the results of this study. For example, Metformin was proven to possess anti-cancer ability in acute myeloid leukemia (AML)  and lung cancer  by inhibiting glycolytic process. Ashton TM et al. suggested that Metformin may be a potential antitumor drug for its antagonistic action against oxidative phosphorylation (OXPHOS) and glycolysis . Regarding FBG, Ayush Sharma et al. confirmed that high level of FBG was associated with large tumor volume and undesirable histological grade in pancreatic cancer . In gastric cancer, elevated FBG and high SNHG8 expression resulted in poor survival outcomes after radical gastrectomy . Unfortunately, TCGA database does not provide information on FBG and anti-diabetes drugs, which we believe is indeed one of the limitations of the present study.
To the best of our knowledge, we constructed a novel glycolysis-related risk signature in RCC for the first time. The glycolytic risk signature played a crucial role in progression, prognosis and immune microenvironment of RCC. Moreover, three core glycolytic risk genes, CD44, PLOD1 and PLOD2, were capable of promoting the proliferation and invasion of renal cancer cells. Our findings can provide new inspirations for RCC prognostic analysis and treatment.
Availability of data and materials
The datasets used and/or analyzed in the current study are available from the corresponding author upon reasonable request.
Renal cell carcinoma
The Cancer Genome Atlas
Gene-set enrichment analysis
Overall survival rate
The molecular signatures database
Normalized enrichment score
False discovery rate
Differentially expressed genes
Kyoto encyclopedia of genes and genomes
Procollagen-lysine, 2-oxoglutarate 5-dioxygenases
Hyaluronan mediated motility receptor
DEP domain containing 1
Kinesin family member 20A
Cytotoxic T lymphocytes
The Human Protein Atlas
Decision curve analysis
Receiver operating characteristic curve
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68:7–30.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68:394–424.
De P, Otterstatter MC, Semenciw R, Ellison LF, Marrett LD, Dryer D. Trends in incidence, mortality, and survival for kidney cancer in Canada, 1986-2007. Cancer Causes Control. 2014;25:1271–81.
Gupta K, Miller JD, Li JZ, Russell MW, Charbonneau C. Epidemiologic and socioeconomic burden of metastatic renal cell carcinoma (mRCC): a literature review. Cancer Treat Rev. 2008;34:193–205.
Gray RE, Harris GT. Renal Cell Carcinoma: Diagnosis and Management. Am Fam Physician. 2019;99:179–84.
Greef B, Eisen T. Medical treatment of renal cancer: new horizons. Br J Cancer. 2016;115:505–16.
Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144:646–74.
Torresano L, Nuevo-Tapioles C, Santacatterina F, Cuezva JM. Metabolic reprogramming and disease progression in cancer patients. Biochim Biophys Acta Mol Basis Dis. 1866;165721:2020.
Dobbelstein M, Moll U. Targeting tumour-supportive cellular machineries in anticancer drug development. Nature reviews. Drug Discov. 2014;13:179–96.
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.
Coelho RG, Calaça IC, Celestrini DM, et al. Hexokinase and phosphofructokinase activity and intracellular distribution correlate with aggressiveness and invasiveness of human breast carcinoma. Oncotarget. 2015;6:29375–87.
Gao Y, Xu D, Yu G, Liang J. 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.
Bartrons R, Simon-Molas H, Rodríguez-García A, et al. Fructose 2,6-Bisphosphate in Cancer cell metabolism. Front Oncol. 2018;8:331.
Peng F, Li Q, Sun JY, Luo Y, Chen M, Bao Y. PFKFB3 is involved in breast cancer proliferation, migration, invasion and angiogenesis. Int J Oncol. 2018;52:945–54.
Czarnecka AM, Kukwa W, Kornakiewicz A, Lian F, Szczylik C. Clinical and molecular prognostic and predictive biomarkers in clear cell renal cell cancer. Future oncology (London). 2014;10:2493–508.
Wang ZH, Zhang YZ, Wang YS, Ma XX. Identification of novel cell glycolysis related gene signature predicting survival in patients with endometrial cancer. Cancer Cell Int. 2019;19:296.
Zhang L, Zhang Z, Yu Z. Identification of a novel glycolysis-related gene signature for predicting metastasis and survival in patients with lung adenocarcinoma. J Transl Med. 2019;17:423.
Lu C, Fang S, Weng Q, et al. Integrated analysis reveals critical glycolytic regulators in hepatocellular carcinoma. Cell Commun Signal. 2020;18:97.
Zhang C, Gou X, He W, Yang H, Yin H. A glycolysis-based 4-mRNA signature correlates with the prognosis and cell cycle process in patients with bladder cancer. Cancer Cell Int. 2020;20:177.
Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (Oxford). 2011;27:1739–40.
Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–25.
Liu J, Li S, Feng G, et al. Nine glycolysis-related gene signature predicting the survival of patients with endometrial adenocarcinoma. Cancer Cell Int. 2020;20:183.
Liu C, Li Y, Wei M, Zhao L, Yu Y, Li G. Identification of a novel glycolysis-related gene signature that can predict the survival of patients with lung adenocarcinoma. Cell Cycle (Georgetown). 2019;18:568–79.
Liu Y, Yin S. A novel prognostic index based on the analysis of glycolysis-related genes in head and neck squamous cell carcinomas. J Oncol. 2020;2020:7353874.
Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–50.
Franceschini A, Szklarczyk D, Frankild S, et al. STRING v9.1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41:D808–15.
Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.
Huang DW, Sherman BT, Tan Q, et al. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8:R183.
Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Med Decis Making. 2006;26:565–74.
Gao J, Aksoy BA, Dogrusoz U, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6:l1.
Cerami E, Gao J, Dogrusoz U, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2:401–4.
Uhlén M, Fagerberg L, Hallström BM, et al. Proteomics. Tissue-based map of the human proteome. Science (New York). 2015;347:1260419.
Vinuesa CG, Linterman MA, Yu D, MacLennan IC. Follicular helper T cells. Annu Rev Immunol. 2016;34:335–68.
Juang CM, Hung CF, Yeh JY, et al. Regulatory T cells: potential target in anticancer immunotherapy. Taiwan J Obstet Gynecol. 2007;46:215–21.
Italiani P, Boraschi D. From monocytes to M1/M2 macrophages: Phenotypical vs. Funct Differ Front Immunol. 2014;5:514.
Brossart P, Wirths S, Brugger W, Kanz L. Dendritic cells in cancer vaccines. Exp Hematol. 2001;29:1247–55.
Dyduch G, Kaczmarczyk K, Okoń K. Mast cells and cancer: enemies or allies? Pol J Pathol. 2012;63:1–7.
Heiden MGV, Cantley LC, Thompson CB. Understanding the Warburg effect: the metabolic requirements of cell proliferation. Science (New York). 2009;324:1029–33.
Opdenakker G, Rudd PM, Ponting CP, Dwek RA. Concepts and principles of glycobiology. FASEB J. 1993;7:1330–7.
Patra KC, Hay N. The pentose phosphate pathway and cancer. Trends Biochem Sci. 2014;39:347–54.
Pinthus JH, Whelan KF, Gallino D, Lu JP, Rothschild N. Metabolic features of clear-cell renal cell carcinoma: mechanisms and clinical implications. Can Urol Assoc J. 2011;5:274–82.
Tindall MJ, Dyson L, Smallbone K, Maini PK. Modelling acidosis and the cell cycle in multicellular tumour spheroids. J Theor Biol. 2012;298:107–15.
Smallbone K, Gavaghan DJ, Maini PK, Brady JM. Quiescence as a mechanism for cyclical hypoxia and acidosis. J Math Biol. 2007;55:767–79.
Yeh AC, Ramaswamy S. Mechanisms of Cancer cell dormancy--another Hallmark of Cancer? Cancer Res. 2015;75:5014–22.
Christofk HR, Vander Heiden MG, Harris MH, et al. The M2 splice isoform of pyruvate kinase is important for cancer metabolism and tumour growth. Nature. 2008;452:230–3.
DeBerardinis RJ, Lum JJ, Hatzivassiliou G, Thompson CB. The biology of cancer: metabolic reprogramming fuels cell growth and proliferation. Cell Metab. 2008;7:11–20.
Panayiotou C, Solaroli N, Karlsson A. The many isoforms of human adenylate kinases. Int J Biochem Cell Biol. 2014;49:75–83.
Godinot C, de Laplanche E, Hervouet E, Simonnet H. Actuality of Warburg's views in our understanding of renal cancer metabolism. J Bioenerg Biomembr. 2007;39:235–41.
Zhang L, Du Y, Xu S, et al. DEPDC1, negatively regulated by miR-26b, facilitates cell proliferation via the up-regulation of FOXM1 expression in TNBC. Cancer Lett. 2019;442:242–51.
Nickerson ML, Jaeger E, Shi Y, et al. Improved identification of von Hippel-Lindau gene alterations in clear cell renal tumors. Clin Cancer Res. 2008;14:4726–34.
Gordan JD, Thompson CB, Simon MC. HIF and c-Myc: sibling rivals for control of cancer cell metabolism and proliferation. Cancer Cell. 2007;12:108–13.
Courtnay R, Ngo DC, Malik N, Ververis K, Tortorella SM, Karagiannis TC. Cancer metabolism and the Warburg effect: the role of HIF-1 and PI3K. Mol Biol Rep. 2015;42:841–51.
Gilkes DM, Bajpai S, Wong CC, et al. Procollagen lysyl hydroxylase 2 is essential for hypoxia-induced breast cancer metastasis. Mol Cancer Res. 2013;11:456–66.
Pearlman R, Frankel WL, Swanson B, et al. Prevalence and Spectrum of Germline Cancer susceptibility gene mutations among patients with early-onset colorectal Cancer. JAMA Oncol. 2017;3:464–71.
Yan Y, Zuo X, Wei D. Concise review: emerging role of CD44 in Cancer stem cells: a promising biomarker and therapeutic target. Stem Cells Transl Med. 2015;4:1033–43.
Yu G, Li H, Wang J, et al. miRNA-34a suppresses cell proliferation and metastasis by targeting CD44 in human renal carcinoma cells. J Urol. 2014;192:1229–37.
Qi Y, Xu R. Roles of PLODs in collagen synthesis and Cancer progression. Front Cell Dev Biol. 2018;6:66.
Yamada Y, Kato M, Arai T, et al. Aberrantly expressed PLOD1 promotes cancer aggressiveness in bladder cancer: a potential prognostic marker and therapeutic target. Mol Oncol. 2019;13:1898–912.
Correa AF, Jegede O, Haas NB, et al. Predicting renal Cancer recurrence: defining limitations of existing prognostic models with prospective trial-based validation. J Clin Oncol. 2019;37:2062–71.
Iwahori K. Cytotoxic CD8 lymphocytes in the tumor microenvironment. Adv Exp Med Biol. 2020;1224:53–62.
Borst J, Ahrends T, Bąbała N, Melief CJM, Kastenmüller W. CD4 T cell help in cancer immunology and immunotherapy. Nat Rev Immunol. 2018;18:635–47.
Sudarshan S, Karam JA, Brugarolas J, et al. Metabolism of kidney cancer: from the lab to clinical practice. Eur Urol. 2013;63:244–51.
Gill KS, Fernandes P, O'Donovan TR, et al. Glycolysis inhibition as a cancer treatment and its role in an anti-tumour immune response. Biochim Biophys Acta. 2016;1866:87–105.
Simon AG, Esser LK, Ellinger J, et al. Targeting glycolysis with 2-deoxy-D-glucose sensitizes primary cell cultures of renal cell carcinoma to tyrosine kinase inhibitors. J Cancer Res Clin Oncol. 2020;146:2255–65.
Wang S, Zheng W, Ji A, Zhang D, Zhou M. Overexpressed miR-122-5p promotes cell viability, proliferation, migration and glycolysis of renal Cancer by negatively regulating PKM2. Cancer Manag Res. 2019;11:9701–13.
Li S, Zhu K, Liu L, Gu J, Niu H, Guo J. lncARSR sponges miR-34a-5p to promote colorectal cancer invasion and metastasis via hexokinase-1-mediated glycolysis. Cancer Sci. 2020;111:3938–52.
Chen HL, Ma P, Chen YL, et al. Effect of metformin on proliferation capacity, apoptosis and glycolysis in K562 cells. Zhongguo Shi Yan Xue Ye Xue Za Zhi. 2019;27:1387–94.
Zhou X, Liu S, Lin X, et al. Metformin inhibit lung Cancer cell growth and invasion in vitro as well as tumor formation in vivo partially by activating PP2A. Med Sci Monit. 2019;25:836–46.
Ashton TM, McKenna WG, Kunz-Schughart LA, Higgins GS. Oxidative phosphorylation as an emerging target in Cancer therapy. Clin Cancer Res. 2018;24:2482–90.
Sharma A, Smyrk TC, Levy MJ, Topazian MA, Chari ST. Fasting Blood Glucose Levels Provide Estimate of Duration and Progression of Pancreatic Cancer Before Diagnosis. Gastroenterology. 2018;155:490–500.e492.
Lin Y, Hu D, Zhou Q, Lin X, Lin J, Peng F. The fasting blood glucose and long non-coding RNA SNHG8 predict poor prognosis in patients with gastric carcinoma after radical gastrectomy. Aging (Albany). 2018;10:2646–56.
The authors would like to thank Editage (https://www.editage.cn) for providing linguistic assistance. The first author also would like to thank Mrs. Xin Wang for encouragement.
This work was supported by National Natural Science Foundation of China (No. 81770689). The funder Tie Chong designed the study. The funding body provides support in experimental validation.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Differential expression of eight risk signature genes in RCC samples. RCC, renal cell carcinoma.
The glycolysis-related risk signature could not distinguish the prognostic difference of RCC patients with G1, G4, N1 and T4 stages. RCC, renal cell carcinoma. G, grade.
The distribution of immune abundance of 22 leukocyte subtypes in each RCC samples. RCC, renal cell carcinoma.
The mRNA differential expression of five glycolytic genes in different cell lines. * means p<0.05; * * means p<0.01.
The description of glycolysis-related gene sets.
The detailed parameter settings of GSEA. GSEA, Gene-set enrichment analysis; MSigDB, The Molecular Signatures Database; NA, not applicable; Meandiv, Mean deviation.
All primers for qRT-PCR.
About this article
Cite this article
Xu, F., Guan, Y., Xue, L. et al. The effect of a novel glycolysis-related gene signature on progression, prognosis and immune microenvironment of renal cell carcinoma. BMC Cancer 20, 1207 (2020). https://doi.org/10.1186/s12885-020-07702-7