A Novel Glycolysis and Hypoxia Combined Gene Signature Predicts the Prognosis and Affects Immune Infiltration of Patients with Colon Cancer

Introduction

The prevalence of colon cancer ranks the third worldwide.¹ Compared with patients diagnosed at stage I or II, the 5-year survival rate of those at stage IV declines to 12%, with a decrease of 86.8% and 85.4%, respectively.² The prognosis of colon cancer patients varies at different stages. More efforts are needed to find a reliable and accurate prognostic model or staging system to assess the prognosis and specify treatment strategies of patients with colon cancer.

Emerging evidence indicated that colon cancer should be considered primarily as a metabolic disease involving rewiring of cellular metabolism, a characteristic of most cancer types.³ Cancer cells exhibit metabolic profile, including extensive consumption of glucose and glutamine and extensive glycolysis⁴ and the glycolysis was a preferential choice for energy generation instead of oxidative phosphorylation in tumors, which is known as the “Warburg effect”.⁵ Many more metabolic pathways observed in cancer cells often include glucose, glutamine and lipid metabolism.⁶ Furthermore, hypoxia is a significant cellular metabolic process in neoplasia related to the Warburg effect and the clinical prognosis of cancer patients.^7,8 As a hallmark of tumor microenvironment, hypoxia changes expression profile and affects immune response, which has been noticed to impact immunotherapy and seen as potential targets in cancer treatment.^9,10 Besides, oncogenes and tumor suppressors have been discovered in these reprogrammed energy metabolisms, and these cancer-associated metabolic changes provide cancer cells with energy for initiation, progression and spreading.¹¹

Here, a glycolysis and hypoxia combined gene signature was built up for a more accurate prediction of the prognosis and evaluation of immune infiltration in colon cancer patients. Furthermore, a nomogram was established combining this gene signature with clinicopathological factors for more comprehensive and accurate prognosis prediction.

Materials and Methods

Collection of Clinical Parameters and Gene Expression

The Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/) is a powerful data-driven platform that provides researchers and bioinformaticians with access to search and download cancer data for analysis. We sought and obtained the clinical parameters and RNA expression profiles from The Cancer Genome Atlas (TCGA) database on this platform. The main clinical parameters include gender, age at diagnosis, tumor stage, follow-up time and survival status. Finally, we included 417 samples (378 colon cancer and 39 normal tissues) for our analysis. The general clinicopathological information is detailed in Table 1. The procedure of the analyses in this study is shown in Supplementary Figure 1.

Table 1 Summary of Baseline Clinicopathological Features of Patients with Colon Cancer in This Research

Gene Set Enrichment Analysis (GSEA)

Genes that possess similar chromosomal localization, functions, pathways or other features can be collected into one gene set. GSEA was performed to select glycolysis and hypoxia gene sets enriched among cases of colon cancer patients by GSEA software.¹² Before GSEA analysis, we searched the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb) and obtained five predefined glycolysis-related gene sets (BioCarta glycolysis, hallmark glycolysis, GO glycolytic, reactome glycolysis gene sets and KEGG glycolysis) and six hypoxia-related gene sets (Buffa hypoxia, Harris hypoxia, Mizukami hypoxia up, Mizukami hypoxia down, reactome cellular response to hypoxia, and Winter hypoxia gene sets). For each analysis, gene set files, gene expression, and phenotype labels files were loaded in GSEA software. To reach normalized enrichment score (NES), we conducted the gene set permutations were for 1000 times. The minimum and maximum criteria for gene sets selection were 15 and 500 genes, respectively. Finally, the gene sets with |NES| > 1.6, P < 0.05 and false discovery rate (FDR) <0.1 were proceeded for further subsequent analysis.

Construction of Glycolysis and Hypoxia Gene Signature

