Skip to main content
Download PDF

Identifying prognostic characteristics of m6A-related glycolysis gene and predicting the immune infiltration landscape in bladder cancer

Cancer Cell International volume 23, Article number: 300 (2023) Cite this article

Abstract

Backgrounds

Glucose metabolism is associated with the development of cancers, and m6A RNA methylation regulator-related genes play vital roles in bladder urothelial carcinoma (BLCA). However, the role of m6A-related glucose metabolism genes in BLCA occurrence and development has not yet been reported. Our study aims to integrate m6A- and glycolysis-related genes and find potential gene targets for clinical diagnosis and prognosis of BLCA patients.

Methods

Sequencing data and clinical information on BLCA were extracted from common databases. Univariate Cox analysis was used to screen prognosis-related m6A glucose metabolism genes; BLCA subtypes were distinguished using consensus clustering analysis. Subsequently, genes associated with BLCA occurrence and development were identified using the “limma” R package. The risk score was then calculated, and a nomogram was constructed to predict survival rate of BLCA patients. Functional and immune microenvironment analyses were performed to explore potential functions and mechanisms of the different risk groups.

Results

Based on 70 prognosis-related m6A glucose metabolism genes, BLCA was classified into two subtypes, and 34 genes associated with its occurrence and development were identified. Enrichment analysis revealed an association of genes in high-risk groups with tricarboxylic acid cycle function and glycolysis. Moreover, significantly higher levels of seven immune checkpoints, 14 immune checkpoint inhibitors, and 32 immune factors were found in high-risk score groups.

Conclusions

This study identified two biomarkers associated with BLCA prognosis; these findings may deepen our understanding of the role of m6A-related glucose metabolism genes in BLCA development. We constructed a m6A-related glucose metabolism- and immune-related gene risk model, which could effectively predict patient prognosis and immunotherapy response and guide individualized immunotherapy.

Background

Bladder urothelial carcinoma (BLCA) is the most prevalent malignancy of the urinary system, with an estimated 573,000 new cases and 213,000 deaths in 2020, and is ranked as the 10th most common malignancy worldwide [1]. BLCA occurrence and development is a complex pathophysiological process, involving factors such as genetics and external environment. Currently recognized risk factors for BLCA include smoking and long-term exposure to chemicals [2, 3]. Although the development of treatment strategies for BLCA has rapidly advanced in recent years, the carcinoma is still characterized by a high recurrence rate and poor prognosis. Therefore, studying the molecular mechanisms underlying BLCA is particularly important to identifying new therapeutic targets.

N6-methyladenine (m6A) is the most common modification in eukaryotic cells at the mRNA and non-coding RNA levels [4]. m6A modification is a reversible process catalyzed by methyltransferases (writers), demethylases (erasers), and m6A-binding proteins (readers), and its aberrant regulation is associated with the development of multiple cancers [5, 6]. Reprogramming of energy metabolism is a hallmark of malignant tumors. The most prominent feature of this process is aerobic glycolysis, also known as the Warburg effect, where cancer cells rely on glycolysis for energy even when sufficient oxygen is present [7]. In addition to providing energy for tumor proliferation, this process generates metabolites that reshape the tumor microenvironment (TME) [8]. Therefore, targeting reprogrammed metabolic pathways is a potentially valuable therapeutic strategy. The m6A modification is involved in regulating expression of various glycolysis-related genes; however, the role of m6A-related glucose metabolism genes in BLCA occurrence and development has not been reported.

In this study, we integrated m6A- and glycolysis-related genes, and analyzed transcriptome and clinical data of BLCA from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Based on the expression of m6A-related glycolytic genes associated with BLCA prognosis, BLCA cases were divided into subtypes. We identified critical differential genes related to BLCA occurrence and development, and constructed a prognostic model. Following this, we conducted independent prognostic and clinical correlation analyses of the model, and functional and immune microenvironment analyses of high- and low-risk groups. This study provides potential gene targets for clinical diagnosis and prognosis of BLCA patients.

Methods

Data extraction

