The expressed gene RNA-sequencing dataset (fragments per kilobase of transcript per million [FPKM] value) and corresponding TCGA clinical information of 473 COAD and 71 non-tumor tissue samples were extracted¹. A total of 298 GRGs, of which 226 were DE-GRGs, were identified by gene set enrichment analysis (GSEA) of glycolysis-related pathways (Supplementary Figure 1).

Differential glycolysis-related genes in TCGA COAD tumor tissues and adjacent non-tumor tissues with a false discovery rate (FDR) < 0.05 were identified using the limma R package. The 226 DE-GRGs identified were displayed by Venn diagram (Supplementary Table 4). Biological process (BP) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed by the clusterProfiler R package (Liao et al., 2019) to identify the functionally enriched genes and classify the gene clusters. A q-value < 0.05 was considered statistically significant.

The PPI of the 226 DE-GRGs was analyzed in the String online database² as previously described (Szklarczyk et al., 2019). The network was visualized by Cytoscape 3.7.2 software³. The key modules in the PPI network were filtered according to the score and the criterion of node counts >5 (Bader and Hogue, 2003) using the MCODE plug-in in Cytoscape. A p-value ≤ 0.05 denoted a significant difference.

Key GRGs identified using the survival R package were further subjected to univariate Cox regression analysis to screen and identify significant prognosis candidates. The candidate genes were incorporated into the multivariate Cox risk model. Eight GRGs associated with significant prognosis potential were retained in the process of multiple calculations. The risk score for each patient was calculated as follows:

Colon adenocarcinoma patients in the training group were categorized into the low- and high-risk groups based on their median risk score. The OS rates of the two subgroups were compared by the Kaplan–Meier method and the log-rank test. A receiver operating characteristic (ROC) curve was constructed by SurvivalROC R package to evaluate the prognosis ability of the aforementioned model. In addition, 177 COAD patient samples from the GSE17536 dataset (Li et al., 2018) were used to validate the predictive ability of this prognostic model. Finally, the RMS-R package was used to construct a nomogram with calibration plots to show the likelihood OS rate.

A total of 43 pairs of colon tumor tissues and adjacent non-tumor tissues were collected from patients who underwent colon cancer surgery at Shanghai Tenth People’s Hospital from September 2019 to June 2020. After surgical resection, the colon cancer tissue samples were immediately cleaned in normal saline and cryopreserved. Prior to the surgery, none of the patients had received radiation and chemotherapy. This study was sanctioned by the Institutional Research Ethical Committee (IREC) of Shanghai Tenth People’s Hospital of Tongji University of Medicine, Shanghai. All patients provided their verbal and written consent.

Total RNA from the colon cancer tissue and CRC cell lines was extracted using TRIzol reagent (Invitrogen, Thermo Scientific, Shanghai, China). Total mRNA was reverse-transcribed into cDNA using the prime script Reverse Transcription reagent kit (TaKaRa Bio, Shiga, Japan). Relative expression of eight potentially prognostic GRG markers was determined by qRT-PCR using the SYBR Green reagent kit (TaKaRa Bio) in an ABI 7500 PCR system (Applied Biosystems). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used to normalize gene expression. Primer sequences are listed in Supplementary Table 1.