To sort out the survival-related glycolysis and hypoxia genes, univariate Cox regression analysis was firstly initiated among the identified gene sets followed by the multivariate Cox regression analysis when P value <0.1 observed in the univariate Cox analysis. Then, using a linear combination of gene expression levels, a novel gene signature associated with glycolysis and hypoxia was established to predict prognosis. The risk score is defined as f (s, Z) = f0(s) × exp (expression of gene1 × α1 + expression of gene2 × α2 + ⋯ + expression of genen × αn). Z represents gene expression of each case in the formula and f0(s) is a constant, which was confirmed in the multivariate Cox regression analysis. Besides, number n refers to the number of genes with independent prognostic value and αn represents the corresponding regression coefficient. Exp () denotes an exponential function with the base e. Risks of all patients were calculated, after which they were divided into two risk categories with the median risk value.

Comparisons of Immune Cell Infiltration and Checkpoint Genes Expression Between Two Groups

Immune responses and immunocytes’ infiltration were analyzed within high- and low-risk groups by the single sample gene set enrichment analysis (ssGSEA) with the “GSVA” R package.^13,14 General immune status of each patient was assessed by ESTIMATE score, stromal score, immune score and tumor purity score. Additionally, in each colon cancer sample, we utilized CIBERSORT algorithm and depicted the profile of 22 immunocyte subtypes.¹⁵ Furthermore, expressions of checkpoint genes were compared, such as CD276 and CTLA4, between the high- and low-risk groups to explore the potential significance of this multiple genes combined model in immunotherapy response evaluation.

Establishment and Assessment of the Nomogram

Incorporating clinicopathological factors with the gene signature, a comprehensive nomogram was constructed using the “rms” package of R software to predict the prognosis of colon cancer patients. Harrell’s C-index was then calculated to estimate the overall survival discrimination ability of the nomogram model. A higher Harrell’s C-index demonstrates a better discrimination ability. Besides, the precision of the nomogram model was evaluated through ROC curves and calibration plots.

Statistical Analysis

The prognostic genes were sorted by univariate and multivariate Cox analyses, which were further evaluated by Kaplan–Meier analysis and Log rank test. Student’s t-test was used for comparison of continuous variables in this study. R software was utilized to construct the heatmaps, forest plots, nomogram plot, ROC curves, and calibration plots. A P-value of less than 0.05 is statistically significant. All statistical analyses were conducted by R software (v.4.0.2) and Excel software (Microsoft).

Results

Identification of Gene Sets Enriched in Colon Cancer by GSEA

In colon cancer cases, three glycolysis and two hypoxia gene sets were enriched significantly through GSEA, including GO glycolytic process (NES = 1.91, P = 0.001, FDR = 0.002), hallmark glycolysis (NES = 1.71, P = 0.005, FDR = 0.01), reactome glycolysis (NES = 2.00, P < 0.001, FDR < 0.001), Buffa hypoxia (NES = 1.86, P = 0.004, FDR = 0.008) and reactome cellular response to hypoxia gene sets (NES = 1.71, P = 0.026, FDR = 0.05) (Figure 1).

Figure 1 Enrichment plots of five gene sets between colon cancer and normal tissues performed by GSEA. (A) GO glycolytic process, (B) hallmark glycolysis, (C) reactome glycolysis, (D) Buffa hypoxia and (E) reactome cellular response to hypoxia gene sets.

Identification of Prognostic Glycolysis and Hypoxia Genes in Colon Cancer Patients

Based on gene expression analysis, we found that expressions of 95 genes were downregulated and that of 225 genes were upregulated in colon cancer cases when compared to normal ones (Supplementary Figure 2). Through univariate and multivariate Cox regression analysis, ENO3, STC2, SPAG4, GPC1, and P4HA1 genes were identified as independent indicators associated with the prognosis of colon cancers (Figure 2; Table 2). The risk of each patient was calculated by a signature model composed of the glycolysis and hypoxia genes obtained previously. Based on the median risk value, the whole cohort was stratified into high- or low-risk groups for further validation of the five-gene signature. Then, a heatmap exhibited the expression of the five genes between the two groups (Figure 3). As shown in Figure 4A and B, the distribution of survival status and risk scores of colon cancer patients indicated poorer prognosis of patients with higher risks. And the area under the ROC curve was 0.750 which indicated a good accuracy for prognosis prediction (Figure 4C). It was further confirmed by Log rank test and Kaplan–Meier survival analysis that patients with higher risks were likely to have significantly poorer prognoses (P = 0.002, Figure 4D).