RNA sequencing data, survival information, and clinical information of BLCA were downloaded from TCGA (https://portal.gdc.cancer.gov) and GEO databases (https://www.ncbi.nlm.nih.gov/geo/). The TCGA dataset, containing 411 BLCA samples and 19 normal control (HC) samples, was used as the training dataset. Within the dataset, 406 BLCA samples had both survival and clinical information. GSE32894, containing 308 BLCA samples of which 224 had both survival and clinical information, was used as the validation dataset to verify the availability of the prognostic model. We obtained a total of 753 glucose metabolism-related genes from the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb/), and 24 m6A RNA methylation regulator-related genes from the database of functional variants involved in RNA modification (RMVar, http://rmvar.renlab.org/index.html) (Additional file 1: Table S1) [9,10,11].

Screening of prognosis-related m6A glucose metabolism genes

A Spearman correlation was performed between m6A methylation regulatory-related and glucose metabolism-related genes, to obtain m6A-related glycolytic genes with |Cor|> 0.5 and p < 0.05. A univariate Cox analysis was constructed to determine prognostic value of these genes and identify prognosis-related m6A glucose metabolism genes (p < 0.05).

Survival and clinical characteristics analysis of different subtypes patients

Consensus clustering analysis is commonly used to divide samples into subtypes. The “ConsensusClusterPlus” R package (version 1.58.0) was used to analyze BLCA based on prognosis-related m6A glucose metabolism gene expression [12]. The best clustering method was selected based on cumulative distribution function values, and BLCA was divided into subtypes verified by principal component analysis.

Using the “survival” R package, Kaplan–Meier (KM) survival analysis was performed to study differences in survival rate between subtypes of patients (version 3.2–13) [13]. Clinical information (age, gender, overall survival, pathologic T, N, M, tumor stage) was collected. The chi-square test was used to select clinical characteristics of different subtypes, which were visualized using the “ComplexHeatmap” R package (version 1.14.0) [14].

Screening of genes associated with BLCA occurrence and development

mRNA expression levels between BLCA and HC samples, as well as between patient subtypes, were compared using the “limma” R package (version 3.50.3) (adj. p-value < 0.05 and |log2fold change|> 0.5) [15]. Following this, genes associated with BLCA occurrence and development (target genes) were screened by using the online tool “Venny 2.1.0” to intersect the two groups of differentially expressed genes. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the “clusterprofiler” R package to investigate the functions of target genes (version 4.2.2) (adj.p-value < 0.05) [16].

Construction of the prognostic model of BLCA

Biomarkers of BLCA obtained by Cox regression analysis were screened to construct a survival risk model. The risk score was calculated using the following algorithm, and high- and low-risk groups were delineated according to the median risk value.

$${\mathrm{Riskscore}}_{\mathrm{sample}}={\sum }_{\mathrm{n}=1}^{\mathrm{n}}{(\mathrm{Coef}}_{\mathrm{i}}*{\mathrm{x}}_{\mathrm{i}})$$

A KM survival analysis was performed between the two groups, and the accuracy of the survival risk model was predicted using the risk and receiver operating characteristic (ROC) curves. The GSE32894 dataset was used to verify the applicability of our model.

Independent prognostic and clinical correlation analyses

The risk scores between different subtypes of clinical characteristics (age, gender, tumor stage, pathologic T, pathologic N, and pathologic M) were compared using the “wilcox.test” function. Then, univariate and multivariate Cox analyses were employed to identify the independent predictors for BLCA. The nomogram was constructed based on the independent predictors to predict the 1-, 3-, and 5-years survival rate. Its validity was then verified by drawing a calibration curve as well as a decision curve (DCA) at 1, 3, and 5 years of the prediction model.

Analysis of function and immune microenvironment

Gene set enrichment analysis (GSEA) was used to explore functions of the high- and low-risk groups, and a ridge map was drawn based on the top five enrichment pathways.

The proportions of 22 immune cell types and functions in the high- and low-risk groups were calculated, and “wilcox.test” was used to compare the differences. Following this, Spearman’s correlations were performed between biomarkers and differentially immune cells, biomarkers and differentially immune functions, and risk score and differentially immune cells. The inter-group differences between the eight immune checkpoints, 24 immune checkpoint inhibitors, and 38 immune factors were calculated and compared using “wilcox.test” [17, 18]. Furthermore, immune phenotype scores (IPS) in different risk groups were calculated, and a Spearman’s correlation was performed to determine the relationship between IPS and risk score.

Drug sensitivity analysis

The sensitivity of 138 commonly used drugs between high-risk and low-risk patient groups was analyzed and compared using the “pRRophetic” package (version 0.5) [19] and the “wilcox.test”, respectively. Drug sensitivity is measured by the half-maximal inhibitory concentration (IC50) value, with lower IC50 value indicating higher sensitivity of cells to the drug.

Patient samples

Bladder cancer and paracancerous tissues were obtained from BLCA patients admitted to the Qilu Hospital of Shandong University, between 2021 and 2022. All participants were informed of the study and provided consent before undergoing surgery. This study was approved by the Institutional Review Board of the Qilu Hospital of Shandong University (No. 2020046).

Cell culture

Bladder cancer cell lines 5637 and UM-UC-3 were purchased from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). All cell lines were tested for mycoplasma and resulted negative. 5637 cells were cultured in RPMI 1640 medium (Gibco, USA) supplemented with 10% fetal bovine serum (Gibco, USA). UM-UC-3 cells were cultured in Dulbecco’s modified medium (Gibco, USA) supplemented with 10% FBS (Gibco, USA). All cell lines were cultured in a 5% CO2 incubator at 37 °C.

siRNA transfection

Cells were plated in six-well dishes and transfected with siRNA-PLA2G2F, siRNA-IP6K2, or negative control using Lipofectamine 3000 (Invitrogen, USA). All siRNA sequences are listed in Additional file 2: Table S2.

RNA isolation and quantitative reverse transcription polymerase chain reaction (RT-qPCR)

Total RNA was extracted from tissues and cell lines using TRIzol reagent (Invitrogen, USA) and used to synthesize cDNA using Evo M-MLV RT Premix (Accurate Biology, China). RT-qPCR was performed using Premix Pro Taq HS qPCR Kit (Accurate Biology, China) on a LightCycler 96 (Roche, USA). β-actin was used as an internal control. All assays were performed in triplicate, and data were analyzed using the 2−ΔΔCT method. All PCR primers were purchased from Accurate Biology (Shanghai, China), and are listed in Additional file 3: Table S3.