Paraffin-embedded colon cancer tissue was cut into 5-μm-thin sections. The sections were deparaffinized and rehydrated. Antigen retrieval was performed by exposure of the sections to 3% hydrogen peroxide along with normal goat serum. Antigen retrieval tissue sections were initially treated with antibody to stanniocalcin-2 (STC2, 1:200 dilution), prolyl 4-hydroxylase subunit alpha 1 (P4HA1, 1:200), protein phosphatase 2 catalytic subunit beta (PPP2CB, 1:200 or 1:250), peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PPARGC1A, 1:250), enolase 3 (ENO3, 1:200), chondroitin polymerizing factor 2 (CHPF2, 1:500), phosphoglucomutase 2 (PGM2, 1:200), and phosphomannomutase 2 (PMM2, 1:200) (all from Abcam, Cambridge, United Kingdom). The sections were subsequently treated with biotinylated secondary rabbit antibodies conjugated with streptavidin–horseradish peroxidase (1:200). The binding of secondary antibody was revealed by reaction with added 3,3’-diaminobenzidine (DAB) followed by counterstaining with hematoxylin. The stained slides were visualized and photographed by optical microscopy (Leica, Tokyo, Japan). Image-Pro Plus 6.0 Software (Media Cybernetics, Rockville, MD, United States) was used to analyze protein expression in the stained tissue sections.

HCT116 and SW480 human CRC cell lines were cultured at 37°C in a 5% CO₂ incubator in Roswell Park Memorial Institute (RPMI) 1,640 medium or Dulbecco’s modified Eagle’s medium (both from Gibco, Grand Island, NY, United States) supplemented with 10% fetal bovine serum (FBS). For transient transfection, prior to the day of transfection 2.5 × 10⁵/ml of HTC116 and SW480 cells were seeded in six-well plates and incubated overnight to allow development of 40 to 60% confluent growth. These cells were transfected with small interfering RNA (siRNA) to STC2 (Gene Pharma, Shanghai, China) using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, United States) and incubated further for 48 h. Target sequences for siRNAs were as follows:

Proliferation of STC2-siRNA-treated HCT116 and SW480 cells was determined using the Cell Counting Kit-8 (Dojindo Molecular Technologies, Inc., Rockville, MD, United States) as previously described (Liu et al., 2019). A clonogenic assay was performed to measure and the monitor the colony-forming capabilities of STC2-siRNA CRC cell (Zhen et al., 2019).

Matrigel precoated Transwell inserts were placed into 24-well plate filled with 500 μl of culture medium supplemented with 10% FBS. Aliquots (400 μl containing 5 × 104 SW480 or HCT116 cells) of STC2-siRNA-treated cell suspensions were added to the upper chamber of Transwell units and incubated at 37°C for 20–24 h in a 5% CO₂ incubator. Cells that traversed the membrane separating the upper and lower chambers of each unit (i.e., invading cells) were fixed with 4% paraformaldehyde (PFA) for 15 min at room temperature and stained with 0.1% crystal violet. Cells that did not invade were gently removed from the membrane surface exposed to the upper chamber using a wet cotton swab. Each Transwell unit was allowed to dry and examined by inverted microscopy (Leica). Quantitative analysis of invaded cell was performed using ImageJ software (NIH, Bethesda, MD, United States).

For the cell migration assay, SW480 and HCT116 cells were transfected with STC2-siRNA and cultured on solid medium. Twenty-four hours later, a wound was created horizontally by gently pressing a 200-μl sterile micropipette tip on the surface of the confluent cell monolayer on six-well plates. Each well was washed twice with serum-free medium. Photographs were taken with an inverted microscope (Leica) at 0 and 24 h. The cell migration rate was determined by time-lapse analysis using ImageJ software.

The CIBERSORT⁴ online bioinformatics analytical tool was used to estimate and distinguish 22 commonly infiltrating immune cells that included B cells, B memory cells, plasma cells, CD8+ and naïve CD4+ T cells, natural killer cells, macrophages, dendritic cells, and mast cells in COAD patients with low- and high-risk scores (Newman et al., 2015).

Tumor mutation burden (TMB) on the target gene was calculated by dividing the total number of mutations by the size of the coding region of the target gene. The GenVisR R package was used for the evaluation and analysis of the top 30 frequently mutated genes among the low- and high-risk groups (Skidmore et al., 2016). The copy number variation (CNV) in these groups of patients was calculated using chi-square test (Feng et al., 2020). The Rcircos R package was used to construct Circos plots to determine the CNV in chromosomes (Zhang et al., 2013). Utilization of the mRNA expression-based stemness index (mRNAsi) and epidermal growth factor receptor (EGFR)-mRNAsi data obtained from Tathiane Malta, University of São Paulo, was used to estimate the tumor stemness index to assess the dedifferentiation potential of oncogenic cells (Malta et al., 2018).