Table 2 The Information of Five Prognostic Glycolysis Genes Associated with Overall Survival in Patients with Colon Cancer

Figure 2 Five identified gene expression profiles between colon cancer and normal tissues. (A) P4HA1, (B) ENO3, (C) SPAG4, (D) STC2 and (E) GPC1. *** Indicates P-value <0.001.

Figure 3 Heat map of the five genes that were used to construct the risk score model.

Figure 4 The glycolysis gene signature related to risk score predicts the overall survival of colon cancer patients. (A and B) Patient risk score distribution and survival status. (C) The prognostic power of the gene signature was evaluated using the AUC value of the ROC curve. (D) Kaplan–Meier curves of overall survival in high- or low- risk groups.

Performance of the Prognostic Gene Signature in Subgroups by Age, Gender and AJCC Stage

Furthermore, prognostic factors were selected by univariate and multivariate Cox analyses in colon cancer patients. The univariate Cox analysis concluded that risk score (P < 0.001, HR = 1.280, 95% confidence interval (CI): 1.178–1.391), age (P = 0.018, HR = 1.030, 95% CI: 1.005–1.055), and AJCC stage (P < 0.001, HR = 2.502, 95% CI: 1.851–3.382) were significantly associated with the overall survival of colon cancer (Figure 5A). Multivariate Cox regression showed that risk score (P < 0.001, HR = 1.263, 95% CI: 1.155–1.382), age (P < 0.001, HR = 1.045, 95% CI: 1.019–1.071), and AJCC stage (P < 0.001, HR = 2.680, 95% CI: 1.952–3.680) were independent risk factors of overall survival (Figure 5B). Through the Log rank tests and Kaplan–Meier survival analyses, we found that patients over 65 years old and at stage III–IV had significantly poorer prognoses than the control groups (Figure 5C–E).

Figure 5 Risk score based on the gene signature acts as an independent prognostic marker. (A and B) Univariable and multivariable analysis of risk score and various clinical factors. (C–E) Kaplan–Meier curves survival analysis in different groups classified by age, gender, and AJCC stage.

Then, we classified all patients into different categories according to age (≤65 versus >65 years), gender (female versus male), and AJCC stage (stage I–II versus stage III–IV), and by the median risk value, they were then classified into high- or low-risk groups. High scores indicated worse prognosis in the elderly (>65 years old) subgroup (P = 0.003; Figure 6B), female subgroup (P = 0.007; Figure 6C), and stage I–II subgroup (P = 0.037; Figure 6E), but no striking difference in prognosis was discovered between high and low-risk group in the young (≤65 years old) (P = 0.297; Figure 6A), male (P = 0.168; Figure 6D) and stage III–IV (P = 0.077; Figure 6F) subgroups.

Figure 6 Kaplan–Meier curves for the prognostic value of risk-score signature for the patients with colon cancer stratified by various clinical features in the TCGA dataset. (A and B) Age. (C and D) Gender. (E and F) AJCC stage.

Status of Immune Infiltration and Checkpoint Genes Between High- and Low-Risk Groups

From the results of ssGSEA, immune responses, pathways and immunocytes infiltration were analyzed between the two risk groups and presented in a heat map (Figure 7A). Besides, to assess the general immune status of each case, we calculated the ESTIMATE score, stromal score, immune score and tumor purity score with the “estimate” package by R software, which was also shown in the heatmap. Furthermore, weighed by respectively specific expression profile, 22 immunocyte subtypes in colon cancers were speculated by the CIBERSORT algorithm (Figure 7B). In high-risk group, the proportion of NK cells resting significantly increased (P < 0.05) while that of monocytes (P < 0.01), dendritic cells resting (P < 0.01) and dendritic cells activated (P < 0.05) significantly decreased compared to low-risk group.