Cell counting kit-8 (CCK-8) and colony formation assay

Cells were seeded in 96-well plates, with approximately 2000 cells per well. The Cell Counting kit-8 (CCK-8) (Bioss, China) was used to detect cell proliferation 0, 24, 48, 72, and 96 h after seeding. Absorbance was measured at 450 nm using a spectrophotometer (Tecan, Switzerland).

For colony formation assay, cells were seeded into 6-well plates with approximately 800 cells per well, and incubated at 5% CO2 and 37 °C. Colonies were fixed with 100% methanol and stained with crystal violet. Colonies containing more than 50 cells were considered survivors.

Results

Survival differences existed between patients in cluster 1 and 2

Firstly, the expression of m6A methylation regulatory-related and glucose metabolism-related genes was retrieved in the TCGA-BLCA dataset. A total of 326 m6A-related glycolytic genes were identified by Spearman correlation analysis with |Cor|> 0.5 and p < 0.05 as the threshold of significance in BLCA samples (Additional file 1: Table S1), of which 70 were prognosis-related m6A glucose metabolism genes, as identified by univariate Cox analysis. Among these, IGF2BP3, IGF2BP2, SLC16A1, DSC2, TGFBI, P4HA2, ANXA1, and CALU were negative prognostic factors, and CDK10 was the most positive prognostic factor (Additional file 4: Table S4). Based on the expression pattern of prognosis-related m6A glucose metabolism genes, BLCA were divided into two subtypes by consensus clustering analysis (Fig. 1A, B).

Fig. 1
figure 1

Survival differences existed between patients in cluster 1 and 2. A, B Results of consistent cluster analysis used to divide BLCA patients into two subtypes. C KM curves of the two subtypes (p = 0.00738). D Clinical features of the two subtypes. **P < 0.01, ***P < 0.001

Full size image

The KM curves indicated a significant difference in survival between the two subtypes of patients (p = 0.00738); cluster 2 patients had a better prognosis than that of cluster 1 patients (Fig. 1C). In addition, four clinical features, including gender, pathological T, M, and tumor stage, were significantly different between subtypes (Fig. 1D).

34 target genes were associated with metabolism-related pathways in BLCA

There were 7119 upregulated and 4,216 downregulated genes between 411 BLCA and 19 HC samples (Fig. 2A, B), and 22 upregulated and 66 downregulated genes between cluster 1 and cluster 2 samples (Fig. 2C, D). Subsequently, a total of 34 genes associated with BLCA occurrence and development (target genes) were screened by intersecting the differential expressed genes (DEG)1 between BLCA and HC samples and DEG2 between cluster 1and 2 (Fig. 2E, Additional file 5: Table S5).

Fig. 2.
figure 2

34 target genes were associated with metabolism-related pathways in BLCA. A, B DEGs between 411 BLCA and 19 HC samples. C, D DEGs between cluster 1 and cluster 2. E Screening of m6A-related glycolytic genes by intersecting DEG1 and DEG2. F GO analysis of m6A-related glycolytic genes. G KEGG analysis of m6A-related glycolytic genes

Full size image

GO functional analysis showed that the 34 target genes were related to several terms, including carbohydrate catabolism, glycolysis, ATP metabolism, and pyruvate metabolism (Fig. 2F, Additional file 6: Table S6). As for KEGG pathways, these genes were significantly associated with amino acid biosynthesis, carbon metabolism, glycolysis, and gluconeogenesis (Fig. 2G, Additional file 7: Table S7).

The survival risk model could be used to predict the prognosis of BLCA patients

Firstly, 34 intersection genes were submitted to univariate Cox regression analysis, and seven genes were obtained with p-value < 0.05 (Fig. 3A). Then, stepwise multivariate Cox regression analysis was used to identify two biomarkers (IP6K2 and PLA2G2F), both of which were positive factors (Hazard Ratio [HR] < 1) (Fig. 3B). The risk score was calculated for each patient, and the median risk score (0.961) was used as the cutoff to divide patients (n = 406, with complete survival data) into the high-risk group (n = 203) and low-risk group (n = 203). The results of KM and risk curves showed significant survival differences between the two groups (Fig. 3C, D). The AUC value of the ROC curve at 1, 3, and 5 years was all greater than 0.6, indicating that the survival risk model could be used as a prognostic model effectively (Fig. 3E). In addition, a significant negative correlation was found between disease risk and expression of IP6K2 and PLA2G2F.

Fig. 3
figure 3

The survival risk model could be used to predict the prognosis of BLCA patients. A, B Results of Cox analysis identifying IP6K2 and PLA2G2F as positive factors (Hazard Ratio < 1). C, D Significant survival differences observed between the high- and low-risk groups. E ROC curve revealing the potential use of the survival risk model as a prognostic model. F–H Results of KM, risk and ROC curves in the validation cohort, consistent with the training cohort

Full size image

GSE32894 was used as a validation dataset to verify applicability of the model. The results of the KM, risk, and ROC curves were consistent with those of the training dataset (Fig. 3F–H).

Risk score was an independent risk factor for BLCA patients

