Glycolysis-related five-gene signature correlates with prognosis and immune infiltration in gastric cancer

Abstract

BACKGROUND

Gastric cancer (GC) is one of the most common malignancies worldwide. Glycolysis has been demonstrated to be pivotal for the carcinogenesis of GC.

AIM

To develop a glycolysis-based gene signature for prognostic evaluation in GC patients.

METHODS

Differentially expressed genes correlated with glycolysis were identified in stomach adenocarcinoma data (STAD). A risk score was established through a univariate Cox and least absolute shrinkage and selection operator analysis. The model was evaluated using the area under the receiver operating characteristic curves. RNA-sequencing data from high- and low-glycolysis groups of STAD patients were analyzed using Cibersort algorithm and Spearman correlation to analyze the interaction of immune cell infiltration and glycolysis. Multiomics characteristics in different glycolysis status were also analyzed.

RESULTS

A five-gene signature comprising syndecan 2, versican, malic enzyme 1, pyruvate carboxylase and SRY-box transcription factor 9 was constructed. Patients were separated to high- or low-glycolysis groups according to risk scores. Overall survival of patients with high glycolysis was poorer. The sensitivity and specificity of the model in prediction of survival of GC patients were also observed by receiver operating characteristic curves. A nomogram including clinicopathological characteristics and the risk score also showed good prediction for 3- and 5-year overall survival. Gene set variation analysis showed that high-glycolysis patients were related to dysregulation of pancreas beta cells and estrogen late pathways, and low-glycolysis patients were related to Myc targets, oxidative phosphorylation, mechanistic target of rapamycin complex 1 signaling and G2M checkpoint pathways. Tumor-infiltrating immune cells and multiomics analysis suggested that the different glycolysis status was significantly correlated with multiple immune cell infiltration. The patients with high glycolysis had lower tumor mutational burden and neoantigen load, higher incidence of microsatellite instability and lower chemosensitivity. High glycolysis status was often found among patients with grade 2/3 cancer or poor prognosis.

CONCLUSION

The genetic characteristics revealed by glycolysis could predict the prognosis of GC. High glycolysis significantly affects GC phenotype, but the detailed mechanism needs to be further studied.

Key Words: Glycolysis; Tumor microenvironment; Immune infiltration; Prognosis; Gastric cancer

Core Tip: A five-gene signature consisting of syndecan 2, versican, malic enzyme 1, pyruvate carboxylase, and SRY-box transcription factor 9 was constructed. Based on their risk scores, the patients were stratified into high- or low-glycolysis groups. The patients with high glycolysis were significantly poorer. The high-glycolysis patients revealed dysregulation of pancreas beta cells and estrogen late pathways, while low-glycolysis patients showed associations with Myc targets, oxidative phosphorylation, mechanistic target of rapamycin complex 1 signaling, and G2M checkpoint pathways. A significant correlation between different glycolysis statuses and multiple immune cell infiltrations was observed. The high glycolysis status was associated with lower tumor mutational burden, higher incidence of microsatellite instability-high, as well as reduced chemosensitivity among gastric cancer patients.

INTRODUCTION

Gastric cancer (GC) is a common cancer, with about 1 million new cases in 2020 and approximately 769000 deaths, and has the fifth highest incidence and fourth highest mortality[1]. Currently, surgery and chemotherapy are widely used to treat GC. The 5-year survival of GC remains low and many patients are diagnosed late due to its strong heterogeneity and complex tumor microenvironment (TME)[2]. Hence, how to diagnose GC early, explore the molecular mechanisms and improve survival rate have become the focus of our research.