Based on the Genomics of Drug Sensitivity in Cancer (GDSC)⁵, the efficacy of distinct antineoplastic drugs on COAD TCGA samples was determined using the pRRophetic R package. The half-maximal inhibitory concentration (IC50) of a specific chemotherapeutic drug was obtained by ridge regression analysis. The prediction accuracy was measured by the 10-fold cross-validation of the GDSC cell line expression profile data. In addition, the Broad Institute Connectivity Map (CMap⁶) was used to predict the small candidate molecules based on the top 1,000 differentially expressed genes (Subramanian et al., 2017).

All statistical analyses were performed using R software (version 3.6.1). The Wilcoxon test was used for the comparison of two independent non-parametric samples. The Kruskal–Wallis test was used for multiple independent samples. The Kaplan–Meier survival curves were compared with the log-rank test. Univariate and multivariate Cox proportional hazard regressions were performed to identify independent prognostic factors associated with OS. A p-value < 0.05 was considered statistically significant.

The primary screening and analysis of TCGA samples by GSEA identified 226 DE-GRGs in COAD samples compared with normal tissues (Figure 2A). Of these, 154 were overexpressed and 72 were underexpressed (Figures 2B,C; Supplementary Table 2). The DE-GRGs were significantly augmented in BP terms related to carbohydrate metabolism, nucleotide metabolism, AMP-activated protein kinase, and hypoxia-inducible factor signaling pathway (Figure 2D and Supplementary Table 2). KEGG pathway analysis shows that DE-GRGs are involved in HIF-1 signaling pathway, glycolysis/gluconeogenesis, AMPK signaling pathway, and galactose metabolism (Figure 2E and Supplementary Table 3).

To investigate the role of GRGs in COAD, a PPI network with 177 nodes and 1,216 edges was constructed using Cytoscape String software (Figure 3A). In addition, a co-expression network with 69 nodes and 701 edges was created using the MCODE tool (Figure 3B) to identify the key modules of GRGs. To explore the prognostic value of these GRGs, univariate Cox regression analysis was performed. Seventeen prognostic-associated candidate hub GRGs were revealed (Figure 3C). These candidate hub GRGs were analyzed by multiple stepwise Cox regression to investigate their impact on patient survival and clinical outcomes. Eight hub GRGs were independent predictors in COAD patients (Figure 3D and Table 2).

A prediction model incorporating a prognosis-related signature of the eight hub GRGs was determined with the following calculation of the risk score:

Survival analysis in the low- and high-risk groups was determined using the Kaplan–Meier survival plot. The analysis revealed a relatively low OS rate in patients with high-risk score compared with the rate in the low-risk group of patients (Figure 4A). The risk scores, survival scores, and heatmap of GRG expression among the low- and high-risk patients are presented in Figure 4B. In addition, the area under the curve (AUC) analysis of COAD patients at 1, 3, and 5 years revealed respective OS rates of 73.4, 74.2, and 78.8% (Figure 4C). To verify the effectiveness of this model, a risk score was calculated in the GSE17536 cohort of patients, and they were classified into the high-risk group (n = 71) and low-risk group (n = 106) based on median risk score. Surprisingly, the results obtained from the GSE17536 cohort revealed that the high-risk score patients had lower OS rates than the patients with low-risk score, which was identical to the results obtained in TCGA cohort analysis (Figures 4D–F).