Univariate and multivariate Cox regression analyses were conducted to evaluate the clinical parameters and risk score to assess their prognostic value. The results revealed that gender, tumor stage, and pathological T stage were associated with risk score (Fig. 4A). Four clinical factors, including age, pathology T, N, and risk score, were significantly associated with patient survival (Fig. 4B, C). A nomogram based on these clinical factors was constructed to predict the patients’ 1-, 3-, and 5-year survival rates; the calibration curve showed that the slopes were closest to 1 (Fig. 4D, E). In addition, the nomogram model had a higher benefit rate than those of others, and the “number high risk” curves coincided with “number high risk with event” curves. These results indicated that the nomogram model has accurate predictive ability for BLCA (Fig. 4F, G).

Fig. 4
figure 4

Risk score was an independent risk factor for BLCA patients. A Gender, tumor stage, and pathologic T were associated with risk score. B, C Age, pathology T, N, and risk score were negatively associated with patient survival. D, E The nomogram for predicting survival in BLCA patients (1-, 3- and 5-year). F, G Evaluation of the accuracy of prediction using the DCA curve (F) and Clinical Impact Curve (G). T (Tumor) refers to the extent of the primary tumor; N (Node) refers to the involvement of regional lymph nodes; M (Metastasis) indicates the presence of metastasis

Full size image

Difference in function enrichment between the high- and low-risk groups

To investigate the potential molecular mechanisms of risk model, we performed GSEA enrichment analysis between the high- and low-risk groups. The GSEA results for the two groups are shown in Fig. 5A–D. T-cell activation, antigen receptor-mediated signaling pathway, B-cell differentiation, B-cell-mediated immunity, and B-cell receptor signaling pathway were enriched in the high-risk groups (Fig. 5A), whereas cellular glucuronidation, estrogen metabolism, miRNA-mediated gene silencing by inhibition of translation, sodium ion homeostasis, and translation repressor activity were enriched in the low-risk groups (Fig. 5B). The genes expressed in the high-risk groups were involved in antigen processing and presentation, cell adhesion molecules (CAMs), chemokine signaling, hematopoietic cell lineage, and JAK/STAT signaling (Fig. 5C). Meanwhile, genes expressed in the low-risk groups were involved in ascorbate and aldarate metabolism, drug metabolism, other enzymes, metabolism of xenobiotics by cytochrome p450, porphyrin and chlorophyll metabolism, and steroid hormone biosynthesis (Fig. 5D). The GSEA results revealed potential pathways or mechanisms that were activated during tumorigenesis and development, which can help us evaluate the prognosis of BLCA.

Fig. 5
figure 5

The difference in function enrichment between the high- and low-risk groups A GO analysis showing enrichment of T-cell activation, antigen receptor-mediated signaling pathway, B-cell differentiation, B-cell mediated immunity, and B-cell receptor signaling pathway in the high-risk groups. B GO analysis showing enrichment of cellular glucuronidation, estrogen metabolism, miRNA-mediated gene silencing by inhibition of translation, sodium ion homeostasis, and translation repressor activity in the low-risk groups. C KEGG analysis showing enrichment of antigen processing and presentation and of CAMs in the high-risk groups. D KEGG analysis showing enrichment of ascorbate and aldarate metabolism, drug metabolism, and of other enzymes, in the low-risk groups

Full size image

Significant differences were observed in immune microenvironment between the two risk groups

The CIBERSORT algorithm was used to investigate the differences in immune cell infiltration between the high- and low-risk groups. Seven immune cells, including naive B cells, plasma cells, activated memory CD4 + T cells, regulatory T cells, M1 macrophages, M2 macrophages, and activated dendritic cells, were found to have significantly different expression levels in the high- and low-risk groups. Among these, the immune infiltration of activated memory CD4 + T cells, M1 and M2 macrophages was significantly higher in the high-risk groups (Fig. 6A, B). IP6K2 expression correlated positively with M2 macrophages, while PLA2G2F expression correlated positively with M2 macrophages and negatively with activated dendritic cells (Fig. 6C). Meanwhile, risk score correlated positively with activated memory CD4 + T cells, M1 and M2 macrophages, and negatively with regulatory T cells and activated dendritic cells (Fig. 6D). In addition, 28 immune functions were significantly different between the high- and low-risk groups; among these, IP6K2 expression correlated positively with MHC class 1 and NK cell function (Fig. 6E, F). In summary, we speculated that immune microenvironment was associated with the occurrence and development of BLCA.

Fig. 6
figure 6

Significant differences were observed in immune microenvironment between the two risk groups. A, B Seven differential immune cells were identified in the high- and low-risk groups. C Heat map of correlation between biomarkers and differential immune cells. D Correlation between risk scores and differential immune cells. E Box plot of differential immune functions between the high- and low-risk groups. F Heat map of correlation between biomarkers and differential immune function. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Full size image

The risk score can predict the effectiveness of immunotherapy