TME includes various cells and extracellular matrix and plays different roles in GC growth and progression[3]. The different cellular components including endothelial cells, immune cells and fibroblasts in the TME play important roles in cancer development, immune responses and cell adhesion[4,5]. It has been suggested that TME could regulate immunity[6], signal transduction[7], drug sensitivity[8] and cell metabolism[9]. Similar to other types of tumor cells, GC cells prefer aerobic glycolysis, even with sufficient oxygen[10,11]. This special form of glucose metabolism in GC, also called the Warburg effect, increases lactate production, resulting in a unique pH in the TME[12]. The acidic TME caused by this glycolytic metabolic process eventually leads to immune escape of tumor cells, more invasive and metastatic characteristics, and induces resistance to radiotherapy and chemotherapy[13-16]. In contrast, immune or inflammatory factors produced in the TME can also regulate glycolysis to affect the fate of tumor cells[17,18]. This shows that glycolysis plays an important role for tumor cells, maintaining oncogenicity in the TME. Therefore, in-depth exploration of the relationship between glycolysis and TME in GC will benefit diagnosis and treatment of GC.

As high-throughput sequencing has developed, some clusters of molecular biomarkers screened out in specific environments can provide valuable potential treatment strategies for GC patients with poor survival or high mortality[19]. In contrast, the traditional tumor-node-metastasis (TNM) staging for GC has shown limited predictive ability[19]. Two previous studies have reported glycolysis-related gene signatures, although some limitations, including insufficient sample size for validation or multiomics analysis were suggested[20,21]. In this study, we incorporated different glycolysis status data to construct a novel glycolysis-related gene signature based on comprehensive analysis, which may provide guidance for follow-up research on glycolysis and clinical prediction.

MATERIALS AND METHODS

Data collection

