Cell Glycolysis Related Gene Signature Predicts Survival of Patients with Colorectal Cancer

1. Introduction

Switching the way to generate energy for cellular processes from mitochondrial oxidative phosphorylation to aerobic glycolysis, termed as “Warburg effect” [1], happens in most of tumor cells. It is characterized by increased glycolysis and lactate production even when the oxygen is available [2]. In colorectal cancer (CRC), glycolysis promotes colonocyte transformation by the obtainment of oncogenic mutation and loss of tumor suppressors, providing apositive feedback to enhance glycolysis and other unique metabolic pathways [3]. Therefore, tumor cell aerobic glycolysis might work as a promising target for cancer therapy. A plethora of potential therapy targets related with aerobic glycolysis have been explored in previous studies. For example, CD36 promotes the ubiquitination of Glypcian 4, initiates β-catenin/c-MYC signaling pathway and then impairs the downstream glycolysis-related genes, thus suppressing the tumorigenesis [4]. Furthermore, in vitro and in vivo experiments showed that, dioscin, a natural steroid saponin, diminished the phosphorylation of Skp2 and enhanced its degradation, eventually inhibited Hexokinase 2 (HK2) and glycolysis, presenting an antitumor effect [5,6].

CRC is the third leading cause of cancer-related deaths globally [6,7]. Despite the continuous efforts in prevention, screening and management, the incidence of CRC was still increased by 38% from 2007 to 2017 [7,8,9,10]. Patients with CRC show diverse prognosis due to the inherent genetic heterogeneity. Therefore, establishing a reliable prediction model for the prognosis of CRC is meaningful for clinicians to choose therapeutic method individually. Several glycolysis-related markers have been demonstrated to be associated with the prognosis of CRC [11,12]. However, a single gene cannot predict the outcome precisely. The alteration in expression level between groups of different phenotypes was the traditional clue to identify tumor biomarkers. However, those genes possessing biological significance rather than altered expressions might be omitted in this way.

In the present study, we generated a multigene signature predicting the prognosis of individual CRC patient, focusing on the glycolysis-related gene sets from the Gene Set Enrichment Analysis (GESA) [13]. Noteworthy, a clear differential gene threshold is not required in GSEA, having its advantages by coordinated differential expression of annotated groups of gene set. We combined GESA and the Molecular Signatures Database (MSigDB) to discover new prognostic markers in patients with CRC. Accordingly, the complete mRNA expression datasets of 540 CRC patients from the Cancer Genome Atlas (TCGA) was explored to identify genes and pathway signatures that could predict the prognosis of CRC. A total of 257 mRNAs significantly related to glycolysis were identified and a five-gene risk profile for the precise prediction of patient outcomes was established, providing a new insight toward individual treatment of CRC.

2. Materials and Methods

2.1. Clinical information and mRNA expression dataset of patients