As immunotherapy became a promising strategy for the treatment of tumors, we evaluated differences in immunotherapy response between high- and low-risk groups. The expression levels of seven immune checkpoints (CD274, CTLA-4, LAG-3, HAVCR2, PDCD1, PDCD1LG2 and TIGIT), 14 immune checkpoint inhibitors (LAG3, CD40LG, CD86, CD48, CD274, CD244, CD70, CD27, CTLA4, ICOS, CD28, TIGIT, BTLA, TNFRSF4) and 32 immune factors were significantly higher in the high-risk groups, while those of two immune checkpoint inhibitors (BTNL2 and TNFRSF25) and three immune factors (CCR7, CCR10, CCL18) were significantly higher in the low-risk groups (Fig. 7A–C). Furthermore, the risk score significantly negatively correlated with the IPS score (Fig. 7D), indicating that patients in the low-risk group were more likely to benefit from immunotherapy.

Fig. 7
figure 7

The risk score can predict the effectiveness of immunotherapy. A–C Expression of seven immune checkpoints, 14 ICIs and 32 immune factors in the high- and low-risk groups. D IPS scores in the high- and low-risk groups. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

Full size image

The risk score model can predict the drug sensitivity between high- and low-risk groups

Considering that the risk score model could be used to predict the prognosis of BLCA patients, we further investigated the relationship between risk score and chemotherapy resistance. The results revealed that there were significant differences in the IC50 values of 114 drugs between the high- and low-risk groups. (Additional file 8: Table S8), and we further analyze the differences in drug sensitivity between the high- and low-risk groups for 10 drugs, such as parthenolide, sunitinib, cyclopamine, bexarotene, etc. (Additional file 9: Fig. S9). Therefore, we believe that this risk score model can be utilized to predict the sensitivity of BLCA patients to chemotherapy drugs, which could improve the effectiveness of chemotherapy.

IP6K2 and PLA2G2F were downregulated in bladder cancer tissues and could regulate bladder cancer cell proliferation in vitro

IP6K2 and PLA2G2F expression in bladder cancer and adjacent tissues was detected using RT-qPCR. Expression levels of IP6K2 and PLA2G2F were significantly downregulated in cancer tissues compared with their expression levels in adjacent tissues (Fig. 8A).

Fig. 8
figure 8

IP6K2 and PLA2G2F were downregulated in bladder cancer tissues and could regulate bladder cancer cell proliferation in vitro. A IP6K2 and PLA2G2F expression in bladder cancer and matched adjacent normal tissues detected by RT-qPCR. B Efficiency of IP6K2 and PLA2G2F knockdown verified by RT-qPCR. C The proliferative ability of IP6K2 and PLA2G2F knockdown cells was determined using CCK-8 assay. D The colony-forming ability of IP6K2 or PLA2G2F knockdown cells was determined using colony formation assay. *P < 0.05, **P < 0.01, ***P < 0.001

Full size image

Effect of IP6K2 and PLA2G2F on cell proliferation in vitro was explored to determine their role in bladder cancer cell proliferation. Efficiency of IP6K2 and PLA2G2F knockdown was verified using RT-qPCR (Fig. 8B). CCK-8 assays demonstrated that IP6K2 knockdown significantly promoted cancer cell proliferation in 5637 and UM-UC-3 cells (Fig. 8C). Moreover, the colony-forming ability of 5637 and UM-UC-3 cells enhanced after IP6K2 or PLA2G2F knockdown (Fig. 8D).

Discussion

As bladder cancer has a high risk of recurrence after treatment, it is important to explore its molecular mechanisms and discover novel therapeutic targets. Aberrant regulation of m6A modification, a hotspot in cancer research, is closely associated with tumorigenesis and metastasis [5]. Dysregulation of m6A modification may impact the expression of genes related to tumor progression, some of which are associated with glycolysis [6, 20, 21]. The relationship between m6A modification and glycolysis in tumors is an area of active research and has important implications in understanding the molecular mechanisms underlying tumorigenesis and cancer progression. m6A modification can influence the expression of glycolysis-related genes by regulating their mRNAs’ stability and translation [6, 20]. This can activate the glycolytic pathway in tumor cells, thereby providing more energy to meet their rapid growth and proliferation demands. Moreover, tumor cells often rely on glycolysis for their energy needs, and this can create a local environment with increased lactate production and decreased pH value. This acidic microenvironment can influence immune cell function and suppress immune responses, making it easier for tumor cells to evade the surveillance and elimination of immune system [22, 23]. In this study, we identified a group of glucose metabolism genes associated with m6A regulator-related genes, aiming to further explore their roles in the occurrence and development of bladder cancer.

We identified IP6K2 and PLA2G2F as BLCA biomarkers using Cox regression analysis. Two risk groups with different prognoses were demonstrated, which may help develop different therapeutic strategies. The survival risk model was verified using a ROC curve, which could be used to predict the prognosis of BLCA patients. Moreover, the analysis results of validation datasets were consistent with those of the training dataset, suggesting a reliable and effective survival risk model.

IP6K2 converts InsP6 into InsP7/PP-InsP5 [24], all of which are important members of the IPK family. Additionally, it mediates apoptosis by facilitating the DNA damage-induced cell death pathway [25, 26]. We found that IP6K2 knockdown promotes bladder cancer cell proliferation and colony formation. Aberrant PLA2G2F expression is associated with several malignancies [27, 28]. In addition, PLA2G2F knockdown also promotes proliferation and colony formation of bladder cancer cells.