Figure 7 Immune status and immunocyte subtypes in high- and low-risk groups. (A) Immune responses, pathways and general immune status in all cases of colon cancer. (B) Immunocyte subtypes differences between high- and low-risk groups. *indicates P-value < 0.05, ** indicates P-value < 0.01.

Additionally, expression profiles of acknowledged immune checkpoint genes were compared between the two groups and altogether 22 differentially expressed checkpoint genes were observed (P < 0.05; Figure 8A). Most of them tended to overexpress in the high-risk group. Then, we explored the relationship between the prognostic gene signature and checkpoint genes, and the expression of most checkpoint genes were positively correlated with estimated risk scores (P < 0.05). We depicted part of the conclusion in scatter plots (Figure 8B) and the rest of them were shown in Supplementary Figure 3.

Figure 8 Relationship between immune checkpoint genes and the gene signature. (A) Differential expression of checkpoint genes between high- and low-risk groups. (B) Scatter plots of the checkpoint genes and signature-based risk scores. *indicates P-value < 0.05, ** indicates P-value < 0.01. *** indicates P-value < 0.001.

Construction of the Nomogram Model Based on Glycolysis and Hypoxia Gene Signature and Other Clinicopathological Parameters

Based on the gene signature and independent prognostic parameters, an integrative nomogram was built to provide a new tool to predict the survival probability for colon cancer patients (Figure 9A). The survival probability of each individual can be simply calculated and intuitively presented through the cumulative score of all identified values (Supplementary Table 1). Patients were assigned a nomogram risk score and separated into high and low-risk groups. Patients at higher risks were validated to have significantly worse prognoses (Figure 9B). The area under the ROC curves of 3- and 5-year survival probability prediction were 0.796 and 0.803, respectively, showing the good performance of the nomogram in prognosis prediction of colon cancer (Figure 9C and D). Besides, by calculating the Harrell C-index, the predictive accuracy of the nomogram model was measured and the C-index is 0.771 (95% CI: 0.712–0.830) in the internal verification of the nomogram, which indicates good discrimination ability. Calibration plots of 3- or 5-year survival probability prediction also showed the good performance of the nomogram in prognosis prediction of colon cancer patients (Figure 9E and F).

Figure 9 The nomogram for 3- and 5-year overall survival (OS) prediction. (A) Nomogram combining the gene signature and clinical features to predict the 3- and 5-year OS of colon cancer patients. (B) Kaplan–Meier curve of OS in high- or low-risk groups. (C) ROC curve for 3-year OS prediction. (D) ROC curve for 5-year OS prediction. (E) Calibration curve of the nomogram for the prediction of 3-year OS. (F) Calibration curve of the nomogram for the prediction of 5-year OS.

Discussion

It is a distinguishing characteristic in cancer cells that they exhibit remarkably increased metabolic requirements in comparison to normal cells. To maintain their rapid proliferation, cancer cells behaved more glucose consumption and lactic acid production rather than catabolizing glucose via the tricarboxylic acid cycle.^16,17 Along with the aggressive growth, oxygen levels drop down to hypoxic levels in the microenvironment and more cells in return change the metabolic pathway towards glycolysis for essential energy supply.^9,18 Glycolytic intermediates under hypoxia provide continuous energy and nutrients to support uncontrolled proliferation, tumor invasion, and metastasis, and eventually contribute to tumor growth and survival.^19,20 Besides, glycolytic and hypoxic pathway metabolites are widely acknowledged to promote immune cell activation and inflammatory infiltration in tumors.²¹ Thus, it is the glycolysis and hypoxia status that may serve as emerging hallmarks in cancer proliferation, infiltration and metastasis and a potential biomarker for prognosis prediction of cancers. Importantly, identifying glycolysis and hypoxia-related genes may have prognostic value for patients with colon cancer.