The data including mRNA expression profile and clinical information of 540 CRC patients were collected from TCGA database (https://cancergenome.nih.gov/). Clinical information including gender, age, tumor histology and clinical stage are presented in Table S1. The 5 glycolysis-related gene sets were downloaded from the Molecular Signatures Database v7.1 (http://www.broadinstitute.org/gsea/msigdb/index.jsp)

2.2. Gene Set Enrichment Analysis

The difference of 5 glycolysis-relate gene set between non-tumorous samples and tumor samples was explored using the GSEA tool (http://www.broadinstitute.org/gsea/index.jsp). Normalized p value (p < 0.05) and normalized enrichment score (NES) were used as evaluation criteria to determine the necessity of further analysis.

2.3. Statistical analysis

Raw data of 257 mRNAs expression profile was transformed by log2 and quantile-normalized. Then we identified genes relative to overall survival (OS) using univariate Cox regression. Specific alterations of candidate genes were analyzed in the cBioPortal database (http://www.cbioportal.org/). Next, multivariate Cox regression analysis was carried out to evaluate the survival risk of patients. We calculated the risk scores for each patient using candidate mRNA signatures and ranked them accordingly (Risk score= expression of gene 1× β1+ expression of gene 2× β2+⋯+ expression of gene n× βn). A total of 540 patients were classified into high-risk group (n =269) and low-risk group (n=271) according to the risk score formula. We used Kaplan–Meier survival curves and the log rank method to estimate the prognostic significance. P value < 0.05 was considered as statistically significant. Statistical analyses were performed using GraphPad Prism8 software and SPSS 20.0 software.

3. Results

3.1. Initial screening of the genes using GSEA

Gene expression profiles and initial pathologic diagnosis of 540 CRC patients were collected from the TCGA database. We searched all glycolysis-related gene sets on the Molecular Signatures Database v7.1 (http://www.broadinstitute.org/gsea/msigdb/index.jsp), namely five different gene sets (GO_GLYCOLYTIC_PROCESS, BIOCARTA_GLYCOLYSIS_PATHWAY, HALLMARK_GLYCOLYSIS, REACTOME_GLYCOLYSIS, KEGG_GLYCOLYSIS_GLUCONEOGENESIS). GSEA was used to explore whether these five glycolysis-related gene sets differ between tumor samples and adjacent non-tumorous samples (Figure S1). Results showed that three gene sets, including GO_GLYCOLYTIC_PROCESS, HALLMARK_GLYCOLYSIS, and REACTOME_GLYCOLYSIS were significantly enriched (p <0.01) and selected for further analysis.

3.2. Identification of glycolysis-related genes associated with survival in CRC patients

First, we employed univariate Cox regression analysis of the 236 genes from three gene sets aforementioned for preliminary screening and found that 6 genes were significantly associated with the survival of patient (p<0.05). Then these 6 genes entered the stepwise multivariate Cox regression analysis to identify the most valuable mRNAs combinations. Eventually, five genes (ENO3, GPC1, P4HA1, IDUA, ANKZF1) with HR > 1 were verified, as shown in Table S2.

Then we analyzed the alteration of these five selected genes in 526 clinical CRC samples from the cBioPortal database, and found that the genes were altered in 33 (6.3%) cases (Figure 1(a)). Among those patients, 0.19% had mutations and 0.38% had deep deletions in ENO3; 0.19% had amplifications and 0.38% had deep deletions in GPC1; 0.19% had mutations and 0,19% had deep deletions in P4HA1; 0.19% had mutations and 0.76% amplifications in IDUA; 0.38% mutations and 0.19% amplifications in ANKZF1(Figure 1(b)). The expressions of these five genes were significantly upregulated in CRC tissues as compared to according adjacent normal tissues (Figure 2).

3.3. A five genes signature for predicting patients’ outcome

A gene-based prognostic model using a multivariate Cox regression analysis was established to evaluate the survival risk of each patient as follows: risk score = (0.4258 x expression level of ENO3) + (0.378 x expression level of GPC1) + (0.6044 x expression level of P4HA1) + (0.3275 x expression level of IDUA) + (0.4799 x expression level of ANKZF1). We calculated the risk scores for each patient using the five-gene signature. 540 patients were assigned to high-risk group (n =269) and low-risk group (n=271) according to their risk scores. The Kaplan-Meier analysis displayed a significant difference in the outcome of the patients between the high-risk group and the low-risk group (log-rank test p <0.001; Figure 3(b)). Results demonstrated that the lower survival outcome of CRC patients was related to higher risk score (Figure 3(a) and 3(c)). Notably, the relative expression level of those five genes (ENO3, GPC1, P4HA1, IDUA, ANKZF1) were remarkably upregulated when the risk scores were increased (Figure 3(d)).

In order to evaluate the accuracy of this prognostic model, ROC curve analysis was performed. The AUC value of the five genes signature was 0.754, which indicated a favorable sensitivity and specificity in predicting survival of CRC patients. The results showed that the diagnostic significance of the five genes signature was better than a single biomarker (Figure 4).

3.4. Risk score derived from the signature of five genes

To evaluate whether the prognostic ability of the five genes markers is independent of other conventional clinical factors including age, sex, TNM stage, T stage, M stage, N stage, personal neoplasm cancer status, and residual tumor, both univariate and multivariate Cox regression analyses were performed in the entire dataset to assess the importance of these five genes signature. We found that the risk score of these five genes was extremely related to patients’ survival. The results of univariate Cox regression analysis showed that risk score (hazard ratio [HR]: 1.723; 95% confidence interval [CI]: 1.416–2.096; P < 0.001), T stage (HR: 3.769; 95% CI: 2.154–6.597; P < 0.001), N stage (HR: 1.867; 95% CI: 1.336–2.610; P < 0.001), M stage (HR: 5.067; 95% CI: 2.790–9.204; P < 0.001) and TNM stage (HR: 2.382; 95% CI: 1.718–3.304; P < 0.001) were associated with patient survival in the univariate analysis. In addition, from the data set above, we identified risk parameters and T stage as independent prognostic indicators, not only in the univariate but also in the multivariate analyses (Table S3). Notably, the results indicated that the risk parameters had a significant value in predicting the prognosis of CRC patients.

3.5. Validation of five-gene signature for survival prediction

The Kaplan–Meier curves showed that TNM stage (III and IV) (P< 0.001), T stage (T3 and T4) (P < 0.05), N stage (N1 and N2) (P < 0.001), M stage (M1) (P < 0.001), personal neoplasm cancer status (with tumor) (P < 0.001), and presence of residual tumor (P < 0.0001) were significantly associated with poor OS (Figure S2). Moreover, the risk score was an independent prognostic marker for CRC patients grouped by age (<=65 years or >65 years), gender (male or female), N stage (N0 or N1-3), M stage (M0 or M1), and personal neoplasm cancer status (tumor free or with tumor), residual tumor (R0 or R1) (Figure 5). However, the risk score failed to work as prognostic indicators independently according to T stage. The patients in the high-risk group had a worse OS than those in the low-risk group in the TNM stage III-IV (P < 0.001), but no difference was observed in the stage I-II subgroup. Similar results were observed for T stage (T 1-2: P =0.874; T 3-4: P < 0.001), indicating that risk parameter may be a more effective prognostic marker in the advanced stage of CRC.

4. Discussion

An increasing list of genes have been identified as molecular biomarkers to predict the prognosis of CRC [14,15]. For example, in the study in the study conducted by Kubota et al, microRNA-31 worked as a favorable biomarker and a promising therapeutic target in CRC patients [16,17]. Hepatitis B X-interacting protein (HBXIP) was overexpressed in CRC tissues, and assumed to be a prognostic predictor and a new therapeutic target for CRC [18]. Moreover, a previous study showed that ITGA6 and ITGB4 were useful biomarkers for the early detection of CRC and acted as prognostic indicators for the survival of CRC patients [19]. However, a single gene can be affected by various factors. Thus, it’s meaningful to build a model based on combined biomarkers to predict cancer outcome. Previous data indicated a combined signature including 4 genes for the prognosis prediction of CRC with lymphatic metastasis [20]. Furthermore, a gene signature for the prediction of early relapse in stage I-III colon cancer was also discovered recently [21,22]. These signatures showed high efficiency, which indicated that a robust gene signature could offer a more accurate prediction than a single gene.

Ample evidence exists that tumor cells could adjust their energy source from oxidative phosphorylation to glycolysis in order to proliferate efficiently even in a low-oxygen condition. However, researches on glycolysis-related biomarkers of cancer are relatively limited. In this study, we used the CRC dataset in TCGA and performed GSEA analysis to compare glycolysis-related genes of CRC tissues and adjacent normal tissues, in order to identify glycolysis-related biomarkers. Furthermore, a combination of 5 genes (ENO3, GPC1, P4HA1, IDUA, ANKZF1) with prognostic value were identified for CRC patients. Comparing with other biomarkers reported previously, the present risk signature can predict the prognosis of CRC patients with a higher efficiency. Concretely, the present risk parameter may be a more effective prognostic marker in the advanced stage of CRC, which probably attributes to limited sample size, or the variability of glycolysis-related genes may appear mainly in the advanced stage of CRC, making our signature more practical in those stage III-IV patients at who hesitated to take conservative or more aggressive treatment. In addition, the prognostic genes showed a good performance in predicting survival, which was better than TNM staging (AUC for 5-year survival: 0.754 vs 0.686, Figure S3).

Among the five selected biomarkers, β-enolase (ENO3) was involved in the subpathway that synthesizes pyruvate from D-glyceraldehyde 3-phosphate, acting as a lyase in glycolysis. Mutation in ENO3 caused exercise intolerance and rhabdomyolysis [23]. Moreover, the expression of ENO3 was induced as a consequence of serine/threonine kinase 11 (STK11) loss-of-function in lung cancer. STK11 is a major inactive tumor suppressor in 30% lung adenocarcinoma, therefore, ENO3 may be a promising therapeutic target for STK11-mutant lung cancer patients [24]. Glypican-1 (GPC1) is a proteoglycan localizing mainly at the outside surface of cell membrane and can shed into the extracellular environment, regulating vital signal pathways such as Wnts, Hedgehogs, et al [25]. It has already been chosen as one of the genes in the glycolysis-related signature for predicting survival in patients with lung adenocarcinoma or endometrial cancer [26,27]. As to CRC, Jian Li and his team suggested that increased level of GPC1-positive exosomes in plasma promoted cancer progress. It also indicated a poor prognosis and might act as a specific diagnostic biomarker and therapeutic target [28]. Prolyl 4-hydroxylase subunit alpha-1 (P4HA1), encoding the catalytic subunit of collagen prolyl 4-hydroxylases [29], has been studied in a list of cancers, including melanoma, pancreatic cancer, CRC, et al. Consistent with our results, P4HA1 was elevated in CRC tissues, and also promoted tumor growth and metastasis in mouse models. Furthermore, diethyl-pythiDC, a small molecular inhibitor targeting P4HA1, has been proved to be effective in CRC patient-derived xenografts models [30,31].

The researches mentioned above imply that several genes involved in the glycolysis-related gene signature of CRC have already been studied in laboratories and clinical specimens, showing potential for targeted therapy. Whereas the others remain poorly understood. For instance, Alpha-L-iduronidase (IDUA) belongs to the glycoside hydrolase 39 family, and its mutation may cause type I mucopolysaccharidosis, a rare lysosomal storage disease [32,33]. However, the role of IDUA in cancer has rarely been explored. Ankyrin repeat and zinc-finger domain-containing 1 (ANKZF1) was reported to play an integral role in the cellular response to hydrogen peroxide and in the maintenance of mitochondrial integrity under cellular stress conditions [34]. A recent study showed that high expression of ANKZF1 was associated with a poor overall survival and recurrence-free survival in CRC by in silico analysis [35]. The role of IDUA and ANKZF1 in CRC cellular function awaits further proof, as well as the correlation with clinical characteristics.

There are some limitations in this study. Firstly, as the data were derived from public databases, it’s difficult to assess their quality effectively. Secondly, due to bioinformatics analysis in the current study, biological experiments are needed for further validation. In conclusion, we reported a five-gene risk signature related to cellular glycolysis, which could predict the prognosis of CRC patients, and a high-risk parameter indicates a poor prognosis. These findings provide insights into the mechanisms of cellular glycolysis and indicate new potential glycolysis-related biomarkers for the diagnosis and targeted therapy of CRC.

5. Conclusions

A five-gene signature associated with cellular glycolysis for predicting the outcome of CRC patients was generated, which was a valuable prognosis model with high efficiency, providing new perspectives in targeted therapy. These findings could provide insight into the identification of patients with poor prognosis in CRC.