Aerobic glycolysis provides energy for tumor proliferation; in addition, it generates various tumor metabolites that reshape the tumor microenvironment (TME). Specifically, tumor aerobic glycolysis can lead to nutrient deprivation in the microenvironment and induce hypoxia with extracellular lactic acid accumulation, thus forming an inhibitory immune microenvironment [8]. The TME consists of dynamically balanced tumor cells, immune-infiltrating cells, fibroblasts, and cytokines. Tumor-infiltrating immune cells (TIICs) play important roles in tumor growth, invasion, and metastasis, making them effective targets for immunotherapy [29, 30]. TIICs are closely related to clinical outcomes of BLCA and can serve as effective targets for the development of new immunotherapies for BLCA patients [31,32,33]. To explore immune cell infiltration in the high- and low-risk groups, we conducted immune cell infiltration analysis, using the CIBERSORT algorithm to obtain the infiltration of each type of immune cell. Significant differences in the infiltration of 28 immune functions were found between the high- and low-risk groups. Among them, IP6K2 was significantly positively correlated with MHC class 1 and NK cell function.

After analyzing the expression of immune checkpoints in the high- and low-risk groups, we found significant upregulation of 7 immune checkpoints, 14 ICIs, and 32 immune factors in the high-risk groups, and of two ICIs and three immune factors in the low-risk groups. Moreover, IPS was significantly higher in the low-risk groups, indicating greater sensitivity of these patients to immunotherapy. Immune checkpoint receptors are located on the surface of immune cells. After binding to ligands, they play a role in regulate immune response and maintain self-tolerance [34]. Cancer cells use this feature to evade immune system surveillance, leading to immune escape [35, 36]. Immune checkpoint inhibitors (ICIs) are a novel class of cancer therapeutic drugs that enhance the anti-tumor response of the immune system by inhibiting immune escape of tumor cells. By altering the m6A modification of specific genes, the expression of tumor-associated antigens can be up-regulated, thereby enhancing the anti-tumor immune response [37, 38]. Most tumor cells exhibit enhanced glycolytic capacity, which enables them to survive in hypoxic and nutrient-poor microenvironment and provides energy for rapid growth. This process leads to the production of large amounts of lactate and other metabolic byproducts, which can inhibit the function of immune cells [39]. One of the goals of immunotherapy is to reverse the suppressive effect of tumor cells on immune cells by interfering with their glucose metabolism. ICIs eliminate tumor cells by restoring the activity of T cells. These drugs can alter the metabolism of T cells, increasing their demand for glucose in order to enhance their tumor killing activity [39]. Simultaneously, it is necessary to reduce the inhibitory effect on the immune system by intervening in the process of tumor cells’ glucose metabolism [40, 41]. Therefore, it is crucial to fully understand and utilize these relationships in order to optimize immunotherapy strategies and improve the effectiveness of immunotherapy.

Conclusions

In conclusion, our study identified two biomarkers (IP6K2 and PLA2G2F) associated with BLCA prognosis. We constructed a risk model which could effectively predict patient prognosis and immunotherapy response and guide individualized immunotherapy. However, further validation using basic experiments is required to confirm our findings.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Abbreviations

BLCA:

Bladder urothelial carcinoma

m6A:

N6-methyladenine

TME:

Tumor microenvironment

TCGA:

The Cancer Genome Atlas

GEO:

Gene Expression Omnibus

K-M curves:

Kaplan–Meier (K–M) curves

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

ROC:

Receiver operating characteristic

DCA:

Decision curve analysis

GSEA:

Gene set enrichment analysis

IPS:

Immune phenotype scores

IC50:

Half-maximal inhibitory concentration

CCK-8:

Cell counting kit-8

DEGs:

Differentially expressed genes

TIICs:

Tumor-infiltrating immune cells

ICIs:

Immune checkpoint inhibitors