Increasing evidence shows that more information is needed to accurately predict prognosis, not just traditional clinicopathological factors. Thus, tumor markers that can accurately predict the oncological outcomes when incorporated with clinicopathological features are essential for predicting the outcomes of colon cancers.^22–24 More and more evidences show that glycolysis and hypoxia-related genes have good performance in predicting the clinical prognosis of various kinds of solid tumors. For example, a four-gene signature derived from the REACTOME glycolysis gene set had a significant correlation with metastasis and survival of bladder cancer patients.²⁵ High glycolysis scores integrating 10 glycolytic genes profiles predicted poor prognosis and poor chemotherapy response, closely associated with IDH1 mutation in glioblastoma.²⁶ Another four-gene signature based on hypoxia and immune status performed promising prognostic assessing ability in gastric cancer²⁷ and similar hypoxia-immune prognostic signature was developed and validated in head and neck cancer.²⁸ Here, we provide evidence that glycolysis and hypoxia gene signature is an independently effective method for prognostic assessment strategy for colon cancer patients.

In this research, we established a five-gene signature (ENO3, STC2, SPAG4, GPC, and P4HA1) using data from the TCGA database to predict prognosis of colon cancer. Risk of each patient was calculated and survival analysis revealed that patients with higher predicted risk had worse clinical outcome. In further subgroup analysis, the gene signature could classify patients with different prognosis in elderly (>65 years old), female and stage I-II subgroups. However, in the other subgroups, low-risk patients showed the tendency of better overall survival without statistical significance, which might result from insufficient patient cohort and limited follow-up time. Besides, patients under 65 years old or at stage III-IV were more likely to have similar good or bad survival rates that this risk model had difficulty in distinguishing them. As for the survival similarity between the two subgroups in male patients, more investigation is needed. Anyway, the five-gene signature is more suitable for survival prediction of the whole patient cohort. Immune responses, infiltration and checkpoint genes were explored in patients with different risk scores, revealing the impact of glycolysis and hypoxia signature on immune status and indicating more potential targets for immunotherapy. Additionally, a nomogram was built up by integrating the gene signature with clinicopathological parameters (age, gender, and AJCC stage) to efficiently predict the clinical outcome in patients with colon cancer.

ENO is one of the key enzymes in glycolysis.^29,30 It can catalyze the conversion of 2-phosphoglycerate to the high-energy compound phosphoenolpyruvate, which is present in almost all adult tissues.³¹ ENO3, also known as γ-enolase, is acknowledged to be upregulated in multiple cancer cases in humans, such as gastric cancer, glioma, breast cancer, head and neck cancer, lung cancer, endometrial cancer, pancreatic adenocarcinoma, and non-Hodgkin’s lymphoma.³² The studies showed that glucose metabolism significantly decreased after knocking down ENO3 and it was validated strongly that in lung cancer cells, loss-of-function mutation of the gene STK11 directly stimulated the overexpression of ENO3.³⁰ Functional annotation also revealed that ENO3 was enriched in the hypoxia-inducible factor 1 pathway.³³ It was known that ENO3 is prognostic and high expression is unfavorable in colorectal cancer. In accordance, it was also noticed that high expression of ENO3 had a close correlation with tumor metastasis, indicating that patients with colon cancer have a poor prognosis.

