- Open Access
Alteration in glycolytic/cholesterogenic gene expression is associated with bladder cancer prognosis and immune cell infiltration
BMC Cancer volume 22, Article number: 2 (2022)
Oncogenic metabolic reprogramming contributes to tumor growth and immune evasion. The intertumoral metabolic heterogeneity and interaction of distinct metabolic pathways may determine patient outcomes. In this study, we aim to determine the clinical and immunological significance of metabolic subtypes according to the expression levels of genes related to glycolysis and cholesterol-synthesis in bladder cancer (BCa).
Based on the median expression levels of glycolytic and cholesterogenic genes, patients were stratified into 4 subtypes (mixed, cholesterogenic, glycolytic, and quiescent) in an integrated cohort including TCGA, GSE13507, and IMvigor210. Clinical, genomic, transcriptomic, and tumor microenvironment characteristics were compared between the 4 subtypes.
The 4 metabolic subtypes exhibited distinct clinical, molecular, and genomic patterns. Compared to quiescent subtype, mixed subtype was more likely to be basal tumors and was significantly associated with poorer prognosis even after controlling for age, gender, histological grade, clinical stage, and molecular phenotypes. Additionally, mixed tumors harbored a higher frequency of RB1 and LRP1B copy number deletion compared to quiescent tumors (25.7% vs. 12.7 and 27.9% vs. 10.2%, respectively, both adjusted P value< 0.05). Furthermore, aberrant PIK3CA expression level was significantly correlated with those of glycolytic and cholesterogenic genes. The quiescent subtype was associated with lower stemness indices and lower signature scores for gene sets involved in genomic instability, including DNA replication, DNA damage repair, mismatch repair, and homologous recombination genes. Moreover, quiescent tumors exhibited lower expression levels of pyruvate dehydrogenase kinases 1-3 (PDK1-3) than the other subtypes. In addition, distinct immune cell infiltration patterns were observed across the 4 metabolic subtypes, with greater infiltration of M0/M2 macrophages observed in glycolytic and mixed subtypes. However, no significant difference in immunotherapy response was observed across the 4 metabolic subtypes.
This study proposed a new metabolic subtyping method for BCa based on genes involved in glycolysis and cholesterol synthesis pathways. Our findings may provide novel insight for the development of personalized subtype-specific treatment strategies targeting metabolic vulnerabilities.
Bladder cancer (BCa) is one of the most common tumors and the thirteenth leading cause of cancer-related deaths worldwide . Globally, approximately 573,000 patients were diagnosed with BCa in 2020, with 212,000 BCa-related deaths in that year, posing an enormous threat to human health . In recent times, BCa treatment has been revolutionized by emerging therapies such as immunotherapy and molecular-targeted therapy; however, the 5-year survival rate for muscle-invasive BCa remains unsatisfactory [3, 4]. Therefore, continuous understanding of tumor subtyping is desirable to improve the prognostic stratification of BCa to achieve personalized treatment.
Oncogenic metabolic reprogramming is a major hallmark of cancers and allows cancer cells to survive and thrive in harsh conditions [5, 6]. Tumor cells, by an effect known as the Warburg effect, can shift glucose metabolism toward aerobic glycolysis, providing cancer cells with energy and biosynthetic raw materials to promote tumor growth, invasion, and metastasis [5,6,7,8]. Metabolic reprogramming also plays a pivotal role in maintaining genomic instability and stemness in cancer cells to allow for self-expansion and resistance to chemotherapy. In addition, glycolytic reprogramming modifies the tumor microenvironment (TME) into a hypoxic, acidic, and nutritionally deficient environment that facilitates cancer cell growth and inhibits immune cell function [5,6,7,8]. Therefore, targeting metabolic vulnerability is a promising strategy for cancer therapy.
Pyruvate is the terminal product of glycolysis and servers as a precursor for different biosynthetic pathways. Mitochondrial pyruvate complex (MPC) comprised of pyruvate carriers 1 and 2 (MPC1/MPC2) is the entry point for pyruvate to the mitochondrial matrix for oxidative metabolism, and MPC deficiency contributes to tumor initiation and progression by enhancing glycolysis [9, 10]. In addition to MPC, the pyruvate dehydrogenase complex (PDC) also serves as a gatekeeper in maintaining the balance between anaerobic and aerobic glucose metabolism by catalyzing the conversion of pyruvate to acetyl-CoA for entry into the tricarboxylic acid cycle (TAC) . PDC activity is tightly regulated by pyruvate dehydrogenase kinases (PDKs, isoform 1-4), which can phosphorylate and inactivate the PDC, resulting in attenuated pyruvate oxidative metabolism and increased glycolysis . PDK inhibition has been reported to suppress tumor growth and is a promising therapeutic target for several diseases, including cancers [13, 14].
Recent studies have revealed the remarkable cancer prognosis-determining potential of metabolic heterogeneity . For instance, using glycolytic and cholesterogenic genes, Karasinska et al. classified pancreatic cancer into 4 distinct metabolic phenotypes, which were related to patient survival, and established molecular subtypes . BCa is a heterogeneous malignancy , and it is unclear whether it can be stratified into different subtypes through heterogeneity in distinct metabolic pathways to enhance personalized therapy.
In this study, we stratified BCa into 4 distinct metabolic subtypes based on the expression levels of glycolytic and cholesterogenic genes. We aim to clarify the prognostic value of heterogeneity in glycolysis and cholesterol synthesis, determine its association with genomic instability, stemness, and the immune microenvironment in BCa, and provide a research basis for further personalized treatment.
Materials and methods
Study datasets and participants
The Cancer Genome Atlas (TCGA) datasets, including RNA-seq expression, single nucleotide variants (SNV) /indels, copy number variation (CNV), and corresponding clinical profiles were downloaded from the GDC portal (https://portal.gdc.cancer.gov/). The GSE13507 microarray dataset was downloaded from the GEO portal (https://www.ncbi.nlm.nih.gov/geo/). The IMvigor210 study, which evaluated the efficacy and safety of PD-L1 inhibitors in locally advanced or advanced urothelial cancer [17, 18], was obtained from http://research-pub.gene.com/IMvigor210CoreBiologies/. A total of 760 primary BCa samples with survival data (TCGA, n = 400; GSE13507, n = 165; IMvigor210, n = 195) were used for this study. The gene expression values for RNA-seq were transformed into TPM (transcripts per kilobase million). RNA-seq and microarray gene expression data were log2 transformed for analysis. Batch effects caused by non-biological technical biases were corrected using the “ComBat” algorithm of “sva” R package. Principal component analysis (PCA) was used to evaluate the batch effect between samples before and after correction. As shown in Supplementary Fig. S1, the PCA confirmed a reduction in batch effects after normalization. The detailed clinical-pathological features of the 3 datasets are shown in Supplementary Table S1.
To stratify BCa based on the relative expression levels of genes involved in glycolysis and cholesterol biosynthesis, genes belonging to Reactome gene sets, ‘glycolysis’ (n = 72), and ‘cholesterol biosynthesis’ (n = 25), were extracted from ‘MsigDB’ (supplementary Table S2). The batch effect-corrected expression values for these genes were standardized by Z-score, and then subjected to consensus clustering (parameters: rps = 1000, pItem = 0.8, pFeature = 1, clusterAlg = hc, distance = euclidean) using the ‘ConsensusClusterPlus’ R package to identify co-expressed glycolysis and cholesterol synthesis genes. The number of clusters was determined according to the criteria of consensus Cumulative Distribution Function (CDF) and the relative change in area under the CDF curve. Samples were divided into 4 metabolic subtypes based on the median expression values of co-expressed glycolytic and cholesterogenic genes i.e., quiescent (glycolytic ≤0, cholesterogenic ≤0); glycolytic (glycolytic > 0, cholesterogenic ≤0); cholesterogenic (glycolytic ≤0, cholesterogenic > 0); mixed (glycolytic > 0, cholesterogenic > 0) subtypes.
Molecular phenotype classification
Bladder cancer samples were classified into five molecular phenotypes (basal/, luminal, luminal-infiltrated, luminal-infiltration, and neuronal) using the TCGA molecular classifier , or into three molecular phenotypes (basal, luminal, and P53-like) using the MD Anderson molecular classifier , according to the “BLCAsubtype” R package.
Tumor microenvironment immune cell infiltration
The CIBERSORT algorithm was used to estimate the relative abundance of 22 tumor-infiltrating immune cell types with the LM22 reference gene signature and 1000 permutations based on TCGA RNA-seq data . The Kruskal-Wallis test was used to compare the degree of immune cell infiltration between the different metabolic subtypes.
Single-sample gene-set enrichment analysis and stemness index
We performed single-sample gene-set enrichment analysis (ssGSEA) on TCGA samples to estimate the enrichment scores of gene sets involved in maintaining genomic instability . The gene sets included (1) ‘DNA replication,’ (2) ‘mismatch repair,’ (3) ‘base excision repair,’ (4) ‘nucleotide excision repair,’ (5) ‘DNA damage repair,’ (6) ‘homologous recombination,’ and (7) ‘cell cycle’ genes. The stemness index, which was estimated using RNA expression (all available gene sets), was obtained from the UCSC Xena Pan-Cancer Atlas Hub (https://xenabrowser.net/). ANOVA was performed to compare differences in these signature scores between the 4 metabolic subtypes.
Pearson correlation analysis was performed to identify genes significantly correlated with PDK1-3 expression. Genes with |R| > 0.3 and P value< 0.001 were considered to be significant. Then, genes significantly correlated with all PDK1, PDK2, and PDK3 were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the “clusterProfiler” R package.
R (version 4.0.2) was used for statistical analyses. Kaplan-Meier plots with log-rank test were used to test differences in overall survival using the ‘survival’ and ‘survminer’ tools in the R software package. Cox regression was used to evaluate differences in overall survival after adjusting for potential confounders such as age, gender, histological grade, and cancer stage. The Kruskal-Wallis test, Chi-square test, ANOVA, or Fisher exact test was used for between-group comparisons where appropriate. Values of P < 0.05 were considered statistically significant.
Distinction of the 4 bladder cancer metabolic subtypes based on the expression levels of glycolysis and cholesterol synthesis genes
As shown in Fig. 1A, ten clusters of co-expressed glycolysis- and cholesterol-related genes were obtained by consensus clustering based on the consensus CDF and change in area under the CDF curve. We identified 2 groups of robustly co-expressed glycolysis (n = 12, F4 cluster) and cholesterol synthesis (n = 8, F2 cluster) genes which were used for metabolic subtyping (Fig. 1A, Supplementary Table S3). Patients were assigned one of the 4 metabolic subtypes, i.e., quiescent, glycolytic, cholesterogenic, and mixed (Fig. 1B), based on the median expression levels of glycolytic and cholesterogenic genes. The largest proportion of samples exhibited the quiescent subtype (31.1%), followed by the mixed (24.3%), cholesterogenic (22.5%), and glycolytic (22.2%) subtypes. The expression levels of these co-expressed glycolytic and cholesterogenic genes were shown in Fig. 1C.
Clinical significances of the 4 metabolic subtypes
As shown in Fig. 2A, the 4 metabolic subtypes demonstrated significant differences in overall survival (log-rank P value = 0.012). The quiescent subtype exhibited the best survival, while the mixed exhibited shorter survival. The survival benefit of quiescent tumors in comparison to mixed tumors was observed in muscle-invasive BCa (overall log-rank P = 0.011); however, the difference was marginally significant in non-muscle-invasive BCa (log-rank P = 0.09), probably due to the limited sample size. (Fig. 2B-C). In addition, significant differences in the histological grade were observed across the 4 metabolic subtypes (P = 8.95e-5) and quiescent tumors tended to exhibit lower histological grade and less advanced pathological stage, although the latter did not reach statistical significance (P = 0.087) (Fig. 2D-E). Moreover, significant differences in molecular phenotypes were also observed across the 4 metabolic subtypes (P < 2.2e-16). Using the MDA molecular classifier, basal tumors were more common in the mixed and glycolytic subtypes (Fig. 2F). Similar findings were observed using the TCGA molecular classifier. Luminal (luminal infiltrated, luminal papillary, and luminal) tumors were more common in the quiescent and cholesterogenic subtypes, while basal-squamous tumors were more common in the mixed and glycolytic subtypes (Fig. 2G). Furthermore, multivariable cox regression revealed that mixed tumors remained an independent predictor for poorer prognosis after controlling for age, gender, histological grade, clinical stage, and molecular phenotypes (Fig. 2H-I and supplementary Table S4). These data indicate that tumors with higher rates of glycolysis and cholesterol synthesis may be more aggressive than tumors with a quiescent subtype, and the metabolic subtype based on glycolytic and cholesterogenic genes may be a promising classifier for prognostic stratification of BCa.
Genetic alterations of glycolytic and cholesterogenic genes across the 4 metabolic subtypes
To determine the genetic alteration events associated with metabolic subtypes, we investigated the frequency of SNVs, indels, and CNVs of the 20 co-expressed glycolytic and cholesterogenic genes in the TCGA cohort. As shown in Fig. 3A, CNVs commonly occurred while SNV and indel mutations were rare. The SNV and indel mutation frequencies of the genes were comparable across the four metabolic subtypes (Supplementary Table S5). Of note, CNV gain was more frequently observed than CNV loss across the 4 metabolic subtypes (Fig. 3B), consistent with the upregulated expression levels of these genes. Significant differences in frequencies of CNV alterations were observed in glycolytic genes ENO1 (p = 0.009), PFKP (P = 0.004), and cholesterogenic genes HMGCR (P = 0.016), HMGCS1 (P = 0.003), and IDI1 (P = 0.008) across the 4 metabolic subtypes (Fig. 3A and Supplementary Table S6), indicating that aberrant expression of these genes may be related to the development of BCa.
Association between metabolic subtypes and tumor genomic alterations
Genomic alterations are capable of driving tumor metabolic reprogramming [22, 23]. For instance, oncogenic PIK3CA mutations have been reported to reprogram metabolism in a variety of cancers including BCa [23, 24]. In this study, we investigated the genomic alteration (SNV, indel, and CNV) landscapes associated with metabolic subtypes in the TCGA dataset. Of the 30 most frequently altered genes in BCa, the highest frequencies of alteration (SNVs, indels, and CNVs) were observed in TTN (58%), TP53 (45%), and MUC16 (39.0%). Significant differences in the frequencies of alteration in RB1 (P < 0.001), ARID1A (P = 0.008), LRP1B (P = 0.018), CSMD3 (P = 0.011), and PIK3CA (P = 0.018) were observed across the 4 metabolic subtypes (Fig. 4A, Supplementary Table S7). Of note, CNV loss was frequently observed in RB1 and LRP1B, both of which exhibited the lowest frequencies in the quiescent subtype (Fig. 4B-C). Furthermore, correlation analysis revealed a significant correlation between the expression levels of these genes and the median expression of glycolytic and cholesterogenic genes (Supplementary Table S8). Of these genes, PIK3CA expression exhibited the strongest correlation with the expression of glycolytic and cholesterogenic genes (R = 0.43 and 0.52, respectively, and both P values were < 0.001) (Fig. 4D-E). These findings are compatible with the notion that PIK3CA drives tumor metabolic reprogramming and promotes tumor progression . In summary, this study revealed a distinct genomic alteration landscape across the 4 metabolic subtypes, suggesting the importance of crosstalk between genome instability and metabolic reprogramming in the development of BCa.
Association between genomic instability or cancer stemness index and metabolic subtypes
To further investigate the relationship between genome instability and metabolic subtypes, the signature scores of gene sets involved in genomic instability, including DNA replication, mismatch repair, base excision repair, nuclear excision repair, DNA damage repair, homolog recombination, and cell cycle genes, were compared. As shown in Fig. 5A-G, significant differences were observed between the 4 metabolic subtypes, with the lowest scores generally observed in the quiescent and mixed subtypes, indicating a close relationship of glycolysis and cholesterol biosynthesis with genomic instability that drives tumorigenesis and therapeutic resistance. Furthermore, using the mRNA-based stemness indices derived from TCGA RNA-seq data, we observed the lowest and highest stemness indices in the quiescent and mixed subtypes, respectively (Fig. 5H). In summary, these data indicate a close relationship between metabolic reprogramming and genomic instability and identify glycolysis and cholesterol synthesis as potential targets for controlling BCa.
TME infiltrating cells and metabolic subtypes
Infiltrating immune cells are an important component of the TME and play a critical role in carcinogenesis and tumor therapeutic response . In this study, using the CIBERSORT algorithm, distinct infiltration patterns were observed across the 4 metabolic subtypes in the TCGA dataset (Fig. 6A-B). Among the 22 immune cell types, naive B cells, plasma cells, CD4 memory-activated T cells, regulatory T cells (Tregs), resting NK cells, monocytes, macrophages (M0/M2), and resting mast cells demonstrated significant differences in infiltrating abundances across the 4 subtypes. Of note, the quiescent and cholesterogenic subtypes exhibited lower M0/M2 macrophage infiltration than the glycolytic and mixed subtypes, and high M0/M2 macrophage levels were associated with poor prognosis (Supplementary Fig. S2). These data indicate that while both glycolysis and cholesterol synthesis contribute to shaping the TME to facilitate tumor growth, glycolysis may have a more significant effect than cholesterol synthesis. We also investigated the association between metabolic subtype and immunotherapeutic response in the IMvigor210 study. However, no significant differences in response rate and survival benefit were observed between the 4 metabolic subtypes (Fig. 6C-D). Further studies with larger sample sizes are needed to confirm these findings.
Association between MPC1/2 or PDK1-4 alteration and metabolic subtypes
Because pyruvate is the terminal product of anaerobic glycolysis and acts as a precursor for different biosynthetic pathways, including cholesterol biosynthesis, we further investigated MPC and PDKs, both are critical players involved in pyruvate processing. The mitochondrial pyruvate carrier (MPC) complex, which consists of MPC1 and MPC2, is required for efficient glucose production, and decreased MPC activity in tumors enhances glycolytic activity, resulting in tumor progression . In TCGA BCa samples, MPC1 and MPC2 mutations were rare, with only 1 mutation affecting MPC1 observed in a glycolytic subtype patient. In contrast, CNVs were common, with the majority of CNVs being deletions in MPC1 and amplifications in MPC2, although no significant differences in MPC1/2 CNV frequencies were observed between the 4 metabolic subtypes (Fig. 7A). Likewise, no significant differences in the expression levels of MPC1 were observed between the subtypes; however, the expression levels of MPC2 were lower in the glycolytic and mixed phenotypes as compared to the cholesterogenic subtype (Fig. 7B-C).
Similarly, SNV and indel mutations were rare but CNVs were common in PDK1-4; however, no significant differences in rates of genetic alteration (CNV or mutation) were observed between the 4 subtypes (Fig. 7A). Nevertheless, PDK gene expression levels were significantly different between the 4 metabolic subtypes. Overall, quiescent subtypes exhibited the lowest PDK1-3 expression levels (Fig. 7D-G). Because PDK1-3 expression levels demonstrated the most significant correlations with metabolic subtypes and were significantly lower in quiescent tumors compared to mixed tumors, we further performed Pearson correlation analysis to identify genes significantly correlated with PDK1, PDK2, and PDK3 expression. As shown in Fig. 7H-J, a total of 1374 shared genes were identified to be significantly correlated with PDK1, PDK2, and PDK3 (all Pearson correlation R > 0.3, and P value< 0.001), and these genes were significantly enriched in spliceosome, RNA transport, nucleotide excision repair, DNA replication, RNA splicing, et al., in consistence with our findings that metabolic subtypes were tightly associated with genomic instability. In summary, these data suggest that decreased PDC activity resulting from increased PDK mRNA levels may contribute to the development of the malignant features associated with metabolic subtypes, presenting PDKs as potential targets for controlling BCa through the manipulation of metabolic vulnerability.
In this study, based on the expression levels of genes involved in glycolysis and cholesterol synthesis, we identified 4 metabolic subtypes, which demonstrated distinct clinical and molecular characteristics, and different TME immune cell infiltration patterns, although no significant differences were observed in immunotherapy response rate. In summary, our study unveils the importance of metabolic reprogramming in BCa and presents metabolic heterogeneity-based subtyping as a potential prognosis biomarker for personalized therapy.
Glucose metabolic reprogramming is essential for tumor growth and therapeutic resistance . Studies have shown that glycolysis is closely related to BCa development. For instance, upregulation of pyruvate kinase M2 is closely related to tumor growth and chemo-resistance and serves as a potential tumor marker for BCa monitoring [27, 28]. Lactic acid dehydrogenase (LDHA) is a key enzyme in glycolysis and its upregulation in BCa promotes glycolysis, thereby facilitating tumor growth and immune evasion [29, 30]. Aside from glycolysis, increasing evidence shows that cholesterol metabolites also play critical roles in cancer development [31, 32]. Increased cholesterol biosynthesis is a hallmark of many cancers, promoting cancer cell growth and immune evasion by activating cellular signalings such as sonic hedgehog, Notch and receptor tyrosine kinases, and LXR-α signaling . Based on the expression levels of glycolytic and cholesterogenic genes, Karasinska et al. identified 4 distinct metabolic phenotypes with remarkable prognostic significance . As with previous studies, this study also revealed a close relationship between metabolic subtypes and BCa prognosis. Notably, the quiescent and mixed subtypes were associated with the best and worst outcomes, respectively. These results indicate that glycolysis and cholesterol synthesis might act synergistically to accelerate tumor progression in BCa.
It has been well recognized that the molecular characteristics of BCa are associated with prognosis and therapeutic responses. Compared to luminal tumors, basal tumors demonstrated unfavorable survival and suboptimal therapeutic response [19, 20]. In this study, basal tumors were more common in the mixed subtype using both TCGA molecular classifier and MDA classifier, consistent with the unfavorable survival associated with mixed tumors, suggesting a close relationship between metabolic reprogramming and molecular phenotypes. More importantly, metabolic subtype remained as a significant predictor for overall survival after controlling for major confounders including molecular phenotypes. Taken together, our study highlights the prognostic value of metabolic subtypes in guiding personalized therapy.
Genomic instability resulting from mutations in DNA repair genes is a hallmark of most cancers and plays a central role in tumor initiation and progression [34, 35]. Recent studies have revealed the close relationship between metabolic reprogramming and cancer genomic instability . For instance, glycolysis was found to contribute to DNA metabolism by providing metabolites essential for the biosynthesis of nucleotides. Some glycolytic products (like L- and D-lactate) and key glycolytic enzymes (like PGAM1 and PKM2) are involved in DNA damage repair [37,38,39]. Furthermore, the importance of cholesterol biosynthesis in maintaining genome instability has also been reported [40, 41]. Previous studies have reported that genetic defects of the nucleotide excision repair pathway, including ERCC1 deficiency, result in the suppression of cholesterol biosynthesis . In line with previous reports, our study also revealed a close relationship between glycolysis and cholesterol biosynthesis and genomic stability. Of the 4 metabolic subtypes, the quiescent subtype exhibited the lowest activities in mismatch repair, base excision repair, nucleotide excision repair, DNA damage repair, and DNA replication, while the mixed subtype exhibited relatively higher activities than the other subtypes. Moreover, we observed significant differences in genomic alteration patterns between the 4 metabolic subgroups. RB1 is a well-known tumor suppressor gene, and its deletion can enhance glycolytic metabolism and driving tumor progression . A recent study also identified LRP1B, which encodes low-density lipoprotein receptor-related protein 1B, as a novel tumor suppressor, and LRP1B deletion was associated with chemotherapy resistance in ovarian cancer . In this study, although the quiescent tumors harbored similarly high mutation rates in most of the interested genes, consistent with findings reported in a previous study , it is worthy of note that the quiescent subtype exhibited lower frequency of CNV loss in RB1 and LRP1B than mixed subtypes, in line with the survival benefit of the quiescent subtype. In addition, we observed frequent mutations in PIK3CA, which were positively associated with glycolysis and cholesterol synthesis. In consonance with our findings, previous studies have reported that PIK3CA mutation promotes tumor progression partly by enhancing glycolysis . In summary, our findings suggest a close relationship of glycolysis and cholesterol synthesis with genomic instability and present glycolytic and cholesterogenic metabolism targeting as a potential therapeutic strategy for BCa treatment.
Accumulating evidence revealed that metabolic reprograming can contribute to tumor progression by creating a hypoxic, acidic, and nutritionally deficient TME [5,6,7,8]. However, the TME cell-infiltrating characteristics of distinct metabolic subtypes remain unclear. In this study, the 4 metabolic subtypes exhibited a distinct immune cell infiltration pattern, demonstrating that the glycolytic and cholesterogenic reprogramming is critical in shaping different TME landscapes. Of note, the quiescent subtype exhibited significantly lower M0/M2 macrophage levels than the other subtypes. A recent study also demonstrated that macrophage-promoted tumor growth by regulating tumor cell metabolism, in support of our findings . Furthermore, macrophages have been proved to be critical players in driving cancer cell immune evasion and are closely related to poor prognosis . Therefore, a comprehensive assessment of the metabolic patterns may enhance our understanding of TME cell-infiltrating characteristics. In this study, we speculated that metabolic subtype could predict immunotherapeutic response, however, we did not observe significant differences in the response to PD-L1 immunotherapy between the 4 subtypes in IMvigor210 study. Therefore, further studies are warranted to evaluate the roles of glycolysis and cholesterol synthesis in immunoregulation in BCa. Furthermore, recent studies have reported that as with cancer cells, TME immune cells also undergo metabolic reprogramming to facilitate tumor growth and immune evasion [48, 49]. In this study, although we demonstrated a distinct TME cell-infiltration pattern associated with heterogeneity in glycolysis and cholesterol synthesis, we did not investigate the metabolic alteration in TME immune cells, which may also contribute to the alteration in the transcriptome profiles of tumor tissues. Further studies, such as studies with single-cell RNA-sequencing, are needed to investigate the crosstalk between metabolic alterations in tumors and TME cells during tumor initiation and progression.
The MPC deficiency has been linked to tumorigenesis by enhancing glycolysis and may serve as a potential target of anticancer therapy by manipulating glycolytic activity [9, 10]. In line with previous reports on other cancer types , this study revealed that in BCa, MPC1 was frequently deleted, while MPC2 was mostly amplified. Although no significant differences in MPC1 expression levels were observed between the 4 metabolic subtypes, the decreased MPC2 expression in the glycolytic and mixed subtypes suggests that MPC2 may contribute to mitochondrial pyruvate uptake in BCa. In addition to MPC, the PDC also plays a pivotal role in regulating energy homeostasis [50, 51]. Studies have suggested that metabolic reprogramming in cancer cells is associated with PDC inhibition due to phosphorylation of its E1a subunit by PDKs [14, 52, 53]. Thus, inhibition of PDK has been recognized as an attractive strategy in anticancer therapy [54, 55]. While the roles and mechanisms of PDK1-3 in BCa remain unknown, enhanced PDK4 expression in BCa has been reported in previous studies, and PDK4 inhibitors were found to suppress BCa cell invasiveness and to potentiate cisplatin-induced cell death [56, 57]. In this study, we observed significantly lower PDK1-3 expression levels in the quiescent subtype as compared to the other subtypes, consistent with the notion that PDKs contribute to tumorigenesis by inhibiting PDC activity [54, 55]. Furthermore, GO and KEGG enrichment analysis revealed that PDK1-3 was correlational with the critical biological process involved in tumorigenesis, including nucleotide excision repair, DNA replication, RNA splicing, etc., presenting PDK1-3 as potential therapeutic targets for BCa. However, the decreased expression levels of PDK4 in the mixed and glycolytic subtypes did not match our inferences and was inconsistent with other reports [56, 57]. Further studies are needed to investigate the roles played by PDKs in BCa development.
Our study had some limitations. Firstly, we only addressed the intertumoral glycolytic and cholesterogenic heterogeneity, however, the intratumoral metabolic heterogeneity and the comprehensive landscape of tumor metabolism remained uninvestigated. Secondly, the current study is largely based on correlation analysis and there is a lack of experimental validation; therefore, the findings of this study should be interpreted with caution. Thirdly, although we combined 3 cohorts to obtain a large sample size, the sample size in the subgroup analysis of non-muscle-invasive BCa was still small.
Overall, this study identified a metabolic heterogeneity-based classifier with distinct molecular and immune characteristics, and predicted outcomes in BCa, providing novel insights for the development of personalized therapeutic strategies targeting metabolic vulnerabilities.
Availability of data and materials
The current study analyzed publicly available datasets which can be found in https://xena.ucsc.edu/welcome-to-ucsc-xena (TCGA); https://www.ncbi.nlm.nih.gov/gds/?term=GSE13507 (GSE13507); http://research-pub.gene.com/IMvigor210CoreBiologies/ (IMvigor210).
Mitochondrial pyruvate complex
Pyruvate dehydrogenase complex
Tricarboxylic acid cycle
Pyruvate dehydrogenase kinases
Single-sample gene-set enrichment analysis
Principal component analysis
Kyoto Encyclopedia of Genes and Genomes
Richters A, Aben KKH, Kiemeney L. The global burden of urinary bladder cancer: an update. World J Urol. 2020;38(8):1895–904.
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49.
Tran L, Xiao JF, Agarwal N, Duex JE, Theodorescu D. Advances in bladder cancer biology and therapy. Nat Rev Cancer. 2021;21(2):104–21.
Berdik C. Unlocking bladder cancer. Nature. 2017;551(7679):S34–5.
Ward PS, Thompson CB. Metabolic reprogramming: a cancer hallmark even warburg did not anticipate. Cancer Cell. 2012;21(3):297–308.
Faubert B, Solmonson A, DeBerardinis RJ. Metabolic reprogramming and cancer progression. Science. 2020;368(6487):eaaw5473.
Peng X, Chen Z, Farshidfar F, Xu X, Lorenzi PL, Wang Y, et al. Molecular characterization and clinical relevance of metabolic expression subtypes in human cancers. Cell Rep. 2018;23(1):255–269 e254.
DeBerardinis RJ, Lum JJ, Hatzivassiliou G, Thompson CB. The biology of cancer: metabolic reprogramming fuels cell growth and proliferation. Cell Metab. 2008;7(1):11–20.
Herzig S, Raemy E, Montessuit S, Veuthey JL, Zamboni N, Westermann B, et al. Identification and functional expression of the mitochondrial pyruvate carrier. Science. 2012;337(6090):93–6.
Schell JC, Olson KA, Jiang L, Hawkins AJ, Van Vranken JG, Xie J, et al. A role for the mitochondrial pyruvate carrier as a repressor of the Warburg effect and colon cancer cell growth. Mol Cell. 2014;56(3):400–13.
Park S, Jeon JH, Min BK, Ha CM, Thoudam T, Park BY, et al. Role of the pyruvate dehydrogenase complex in metabolic remodeling: differential pyruvate dehydrogenase complex functions in metabolism. Diabetes Metab J. 2018;42(4):270–81.
Zhang S, Hulver MW, McMillan RP, Cline MA, Gilbert ER. The pivotal role of pyruvate dehydrogenase kinases in metabolic flexibility. Nutr Metab (Lond). 2014;11(1):10.
Anwar S, Shamsi A, Mohammad T, Islam A, Hassan MI. Targeting pyruvate dehydrogenase kinase signaling in the development of effective cancer therapy. Biochim Biophys Acta Rev Cancer. 2021;1876:188568.
Stacpoole PW. Therapeutic targeting of the pyruvate dehydrogenase complex/pyruvate dehydrogenase kinase (PDC/PDK) axis in cancer. J Natl Cancer Inst. 2017;109(11):djx071.
Karasinska JM, Topham JT, Kalloger SE, Jang GH, Denroche RE, Culibrk L, et al. Altered gene expression along the glycolysis-cholesterol synthesis axis is associated with outcome in pancreatic cancer. Clin Cancer Res. 2020;26(1):135–46.
Meeks JJ, Al-Ahmadie H, Faltas BM, Taylor JA 3rd, Flaig TW, DeGraff DJ, et al. Genomic heterogeneity in bladder cancer: challenges and possible solutions to improve outcomes. Nat Rev Urol. 2020;17(5):259–70.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.
Rosenberg JE, Hoffman-Censits J, Powles T, van der Heijden MS, Balar AV, Necchi A, et al. Atezolizumab in patients with locally advanced and metastatic urothelial carcinoma who have progressed following treatment with platinum-based chemotherapy: a single-arm, multicentre, phase 2 trial. Lancet. 2016;387(10031):1909–20.
Robertson AG, Kim J, Al-Ahmadie H, Bellmunt J, Guo G, Cherniack AD, et al. Comprehensive molecular characterization of muscle-invasive bladder cancer. Cell. 2017;171(3):540–556 e525.
Choi W, Porten S, Kim S, Willis D, Plimack ER, Hoffman-Censits J, et al. Identification of distinct basal and luminal subtypes of muscle-invasive bladder cancer with different sensitivities to frontline chemotherapy. Cancer Cell. 2014;25(2):152–65.
Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med. 2015;21(8):938–45.
Zhang C, Liu J, Liang Y, Wu R, Zhao Y, Hong X, et al. Tumour-associated mutant p53 drives the Warburg effect. Nat Commun. 2013;4:2935.
Hao Y, Samuels Y, Li Q, Krokowski D, Guan BJ, Wang C, et al. Oncogenic PIK3CA mutations reprogram glutamine metabolism in colorectal cancer. Nat Commun. 2016;7:11971.
Samuels Y, Waldman T. Oncogenic mutations of PIK3CA in human cancers. Curr Top Microbiol Immunol. 2010;347:21–41.
Whiteside TL. The tumor microenvironment and its role in promoting tumor growth. Oncogene. 2008;27(45):5904–12.
Afonso J, Santos LL, Longatto-Filho A, Baltazar F. Competitive glucose metabolism as a target to boost bladder cancer immunotherapy. Nat Rev Urol. 2020;17(2):77–106.
Avayu O, Almeida E, Prior Y, Ellenbogen T. Composite functional metasurfaces for multispectral achromatic optics. Nat Commun. 2017;8:14992.
Liu W, Woolbright BL, Pirani K, Didde R, Abbott E, Kaushik G, et al. Tumor M2-PK: a novel urine marker of bladder cancer. PLoS One. 2019;14(6):e0218737.
Yuan D, Zheng S, Wang L, Li J, Yang J, Wang B, et al. MiR-200c inhibits bladder cancer progression by targeting lactate dehydrogenase A. Oncotarget. 2017;8(40):67663–9.
Mishra D, Banerjee D. Lactate dehydrogenases as metabolic links between tumor and stroma in the tumor microenvironment. Cancers (Basel). 2019;11(6):750.
Ding X, Zhang W, Li S, Yang H. The role of cholesterol metabolism in cancer. Am J Cancer Res. 2019;9(2):219–27.
Kuzu OF, Noory MA, Robertson GP. The role of cholesterol in cancer. Cancer Res. 2016;76(8):2063–70.
Huang B, Song BL, Xu C. Cholesterol metabolism in cancer: mechanisms and therapeutic opportunities. Nat Metab. 2020;2(2):132–41.
Negrini S, Gorgoulis VG, Halazonetis TD. Genomic instability--an evolving hallmark of cancer. Nat Rev Mol Cell Biol. 2010;11(3):220–8.
Shen Z. Genomic instability and cancer: an introduction. J Mol Cell Biol. 2011;3(1):1–3.
Sobanski T, Rose M, Suraweera A, O'Byrne K, Richard DJ, Bolderson E. Cell metabolism and DNA repair pathways: implications for cancer therapy. Front Cell Dev Biol. 2021;9:633305.
Qu J, Sun W, Zhong J, Lv H, Zhu M, Xu J, et al. Phosphoglycerate mutase 1 regulates dNTP pool and promotes homologous recombination repair in cancer cells. J Cell Biol. 2017;216(2):409–24.
Wagner W, Ciszewski WM, Kania KD. L- and D-lactate enhance DNA repair and modulate the resistance of cervical carcinoma cells to anticancer drugs via histone deacetylase inhibition and hydroxycarboxylic acid receptor 1 activation. Cell Commun Signal. 2015;13:36.
Yang W, Xia Y, Hawke D, Li X, Liang J, Xing D, et al. PKM2 phosphorylates histone H3 and promotes gene transcription and tumorigenesis. Cell. 2014;158(5):1210.
Zhang Y, Liu Y, Duan J, Wang H, Qiao K, Wang J. Cholesterol depletion sensitizes gallbladder cancer to cisplatin by impairing DNA damage response. Cell Cycle. 2019;18(23):3337–50.
Enriquez-Cortina C, Bello-Monroy O, Rosales-Cruz P, Souza V, Miranda RU, Toledo-Perez R, et al. Cholesterol overload in the liver aggravates oxidative stress-mediated DNA damage and accelerates hepatocarcinogenesis. Oncotarget. 2017;8(61):104136–48.
Smith SC, Robinson AR, Niedernhofer LJ, Hetman M. Downregulation of cholesterol biosynthesis genes in the forebrain of ERCC1-deficient mice. Neurobiol Dis. 2012;45(3):1136–44.
Conroy LR, Dougherty S, Kruer T, Metcalf S, Lorkiewicz P, He L, et al. Loss of Rb1 enhances glycolytic metabolism in Kras-driven lung tumors in vivo. Cancers (Basel). 2020;12(1):237.
Cowin PA, George J, Fereday S, Loehrer E, Van Loo P, Cullinane C, et al. LRP1B deletion in high-grade serous ovarian cancers is associated with acquired chemotherapy resistance to liposomal doxorubicin. Cancer Res. 2012;72(16):4060–73.
Jiang W, He T, Liu S, Zheng Y, Xiang L, Pei X, et al. The PIK3CA E542K and E545K mutations promote glycolysis and proliferation via induction of the beta-catenin/SIRT3 signaling pathway in cervical cancer. J Hematol Oncol. 2018;11(1):139.
Zhang Y, Yu G, Chu H, Wang X, Xiong L, Cai G, et al. Macrophage-associated PGK1 phosphorylation promotes aerobic glycolysis and tumorigenesis. Mol Cell. 2018;71(2):201–215 e207.
Noy R, Pollard JW. Tumor-associated macrophages: from mechanisms to therapy. Immunity. 2014;41(1):49–61.
Wu S, Kuang H, Ke J, Pi M, Yang DH. Metabolic reprogramming induces immune cell dysfunction in the tumor microenvironment of multiple myeloma. Front Oncol. 2020;10:591342.
Sasidharan Nair V, Saleh R, Toor SM, Cyprian FS, Elkord E. Metabolic reprogramming of T regulatory cells in the hypoxic tumor microenvironment. Cancer Immunol Immunother. 2021;70:2103.
Sutendra G, Kinnaird A, Dromparis P, Paulin R, Stenson TH, Haromy A, et al. A nuclear pyruvate dehydrogenase complex is important for the generation of acetyl-CoA and histone acetylation. Cell. 2014;158(1):84–97.
McFate T, Mohyeldin A, Lu H, Thakar J, Henriques J, Halim ND, et al. Pyruvate dehydrogenase complex activity controls metabolic and malignant phenotype in cancer cells. J Biol Chem. 2008;283(33):22700–8.
Atas E, Oberhuber M, Kenner L. The implications of PDK1-4 on tumor energy metabolism, aggressiveness and therapy resistance. Front Oncol. 2020;10:583217.
Kolobova E, Tuganova A, Boulatnikov I, Popov KM. Regulation of pyruvate dehydrogenase activity through phosphorylation at multiple sites. Biochem J. 2001;358(Pt 1):69–77.
Sradhanjali S, Reddy MM. Inhibition of pyruvate dehydrogenase kinase as a therapeutic strategy against cancer. Curr Top Med Chem. 2018;18(6):444–53.
Woolbright BL, Rajendran G, Harris RA, Taylor JA 3rd. Metabolic flexibility in cancer: targeting the pyruvate dehydrogenase kinase:pyruvate dehydrogenase axis. Mol Cancer Ther. 2019;18(10):1673–81.
Woolbright BL, Choudhary D, Mikhalyuk A, Trammel C, Shanmugam S, Abbott E, et al. The role of pyruvate dehydrogenase kinase-4 (PDK4) in bladder cancer and chemoresistance. Mol Cancer Ther. 2018;17(9):2004–12.
Kim CJ, Terado T, Tambe Y, Mukaisho KI, Kageyama S, Kawauchi A, et al. Cryptotanshinone, a novel PDK 4 inhibitor, suppresses bladder cancer cell invasiveness via the mTOR/betacatenin/Ncadherin axis. Int J Oncol. 2021;59(1):40.
YZ acknowledged Prof. Haiping Long for his continuous supports in her research.
This work was supported by the grants from National Natural Science Foundation of China (No. 81802551, No. 81900688, and No. 81800590), Chinese Postdoctoral Science Foundation (2020 M672593), the Medical Research Foundation of Guangdong Province (B2020011 and A2019473), the Natural Science Foundation of Guangdong Province (2016A030307033 and 2019A1515011107), Guangdong Provincial Joint Fund Youth Project (2020A1515110117), the Science and Technology Foundation of Qingyuan City (2020KJJH009), and the Medical Research Foundation of Qingyuan People’s Hospital (No. 20190206, and No. 20190205).
Ethics approval and consent to participate
The current study investigated the publicly available data, and no ethical approval was required.
Consent for publication
The authors declare no conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Clinical characteristics of TCGA, GSE13507, and IMvigor210 datasets included in this study. Table S2. List of Reactome glycolysis and cholesterol-biosynthesis gene sets. Table S3. List of co-expressed glycolysis and cholesterol biosynthesis genes identified by consensus clustering. Table S4. Multivariate Cox regression determining the association between metabolic subtypes and overall survival. Table S5. Frequency of SNV and indels mutation of the 20 co-expressed glycolytic and cholesterogenic genes in the TCGA dataset. Table S6. Frequency of CNV alterations of the 20 co-expressed glycolytic and cholesterogenic genes in the TCGA dataset. Table S7. Frequency of genetic alterations of the 30 most frequently altered genes across the 4 metabolic subtypes of bladder cancer in the TCGA dataset. Table S8. Correlations between expression of glycolytic and cholesterogenic genes and the 30 most frequently altered genes in TCGA dataset.
The Principal component analysis of samples of TCGA, GSE13507, and IMvigor210 before and after batch effect correction.
Kaplan-Meier curves showing the overall survival of patients stratified by high- and low levels of M0 or M2 macrophages infiltration.
About this article
Cite this article
Zhang, Y., Zhu, B., Cai, Y. et al. Alteration in glycolytic/cholesterogenic gene expression is associated with bladder cancer prognosis and immune cell infiltration. BMC Cancer 22, 2 (2022). https://doi.org/10.1186/s12885-021-09064-0
- Bladder cancer
- Cholesterol synthesis
- Metabolic subtypes
- Immune microenvironment