References

  1. 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.

    Article  PubMed  Google Scholar 

  2. Cumberbatch MGK, Jubber I, Black PC, Esperto F, Figueroa JD, Kamat AM, et al. Epidemiology of bladder cancer: a systematic review and contemporary update of risk factors in 2018. Eur Urol. 2018;74(6):784–95.

    Article  PubMed  Google Scholar 

  3. van Osch FH, Jochems SH, van Schooten FJ, Bryan RT, Zeegers MP. Quantified relations between exposure to tobacco smoking and bladder cancer risk: a meta-analysis of 89 observational studies. Int J Epidemiol. 2016;45(3):857–70.

    Article  PubMed  Google Scholar 

  4. Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m(6)A RNA methylation. Nat Rev Genet. 2014;15(5):293–306.

    Article  CAS  PubMed  Google Scholar 

  5. Sibbritt T, Patel HR, Preiss T. Mapping and significance of the mRNA methylome. Wiley Interdiscip Rev RNA. 2013;4(4):397–422.

    Article  CAS  PubMed  Google Scholar 

  6. Liu Y, Liang G, Xu H, Dong W, Dong Z, Qiu Z, et al. Tumors exploit FTO-mediated regulation of glycolytic metabolism to evade immune surveillance. Cell Metab. 2021;33(6):1221-33 e11.

    Article  CAS  PubMed  Google Scholar 

  7. Garber K. Energy boost: the Warburg effect returns in a new theory of cancer. J Natl Cancer Inst. 2004;96(24):1805–6.

    Article  PubMed  Google Scholar 

  8. Pavlova NN, Zhu J, Thompson CB. The hallmarks of cancer metabolism: still emerging. Cell Metab. 2022;34(3):355–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Zhang C, Wang M, Ji F, Peng Y, Wang B, Zhao J, et al. A novel glucose metabolism-related gene signature for overall survival prediction in patients with glioblastoma. Biomed Res Int. 2021;2021:8872977.

    PubMed  PubMed Central  Google Scholar 

  10. Luo X, Li H, Liang J, Zhao Q, Xie Y, Ren J, et al. RMVar: an updated database of functional variants involved in RNA modifications. Nucleic Acids Res. 2021;49(D1):D1405–12.

    Article  CAS  PubMed  Google Scholar 

  11. Zhong J, Liu Z, Cai C, Duan X, Deng T, Zeng G. m(6)A modification patterns and tumor immune landscape in clear cell renal carcinoma. J Immunother Cancer. 2021;9(2): e001646.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Jung RE, Pjetursson BE, Glauser R, Zembic A, Zwahlen M, Lang NP. A systematic review of the 5-year survival and complication rates of implant-supported single crowns. Clin Oral Implants Res. 2008;19(2):119–30.

    Article  PubMed  Google Scholar 

  14. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9.

    Article  CAS  PubMed  Google Scholar 

  15. Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016;44(8): e71.

    Article  PubMed  Google Scholar 

  16. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3): 100141.

    CAS  PubMed  PubMed Central  Google Scholar 

  17. Lin J, Zhao A, Fu D. Evaluating the tumor immune profile based on a three-gene prognostic risk model in HER2 positive breast cancer. Sci Rep. 2022;12(1):9311.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhang X, Zhang X, Li G, Hao Y, Liu L, Zhang L, et al. A novel necroptosis-associated lncRNA signature can impact the immune status and predict the outcome of breast cancer. J Immunol Res. 2022;2022:3143511.

    PubMed  PubMed Central  Google Scholar 

  19. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9(9): e107468.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Huang J, Sun W, Wang Z, Lv C, Zhang T, Zhang D, et al. FTO suppresses glycolysis and growth of papillary thyroid cancer via decreasing stability of APOE mRNA in an N6-methyladenosine-dependent manner. J Exp Clin Cancer Res. 2022;41(1):42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Xu Y, Song M, Hong Z, Chen W, Zhang Q, Zhou J, et al. The N6-methyladenosine METTL3 regulates tumorigenesis and glycolysis by mediating m6A methylation of the tumor suppressor LATS1 in breast cancer. J Exp Clin Cancer Res. 2023;42(1):10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Guo D, Tong Y, Jiang X, Meng Y, Jiang H, Du L, et al. Aerobic glycolysis promotes tumor immune evasion by hexokinase2-mediated phosphorylation of IkappaBalpha. Cell Metab. 2022;34(9):1312–24.

    Article  CAS  PubMed  Google Scholar 

  23. Hirschhaeuser F, Sattler UG, Mueller-Klieser W. Lactate: a metabolic key player in cancer. Cancer Res. 2011;71(22):6921–5.

    Article  CAS  PubMed  Google Scholar 

  24. Sarmah B, Wente SR. Inositol hexakisphosphate kinase-2 acts as an effector of the vertebrate Hedgehog pathway. Proc Natl Acad Sci USA. 2010;107(46):19921–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Morrison BH, Bauer JA, Kalvakolanu DV, Lindner DJ. Inositol hexakisphosphate kinase 2 mediates growth suppressive and apoptotic effects of interferon-beta in ovarian carcinoma cells. J Biol Chem. 2001;276(27):24965–70.

    Article  CAS  PubMed  Google Scholar 

  26. Chakraborty A, Koldobskiy MA, Sixt KM, Juluri KR, Mustafa AK, Snowman AM, et al. HSP90 regulates cell survival via inositol hexakisphosphate kinase-2. Proc Natl Acad Sci USA. 2008;105(4):1134–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yamamoto K, Miki Y, Sato M, Taketomi Y, Nishito Y, Taya C, et al. The role of group IIF-secreted phospholipase A2 in epidermal homeostasis and hyperplasia. J Exp Med. 2015;212(11):1901–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Mounier CM, Wendum D, Greenspan E, Flejou JF, Rosenberg DW, Lambeau G. Distinct expression pattern of the full set of secreted phospholipases A2 in human colorectal adenocarcinomas: sPLA2-III as a biomarker candidate. Br J Cancer. 2008;98(3):587–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Zhang Y, Zhang Z. The history and advances in cancer immunotherapy: understanding the characteristics of tumor-infiltrating immune cells and their therapeutic implications. Cell Mol Immunol. 2020;17(8):807–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Pan Y, Lu F, Fei Q, Yu X, Xiong P, Yu X, et al. Single-cell RNA sequencing reveals compartmental remodeling of tumor-infiltrating immune cells induced by anti-CD47 targeting in pancreatic cancer. J Hematol Oncol. 2019;12(1):124.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. van Wilpe S, Gerretsen ECF, van der Heijden AG, de Vries IJM, Gerritsen WR, Mehra N. Prognostic and predictive value of tumor-infiltrating immune cells in urothelial cancer of the bladder. Cancers (Basel). 2020;12(9):2692.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Eckstein M, Lieb V, Jung R, Sikic D, Weigelt K, Stohr R, et al. Combination of GP88 expression in tumor cells and tumor-infiltrating immune cells is an independent prognostic factor for bladder cancer patients. Cells. 2021;10(7):1796.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Li F, Guo H, Wang Y, Liu B, Zhou H. Profiles of tumor-infiltrating immune cells and prognostic genes associated with the microenvironment of bladder cancer. Int Immunopharmacol. 2020;85: 106641.

    Article  CAS  PubMed  Google Scholar 

  34. Borcherding N, Kolb R, Gullicksrud J, Vikas P, Zhu Y, Zhang W. Keeping tumors in check: a mechanistic review of clinical response and resistance to immune checkpoint blockade in cancer. J Mol Biol. 2018;430(14):2014–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Topalian SL, Taube JM, Anders RA, Pardoll DM. Mechanism-driven biomarkers to guide immune checkpoint blockade in cancer therapy. Nat Rev Cancer. 2016;16(5):275–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. He W, Xiao K, Xu J, Guan W, Xie S, Wang K, et al. Recurrent sepsis exacerbates CD4(+) T cell exhaustion and decreases antiviral immune responses. Front Immunol. 2021;12: 627435.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Zhang L, Li Y, Zhou L, Zhou H, Ye L, Ou T, et al. The m6A reader YTHDF2 promotes bladder cancer progression by suppressing RIG-I-mediated immune response. Cancer Res. 2023;83(11):1834–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zong X, Xiao X, Shen B, Jiang Q, Wang H, Lu Z, et al. The N6-methyladenosine RNA-binding protein YTHDF1 modulates the translation of TRAF6 to mediate the intestinal immune response. Nucleic Acids Res. 2021;49(10):5537–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Cascone T, McKenzie JA, Mbofung RM, Punt S, Wang Z, Xu C, et al. Increased tumor glycolysis characterizes immune resistance to adoptive T cell therapy. Cell Metab. 2018;27(5):977–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Wang ZH, Peng WB, Zhang P, Yang XP, Zhou Q. Lactate in the tumour microenvironment: from immune modulation to therapy. EBioMedicine. 2021;73: 103627.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Feng Q, Liu Z, Yu X, Huang T, Chen J, Wang J, et al. Lactate increases stemness of CD8 + T cells to augment anti-tumor immunity. Nat Commun. 2022;13(1):4981.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We sincerely thank the editor and reviewers whose comments and suggestions helped improve this article.