The study workflow is shown in Figure 1. RNA sequencing FPKM data and relevant clinical records, including 375 GC samples and 32 controls were downloaded from The Cancer Genome Atlas (TCGA) program (https://portal.gdc.cancer.gov/) (TCGA-STAD). Differentially expressed genes with |logFC| > 1 and P < 0.05 were screened using limma package. The mRNA expression data fles of GSE62254, GSE15459, GSE84437 and GSE26901 were obtained from https://www.ncbi.nlm.nih.gov/geo/ and used as external validation set. We complied with the TCGA publication guidelines and data access policies.

Figure 1 The workflow of construction and validation of the five-gene signature for gastric cancer. GSEA: Gene set enrichment analysis; GC: Gastric cancer; TCGA: The Cancer Genome Atlas; DEA: Differential expression analysis; LASSO: Least absolute shrinkage and selection operator; TMB: Tumor mutational burden; MSI: Microsatellite instability; K-M: Kaplan-Meier; ROC: Receiver operating characteristic.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes functional enrichment analysis

Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed utilizing the clusterProfifiler package in R. The protein-protein interaction network based on differentially expressed genes was shown by Cytoscape software.

Construction and evaluation of glycolysis-related gene signature

Prognosis-affecting glycolysis-related genes were filtered in the training cohort through univariate Cox analysis. Least absolute shrinkage and selection operator (LASSO) analysis was applied to further narrow the gene ranges based on previous studies[22]. The risk scores were calculated as follows based on previous studies[23]: Risk score = ∑(coefficient i × expression of signature gene i). Patients were separated into high- or low-glycolysis groups based on risk scores. Kaplan-Meier analysis was applied to compare survival rates. The performance of this signature was evaluated by the area under the receiver operating characteristic curve. The risk score was also tested as a risk factor using univariate or multivariate Cox regression.

Evaluation of the sensitivity of chemotherapeutic agents

According to the Genomics of Drug Sensitivity in Cancer (www.cancerrxgene.org/) cell line expression profiles and TCGA gene expression profiles, the pRRophetic algorithm was used to predict the half-maximal inhibitory concentration (IC50) of chemotherapy drugs and their sensitivity in different patients based on previous studies[24].

Correlation analysis between immune cell infiltration and glycolysis

Cibersort is widely used to describe the immune cell composition[25]. RNA sequencing data from high- and low-glycolysis groups of STAD patients were analyzed using Cibersort algorithm and Spearman correlation to determine how immune cell infiltration affected glycolysis.

Tumor mutational burden, neoantigen and microsatellite instability

Tumor mutational burden (TMB) measures the total number of mutations per megabyte of tumor tissue. These mutations include gene coding errors, base substitution insertions, or deletions. We also used NetMHCpan v3.0 for neoantigen prediction, as described previously[26]. Microsatellite instability (MSI) data were downloaded from TCGA. The differences in TMB, neoantigen and MSI between high- and low-glycolysis groups were investigated.

Gene set variation analysis

Gene set variation analysis (GSVA) was conducted using GSVA R packages to explore the differences (t score) in activation status of biological pathways between high- and low-glycolysis groups from TCGA STAD. GSVA is widely used for testing the variation in pathway activity in gene expression datasets[27]. mRNA data including glycolysis metabolism were from MSigDB (http://software.broadinstitute.org/gsea/index.jsp).

Nomogram

A prognostic nomogram including clinical features and risk scores was constructed for TCGA cohorts. The prognostic values of age, gender, tumor grade, TNM classification, and risk score were studied using univariate Cox analysis. Overall survival (OS)-related clinical features and risk score were introduced to multivariate Cox analysis to construct a nomogram, which was tested by C-index and calibration plots.

Data analysis

Data analysis was performed using R 3.6.0. The prognostic outcome was assessed by Kaplan-Meier analysis with log-rank test. Univariate or multivariate Cox regression was used to evaluate hazard ratios of prognostic factors. The differences of independent samples were assessed using the Wilcoxon test. P < 0.05 was defined as statistically significant.

RESULTS

Expression and enrichment analysis of glycolysis-related genes in the TCGA-STAD database

Two hundred glycolysis-related genes were downloaded, and their RNA expression data were retrieved. Only 199 glycolysis-related genes with complete expression data were screened out. The Wilcoxon test results showed that 111 differentially expressed glycolysis-related genes were obtained (49 upregulated and 62 downregulated) and they were used as candidate gene sets for subsequent modeling analysis (Figure 2A and B).

Figure 2 Differentially expressed glycolysis genes and their functional analysis. A: Heatmap of the upregulated (yellow) and downregulated (green) genes between control (blue) and tumor (red) tissues of gastric cancer (GC); B: Volcano plot of the upregulated (red) and downregulated (green) genes between the normal (blue) and tumor (red) samples of GC; C and D: Gene Ontology enrichment (C) and Kyoto Encyclopedia of Genes and Genomes (D) analysis of the 111 differential expressed glycolysis genes in GC; E: Interaction network of 111 differential expressed glycolysis genes. FC: Fold change.

To study further the potential function and molecular mechanism of these 199 genes, we applied GO and KEGG analyses. Nineteen enriched pathways and 303 GO terms are listed in Supplementary Table 1. The GO terms revealed that biological processes were mainly enriched in monosaccharide/hexose/glucose metabolic process, and carbohydrate biosynthetic process. Cellular components were mostly enriched in secretory granule and cytoplasmic vesicle lumen. Molecular functions were mostly enriched in carbohydrate/monosaccharide and coenzyme binding. KEGG pathway enrichment analysis of the first five terms demonstrated that those genes affected carbon metabolism, glycolysis/gluconeogenesis, and biosynthesis of amino acids, amino and nucleotide sugar metabolism, and hypoxia-inducible factor-1 signaling pathway (Figure 2C and D). The inter-relation of these genes is visualized in Figure 2E.

Identification and internal validation of survival-affecting glycolysis-related genes

Univariate Cox regression identified that syndecan 2 (SDC2), versican (VCAN), malic enzyme 1 (ME1), pyruvate carboxylase (PC) and SRY-box transcription factor 9 (SOX9) were associated with prognosis (Table 1). LASSO regression was used to choose optimal genes (Figure 3A and B, Table 2). These five genes were developed to build a prognostic model to analyze the glycolysis status of patients: Risk score = (-0.1885 × SOX9 expression) + (0.1062 × PC expression) + (0.1429 × SDC2 expression) + (0.1475 × VCAN expression) + (0.1674 × ME1 expression). We divided the TCGA-STAD randomly into training and validating sets (4:1). According to the scores, patients were separated into high- or low-glycolysis groups. Kaplan-Meier curves and log-rank tests indicated that high-glycolysis patients had poorer OS (Figure 3C and D). Area under the receiver operating characteristic curve at 1-, 3- and 5-year were 0.61, 0.56, and 0.74 in training sets, and 0.79, 0.85, and 0.85 in validation sets (Figure 3E and F), suggesting good sensitivity and specificity of this model in predicting survival prediction for GC patients. Cox regression showed that risk score had the highest hazard ratio, suggesting that risk score might be a definitely independent risk factor for the prognosis of STAD (Figure 3G and H).

Figure 3 Least absolute shrinkage and selection operator Cox regression analysis for the screening of prognostic genes and construction and validation of the risk score in training and test data. A: Choosing optimal lambda in least absolute shrinkage and selection operator; dotted lines were drawn at the optimal values; B: Coefficients of overall survival-related glycolytic genes; C and D: Kaplan-Meier analysis of different glycolysis status in training and internal testing data (yellow: High glycolysis status; blue: Low glycolysis status); E and F: Receiver operating characteristic curves of risk models. Red: 1 year, yellow: 3 years, and blue: 5 years; G and H: Univariate and multivariate Cox regression analysis of clinicopathological parameters including risk score to assess the prognostic value in The Cancer Genome Atlas of Stomach Adenocarcinoma. TCGA: The Cancer Genome Atlas; TPR: True positive rate; FPR: False positive rate; AUC: Area under the receiver operating characteristic curve.

Table 1 Cox proportional risk model of the five genes from the differentially expressed gene analysis based on The Cancer Genome Atlas.

Gene	HR	P value	Lower	Upper
SDC2	1.198408375	0.002623116	1.065124145	1.348371118
VCAN	1.254174878	0.005099287	1.070349846	1.469570562
ME1	1.207011799	0.019435313	1.030825077	1.413312031
PC	1.199767177	0.030505733	1.017278658	1.414992115
SOX9	0.811543291	0.04660378	0.660673598	0.996865191

HR: Hazard ratio; SDC2: Signature comprising syndecan 2; VCAN: Versican; ME1: Malic enzyme 1; PC: Pyruvate carboxylase; SOX9: SRY-box transcription factor 9.

Table 2 Parameters of five genes from least absolute shrinkage and selection operator regression analysis.

Gene	coef	HR	Low CI	Up CI
SOX9	-0.188516073	0.811543291	0.660673598	0.996865191
PC	0.106222641	1.199767177	1.017278658	1.414992115
SDC2	0.14291442	1.198408375	1.065124145	1.348371118
VCAN	0.147527733	1.254174878	1.070349846	1.469570562
ME1	0.167403756	1.207011799	1.030825077	1.413312031

HR: Hazard ratio; CI: Confidence interval; SDC2: Signature comprising syndecan 2; VCAN: Versican; ME1: Malic enzyme 1; PC: Pyruvate carboxylase; SOX9: SRY-box transcription factor 9.

External validation of the five-gene signature

To test the robustness of risk score, four independent datasets including GSE62254, GSE15459, GSE84437, and GSE26901 were used as an external validation. Similar results to the training set were found in the external validation sets that high-glycolysis patients showed worse outcomes. The area under the receiver operating characteristic curves for the five-gene signature base on four Gene Expression Omnibus databases suggested that our model performed well in predicting OS of GC patients (Figure 4A and B).

Figure 4 External validation of risk model in different glycolysis status. A: Kaplan-Meier survival analysis of different glycolysis status in four external validation data (yellow: High glycolysis status; blue: Low glycolysis status); B: Receiver operating characteristic curves of four external validation data for risk models. Red: 1 year, yellow: 3 years, and blue: 5 years. TPR: True positive rate; FPR: False positive rate; AUC: Area under the receiver operating characteristic curve.

Construction and validation of prognostic nomograms

To explore a clinical method for predicting the 3- and 5-year OS in GC patients, a nomogram including clinical pathological parameters and risk scores was constructed (Figure 5A). This nomogram performed well in predicting the 3- and 5-year survival (Figure 5B-D).

Figure 5 Nomogram combining glycolysis-related risk score and clinicopathological parameters for predicting gastric cancer patients’ survival. A: Construction of a nomogram including age, gender, grade, stage, T, N, M and risk score to predict 3- and 5-year overall survival (OS); B and C: A greement between predicted and observed 3- and 5-year OS. Dashed line indicated ideal performance. The actual performances were indicated by blue line (3-year OS) (B) and red line (5-year OS) (C); D: Merged results. OS: Overall survival.

Differences in glycolysis pathways between high- and low-glycolysis patients

We identified functional differences among patients with differently glycolysis status by conducting GSVA. The high-glycolysis patients had greater dysregulation of pancreas beta cell and estrogen pathways. Low-glycolysis patients had greater Myc target V1, oxidative phosphorylation, Myc target V2, mechanistic target of rapamycin complex 1 signaling and G2M checkpoint pathways (Figure 6).

Figure 6 Gene-set variation analysis and correlation between different glycolysis status in gastric cancer. Blue: High glycolysis status group. Green: Low glycolysis status group.

Differences in tumor-infiltrating immune cells and multiomics content between high- and low-glycolysis patients

It has been shown that glycolytic differences between the tumor and immune cells in the TME affect the immunoregulatory properties of the tumor[28]. Therefore, we analyzed tumor immune infiltration under different glycolytic states in GC. High glycolysis state was more relevant to resting of dendritic cells (DCs), mast cells, CD4 memory T cells, monocytes, eosinophils and M2 macrophages. Hypoglycolysis was significantly correlated with M1 macrophages, CD4 memory T cells, DCs, T follicular helper (Tfh) cells and natural killer (NK) cells (Figure 7A).

Figure 7 Correlation between different glycolysis status and immune cell infiltration and multiomics analysis including tumor mutational burden, neoantigen, microsatellite instability and drug sensitivity in patients with different glycolysis status. A: Different immune cell infiltration between different glycolysis status; B: Tumor mutational burden, neoantigen, microsatellite instability between patients with different glycolysis status; C: Drug sensitivity between patients with different glycolysis status. P < 0.05 indicates statistical significance. TMB: Tumor mutational burden; MSI: Microsatellite instability.

The multiomics contents including TMB, neoantigens, MSI and drug sensitivity were analyzed in high- and low-glycolysis patients. Patients with high glycolysis had lower TMB, neoantigen load and higher incidence of MSI compared with low glycolysis patients (Figure 7B). It is worth noting that pRRophetic algorithm prediction of IC₅₀s of chemotherapy drugs (docetaxel, cisplatin, paclitaxel, doxorubicin, cytarabine, bleomycin, gemcitabine, methotrexate, mitomycin C and camptothecin) found that eight drugs had higher IC₅₀ in high glycolysis patients. This indicated that high glycolysis patients were less sensitive to these drugs (Figure 7C).

Correlation of risk of glycolysis status and clinical parameters

We analyzed the correlations between different clinical parameters and the glycolysis score for each GC patients in TCGA-STAD. The results suggested that patients with grade 2/3 GC or poor prognosis (dead) had high glycolysis status (Figure 8).

Figure 8 Correlation between patients with different clinicopathological parameters and different glycolysis status. Gender (0: Male, 1: Female), fustat (0: Live, 1: Dead). P < 0.05 indicates statistical significance. A: Age; B: Survival state; C: Gender; D: Pathological grade; E: M stage; F: N stage; G: Tumor-node-metastasis stage; H: T stage.

DISCUSSION

As a part of metabolomics, glycolysis is closely involved in carcinogenesis of GC. The large glucose consumption, high activity of glycolysis and continuous accumulation of lactic acid caused by glycolysis change the TME in GC, affecting the metabolism, growth, invasion and metastasis of GC cells[11,29]. Compared with genomics and proteomics, the metabolomics such as glycolysis is the last step before phenotype. This shows that it is important for GC phenotypic control to study the mechanism of glycolysis and precisely target the regulation of glycolysis. Considering the poor effectiveness of some biomarkers, such as carcinoembryonic antigen and carbohydrate antigens, more sensitive methods for poor prognosis of GC need to be explored and established. Therefore, glycolysis profile analysis in GC might provide potential biomarkers for early diagnosis and prognosis; in turn providing evidence for identifying drug targets and evaluating treatment response. Based on the above, this study aimed to explore expression characteristics and potential mechanisms of the glycolysis-related genes of GC through in-depth analysis, and evaluate the prognosis of GC by using gene expression status.

By univariate and LASSO regression analysis, five genes, SDC2, VCAN, ME1, PC and SOX9 were selected for evaluation of metabolic status. Abnormal glycolysis mediated by abnormal gene expression may affect the occurrence, development and poor prognosis of GC. Previous data prove that these genes regulate tumor metabolism. ME1, an enzyme converting malate to pyruvate and reducing NADP+ to NADPH, promoted glycolysis by increasing glucose uptake and glycolytic flux and suppression of oxidative phosphorylation in BLBC cells[30]. In nasopharyngeal cancer, downregulated expression of ME1 led to glucose addiction[31]. PC, an anaplerotic enzyme, participated in various cellular metabolic pathways, including gluconeogenesis, amino acid synthesis and glucose-induced insulin secretion. PC also affected the carbon supply of mitochondria from glucose. Many metastatic pathways were closely related to the expression of PC[32]. SOX9, a transcription factor, has also been reported as a critical TME regulator[33]. The involvement of these genes in metabolomics was confirmed in our GO and KEGG analyses. These genes also play roles in the metabolism and prognosis of various tumors[34-37]. Traditional single biomarkers have many limitations for GC prognostic evaluation. For this reason, the five genes screened were integrated to assess the prognosis of GC patients with different glycolytic states. In the training and validation sets, patients with a hyperglycolytic state caused by the five genes had significantly poor prognosis, especially in those with grade 2/3 GC. This was consistent with the results from single gene assessment reported by other researchers. The role of these genes is also critical in the maintenance and promotion of malignant phenotypes. Our GSVA indicated that, in addition to metabolism-related pathways, estrogen pathways, Myc signaling and G2M checkpoint pathways also participated in the anaerobic glycolysis. Similar results have been reported in other studies. For example, syndecan-2 protein promoted epithelial-mesenchymal transition by activating the MAPK pathway in colorectal cancer[38]. VCAN, a chondroitin sulfate proteoglycan in the extracellular matrix, provided antiadhesive and high proliferative activities[39,40]. Similar functions were observed for PC and SOX9 genes[41,42].

In view of the important role of glycolysis and immunity in TME, we analyzed infiltrating immune cells under different glycolysis conditions. High glycolysis state was significantly correlated with immune cell resting and M2 macrophages. Conversely, hypoglycolysis was significantly correlated with M1 macrophages and immune cell activation. These immune cells and TME mutually influence each other[43,44]. Because of a lack of significant stores of nutrients, these immune cells rely heavily on energy supplements, including absorption of glucose and other nutrients from TME to activate and maintain their functions. One method of energy supplementation that can activate immune cells in cancer is upregulation of glycolysis[45]. Unlike the rapid metabolic transitions between glycolysis and oxidative phosphorylation that occur in nonmalignant tissues, immune cells in the malignant tumors are preferred to aerobic glycolysis[45]. In tumor immunity, DCs migrate to lymphatic tissues, where tumor antigens are presented to initial T cells to activate an antitumor response[46]. Although glycolysis is required for the maturation of DCs[47], several studies have indicated that the effect of the acidic TME resulting from glycolysis on the migration and immune function of DCs is one of the important causes of tumor immune escape[48]. This also explains our finding that most mature DCs were in a resting state under hyperglycolysis. For monocytes, previous studies have suggested that they do not inhibit host immune response to tumor cells, and promote tumor angiogenesis, helping cancer cells to avoid clearance by immune cells[49]. We also found that most mast cells were in a resting state under the status of high glycolysis, and did not degranulate and release histamines, inflammatory mediators, cytokines and chemokines, and thus did have an antitumor function[50]. It was notable that three types of immune cells seem to prefer different glycolytic states (eosinophils for hyperglycolysis, Tfh cells and NK cells for hypoglycolysis). The study from Zaynagetdinov et al[51] indicated that in the high glycolysis status, cancer cells enforce the ability of cancer cell metastasis with the help of eosinophilic. Tfh cells are crucial in humoral and cellular immunity by regulating the activity of B and T cells[52,53]. In this study, in the low glycolysis status, a high percentage of NK cells and Tfh cell infiltrations often indicated good prognosis. Similar results were observed in gastric and breast cancer. For instance, low NK cell infiltration was associated with low survival rate and disease worsening in GC[54]. Breast cancer patients with high Tfh cell infiltration in tumor tissues had better prognosis[55]. M1 and M2 are two types of macrophage activation[56]. Although M1 cells preferred to glycolysis[57], our results indicated that patients with hyperglycolysis showed high M2 macrophage levels. Notably, studies suggested that M2 macrophages promoted angiogenesis and tumor development compared with the cytotoxic role of M1 macrophages in tumor immunity[58,59]. To some extent, this shows that M2 macrophages are more inclined to play a carcinogenic role in the high glycolytic status.

Accurate prediction and staging in GC could guide physicians with more treatment strategies to determine prognosis and assess treatment outcomes. According to the vital role of glucose metabolism in GC, we analyzed the multiomics results, including diagnosis and treatment related to glycolysis. The 8^th edition of the American Joint Commission on Cancer TNM staging system suggests a better prognostic stratification for GC patients compared with the 7^th edition, especially for N3 stage GC[60]. However, although staging optimization has largely solved the stratification of prognosis among different GC stages, different outcomes still occur in the same stage. Considering the low predictive sensitivity and specificity from traditional serum markers and single biomarkers, we combined the risk score and some clinicopathological characteristics, including TNM stage, to construct a novel nomogram, which provides a better prediction of prognosis for GC patients. Immune checkpoint inhibitors, especially antibody to programmed cell death 1 or its ligand, have been verified for use in GC treatment[61]. Here, we demonstrated that high glycolysis was negatively correlated with low TMB and neoantigen status in GC patients. TMB has become a useful biomarker in GC for identification of patients that will benefit from immunotherapy[62]. The generation of neoantigens depends heavily on the number of mutations from TMB, promoting T-cell-dependent responses to suppress tumor growth[63,64]. This also explains the consistency between TMB and neoantigens in the same population and indicates that GC patients with low glycolysis might benefit from immunotherapy. Although we did find a high incidence of MSI in high glycolysis patients, a previous study has suggested that the increased glycolysis in the MSI subtype tumor is associated with immunosuppression[65], which indicated that a high glycolysis status weakened the benefits of immunotherapy in MSI-H patients. We also observed some inconsistencies. For instance, Goodman et al[66] indicated that immunoregulation benefited some patients with high TMB. Therefore, the complexity of TME has a great impact on tumor immunotherapy. In addition to TME regulation of immunotherapy, the high glycolysis status also affected the chemotherapeutic drug sensitivity. In our study, many chemotherapeutic drugs, including paclitaxel, also showed low availability in hyperglycolysis status. Further research on the mechanism should be carried out in future.

There were some limitations to this study. First, the risk score and nomogram were derived from a retrospective analysis of the comprehensive data from TCGA and Gene Expression Omnibus databases. Second, external validation lacked our central data result in the limitations of using of the nomogram itself. Third, some mechanisms by which these five prognostic-related genes contributed to abnormal glycolysis of GC need to be verified in vitro and in vivo. Although our study pointed to an association between high glycolysis status and poor prognosis in GC, some important topics still need to be resolved: (1) How to select the preponderant beneficiaries and achieve precision immunotherapy; (2) Why monocytes with elevated glycolysis respond to tumor cells through proinflammatory cytokine production and possibly also through enhanced phagocytosis capacity[67]; and (3) Why some studies reported that high eosinophil counts increased survival in GC[68]. So the TME had the complexity and heterogeneity and multiple factors might affect on the glycolysis and infiltration of immune cells. In this study, we explored how glycolysis affects GC prognosis, and the potential molecular mechanism. The five-gene-mediated mechanism needs to be further studied.

CONCLUSION

In conclusion, high glycolysis status contributed to poor prognosis of GC patients, especially those with grade 2/3 GC and those who had died. The risk score constructed by five genes could distinguish the glycolytic status of GC patients. Patients with hyperglycolysis did not benefit from immunotherapy and chemotherapy. These results provide guidance for the clinical therapy of GC.