Considering the importance of glycolysis in the prognosis of colon cancer, we further analyzed the relationship between the eight GRGs and the clinical characteristics of colon cancer, including age, sex, tumor stage, and TMN stage. The heatmap (Figure 5A) indicated that the expression of CHPF2, ENO3, STC2, and P4HA1 was upregulated in the high-risk group, whereas that of PPP2CB, PGM2, PMM2, and PPARGC1A was downregulated in the high-risk group. Furthermore, a substantial difference in risk scores was evident for the different T, N, and M grades and tumor stages. The risk score increased with tumor progression (Figure 5B and Supplementary Figure 2A). To better predict the prognosis of the model, stratification analysis was used to confirm whether the model retained the ability to predict OS in different clinical subgroups. The OS rate of colon cancer patients with high-risk score, compared with patients in the low-risk group, was worse in the T3 + T4 and M0 groups (Figure 5C). Similarly, the risk model was good at predicting the OS rate of patients >65 or ≤65 years of age, male or female patients, stage I + II or III + IV, grade 1 or grade 2 + 3, and N0 or N1 (Figure 5C and Supplementary Figure 2B). The collective data indicate that the model may be a good prognostic index for patients with COAD.

Gene set enrichment analysis of the GEO and TCGA databases revealed the aberrant activation of multiple signaling pathways. These included oxidative phosphorylation, chemokine receptor interaction, glycerol metabolism, chemokine signal pathway, metabolism or exogenous substances through cytochrome P450, and transforming growth factor-beta signaling pathways in high-risk COAD patients (Figures 6A,B).

CIBERSORT along with the LM22 matrix was used to assess immune cell infiltration in the low- and high-risk groups of COAD patients. The GEO and TCGA and dataset analyses indicated that CD4+ and CD8+ T cells, macrophages, and mast cells were the predominant infiltrating immune cells in COAD (Figures 6C,D). Furthermore, infiltrations of resting CD4+ memory T cells, resting dendritic cells, and activated dendritic cells were significantly downregulated in the high-risk patients. In contrast, the distributions of M0 macrophages and eosinophils were opposite in the low- and high-risk TCGA groups (Figure 6E). Similarly, in the GSE17536 dataset, the infiltrating immune cells comprised significantly low proportions of CD4+ memory resting T cells and resting dendritic cells and high infiltration of M0 macrophages compared with those in the low-risk group (Figure 6F). However, there was no significant difference in the infiltrations of activated dendritic cells and eosinophils between the low- and high-risk groups.

To explore the underlying reasons for the difference in prognosis of patients in the low- and high-risk groups, we analyzed the TMB, CNV, and tumor stemness among the low- and high-risk groups. Waterfall plots were constructed to demonstrate the TMB and to differentiate TMB among the two risk groups. As shown in Figures 7A,B, most genes displayed higher mutation frequencies in the high-risk group. These genes included TP53, TTN, and KRAS. We calculated the TMB for each patient in both risk groups. The mean TMB was not significantly different between the two groups (p > 0.05; Figure 7C). In addition, distinct variations in CNV were evident in 920 genes predominantly located in chromosomes 1, 11, 14, 17, 18, and 19 (Figure 7D). A one-class logistic regression analysis was performed to calculate (mRNAsi) and epigenetic regulation-based stemness index (EREG-mRNAsi). The high-risk group had higher stemness indices compared with the low-risk group (Figures 7E,F). To evaluate the therapeutic effect of the risk model, the pRRophetic R package was used to analyze data downloaded from the GDSC database for the prediction of the sensitivity of distinct chemotherapeutic drugs on GRGs. Estimated IC50 values indicated that the low-risk patients were more sensitive to the anticancer drugs olaparib, veliparib, axitinib, metformin, and rapamycin (p < 0.05; Figure 7G).

To further judge whether the GRG model can be used as an independent factor to predict the prognosis of patients with colon cancer, Cox regression analysis was performed to reveal an association between risk score and OS of the patients.