Funding

This research was supported by the Natural Science Foundation of Shandong Province (ZR2020MH082).

Author information

Author notes

  1. Lijian Gao and Xianzhou Jiang are authors contributed equally to this work.

Authors and Affiliations

  1. Department of Urology, Qilu Hospital of Shandong University, Jinan, 250012, China

    Guanwen Zhou, Yi Li, Xiangguo Ren, Guoliang Qin, Zhaocun Zhang, Haifeng Zhao, Lijian Gao & Xianzhou Jiang

  2. Department of Urology, Qilu Hospital of Shandong University Dezhou Hospital, Dezhou, 253000, China

    Lijian Gao

Authors
  1. Guanwen Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yi Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xiangguo Ren
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Guoliang Qin
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zhaocun Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Haifeng Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Lijian Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Xianzhou Jiang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JXZ and ZGW designed, edited, and led out this study’s experiments. QGL and ZHF provided study material. ZGW performed the in vitro experiments. LY and RXG carried out the statistical analyses. ZGW and GLJ wrote and completed the paper. JXZ and ZGW have access to all the data in this study and take responsibility for the integrity and accuracy of the data analyses. All authors reviewed and approved the manuscript.

Corresponding authors

Correspondence to Lijian Gao or Xianzhou Jiang.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Institutional Review Board of Qilu Hospital (No. 2020046). Informed consent was obtained from all individual participants included in the study.

Consent for publication

Additional informed consent was obtained from all individual participants for whom identifying information is included in this article.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Glycolysis-related genes and m6a-related genes.

Additional file 2: Table S2.

siRNA sequences used in this study.

Additional file 3: Table S3.

Primes sequences for RT-qPCR used in this study.

Additional file 4: Table S4.

Univariate Cox analysis of m6A-related glycolytic genes

Additional file 5: Table S5.

m6A-related glycolytic genes

Additional file 6: Table S6.

GO analysis of m6A-related glycolytic genes.

Additional file 7: Table S7.

KEGG analysis of m6A-related glycolytic genes

Additional file 8: Table S8.

Sensitivity analyses of 138 drugs in patients between high- and low-risk groups

Additional file 9: Figure S1.

Differences in drug sensitivity between high- and low-risk groups.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, G., Li, Y., Ren, X. et al. Identifying prognostic characteristics of m6A-related glycolysis gene and predicting the immune infiltration landscape in bladder cancer. Cancer Cell Int 23, 300 (2023). https://doi.org/10.1186/s12935-023-03160-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12935-023-03160-w

Keywords

Cancer Cell International

ISSN: 1475-2867

Contact us