STC2, a secreted glycoprotein hormone of as yet unknown functions, is thought to be involved in an autocrine or a paracrine pattern.³⁴ Its expression was demonstrated to be induced by estrogen, retinoic acid, progesterone, hypoxia, and/or unfolded protein responses and repressed by anti-estrogen in many human cancer cell lines.^35,36 Under hypoxic conditions, the proliferation and expansion of cells in breast and ovarian cancers were stimulated by the overexpression of gene STC2, showing that STC2 was likely a HIF-1 target gene.³⁷ However, underlying crosstalk and mechanism of STC2 under hypoxic conditions in colorectal cancer has not been fully understood. In our study, we found that higher expression of STC2 was significantly related to a worse prognosis of colon cancer patients, addressing its close relationship with colon cancer cells.

SPAG4, a member of the SUN (Sad-1, UNC-84) protein family, was originally noticed to be a testicular specific gene, participating in sperm tail development, nuclear membrane integrity maintenance and nuclear remodeling.³⁸ Additionally, recent researches had also illustrated that SPAG4 could act as a new tumor marker, which affects the biological behavior of tumors.^39,40 It was found that SPAG4 expression was upregulated in various cancers in humans, including lung cancer,³⁹ renal cell carcinoma⁴¹ and glioblastoma.⁴² So far, studies had hardly reported the SPAG4 expression and its effect on prognosis in colon cancer. The underlying molecular mechanisms played by SPAG4 in tumors are still poorly understood. Given the current evidence, it may appear to include the following actions. Ji et al confirmed that SPAG4 and Nesprin3 interaction have fundamental pathological functions in the process of lung carcinoma cell migration.³⁹ The results of the study in glioblastoma indicated that by repressing the MEK/ERK signaling pathway activation, SPAG4 knockout remarkably suppressed the growth of brain tumor cells including proliferation and migration.⁴² Our research unraveled the independent prognostic significance of gene SPAG4 in colon cancer. Combined with these studies, we further confirmed that glycolysis and hypoxia-related SPAG4 gene risk score model can be accurate in predicting prognosis.

GPC1 was overexpressed in various cancers in humans, including breast cancer,⁴³ pancreatic cancer,^44,45 and glioma.⁴⁶ GPC1 participates in tumor growth and angiogenesis by activating FGF-FGFR signaling.^47–49 Li J et al found that GPC1 overexpression significantly suppressed E-cadherin and increased vimentin protein expression in colorectal cancer cells, GPC1 gene overexpression significantly upregulated Snail protein expression in colorectal cancer cell lines.⁵⁰ Furthermore, overexpression of the GPC1 gene promoted cell invasion, migration, and epithelial-mesenchymal transition. In contrast, silencing of GPC1 gene expression had the opposite effect. Consistently, our data also showed a high GPC1 expression significantly related to a poorer prognosis in patients with colon cancer.

P4HA1 is involved in the formation of proline 4-hydroxylase, maintains the stability of collagen, causes changes in the cancer microenvironment, which therefore participates in the occurrence and metastasis of tumors. In recent decades, the correlation between P4HA1 and different tumors has attracted people’s attention. It was found that hypoxia-inducible factor 1 was widely involved in ECM remodeling by regulating P4HA1 and P4HA2, thereby regulating the occurrence and development of tumors.⁵¹ A new study indicated that miR-122 inhibits the epithelial-mesenchymal transformation process by targeting P4HA1, thereby regulating the invasion and abdominal metastasis of ovarian cancer.⁵² Higher P4HA1 expression was found in some tumor tissues compared with nontumor tissues. Subgroup analyses proved that in normal tissue group, expression of gene P4HA1 was lower than those with stage I in all three types of tumors,⁵³ suggesting that P4HA1 could be an appropriate biomarker for early diagnosis of these cancers. Combined with our findings, P4HA1 could function as a key indicator in monitoring prognosis of patients with colon cancer.

Conclusion

This research constructed a five-gene (ENO3, STC2, SPAG4, GPC1, and P4HA1) signature and a comprehensive nomogram model integrating clinical prognostic features (age, gender and AJCC stage) and this gene signature to predict overall survival and the immune microenvironment of colon cancer. These findings could provide an innovative method in predicting the prognosis of colon cancer patients and specifying treatment strategies.