Univariate Cox regression analysis revealed the significant association of risk score with OS in both TCGA cohort (hazard ratio [HR] = 1.382, 95% confidence interval [CI] = 1.284–1.488, p < 0.001) and GSE17536 cohort (HR = 1.252, 95% CI = 1.125–1.393, p < 0.001); Figures 8A,B). Multivariate Cox regression analysis with correction for distinct confounding parameters still revealed the risk score as an independent predictor for OS (TCGA cohort: HR = 1.382, 95% CI = 1.261–1.515, p < 0.001; GSE 17536 cohort: HR = 1.246, 95% CI = 1.104–1.406, p < 0.001; Figures 8C,D). Subsequently, ROC analysis was used to explore whether the risk score model was more accurate in predicting OS than clinical characters. As shown in Figures 8E,F, the predictive ability of the risk score model was stronger than that of other clinical characteristics. To establish a quantitative tool capable of predicting the clinical application of OS in patients with colon cancer, we established a nomogram integrating the risk score and clinical characters to calculate 1-, 3-, and 5-year OS rates of patients (Figures 8G,H). Calibration plots indicated that the nomogram versus an ideal model showed high consistency in TCGA and GSE17536 cohorts (Figures 8I,J).

Sixty potential small molecule drugs targeting genes for GRGs were identified by CMap database analysis (Figure 9). The mechanisms of the drugs were distinct. Among the 60 drugs, three (bicuculline, NCS-382, and pentylenetetrazol) are gamma aminobutyric acid (GABA) receptor antagonists and three other drugs (prednisolone, diflorasone, and fluocinonide) are glucocorticoid receptor agonists. The findings indicated that drugs targeting the signature GRGs might have therapeutic importance in the treatment of colon cancer.

qRT-PCR was performed from RNA extracted from 43 pairs of colon cancer and adjacent normal tissues to observe the relative expression levels of eight GRGs (P4HA1, STC2, CHPF2, PMM2, PGM2, ENO3, PPARGC1A, and PPP2CB). P4HA1, STC2, CHPF2, PMM2, PGM2, and ENO3 were expressed at higher levels in COAD tissues, with lower expression of PPP2CB and PPARGC1A (Figure 10A). The findings were consistent with the results of analyses of the Gene Expression Profiling Interactive Analysis (GEPIA) and UALCAN online databases (Supplementary Figure 3).

Similarly, immunohistochemistry (IHC) staining of colon cancer tissue with antibodies to P4HA1, STC2, CHPF2, PMM2, PGM2, ENO3, PPARGC1A, and PPP2CB had higher expression of the P4HA1, STC2, PMM2, CHPF2, and ENO3 GRGs in COAD tissue and lower expression of the PPP2CB, PPARGC1A, and PGM2 GRGs in the same tissues (Figures 10B,C). In addition, there were significant differences in the expression of STC2 in the COAD tissues at different stages of COAD (I, II, III, and IV; Figure 10D). These findings suggest that the expression level of STC2 changes with the progression of tumors.

To further explore the potential prognostic value of eight GRGs in COAD patients, we performed the Kaplan–Meier survival analysis to determine the effect of the eight GRGs on OS in TCGA COAD patients. Overexpression of STC2, ENO3, P4HA1, and CHPF2 and low expression of PPP2CB, PPARGC1A, PMM2, and PGM2 are associated with lower OS rate in COAD patients (Supplementary Figure 4).

Considering the potential role of STC2 in CRC progression, we performed in vitro proliferation, colony formation, migration, and invasion assays. The high expression of STC2 was associated with poor outcomes of CRC patients. We further evaluated the effect of STC2 depletion on the cellular process and BP in the CRC cells. Functionally, knockdown of STC2 by siRNA resulted in decreased cellular proliferation and colony formation by HCT116 and SW480 cells (Figures 11A,B). Similarly, knockdown of STC2 in HCT116 and SW480 cells using short hairpin RNA decreased cellular migration and invasion compared with control cells (Figures 11C,D). The findings confirmed the potential role of STC2 in the progression or CRC.