Identification of Hypoxia-Associated Signature in Colon Cancer to Assess Tumor Immune Microenvironment and Predict Prognosis Based on 14 Hypoxia-Associated Genes

Abstract

Purpose

Colon cancer is the main malignant tumor of the digestive tract. Hypoxia is highly related to the occurrence, progression and tumor immune microenvironment (TIME) of cancer. The aim of this study was to identify a hypoxia-associated signature with high accuracy for predicting the prognosis and TIME of colon cancer.

Methods

Download colon cancer data from the GEO and TCGA databases. A novel hypoxia risk model was identified to predict the prognosis of colon cancer patients. Subsequently, GSEA, TIME and mutation analysis were performed in the hypoxia high and low risk score groups. Finally, the signature gene ANKZF1 was selected for functional verification at the cellular level.

Results

A novel hypoxia risk model was identified. The risk score was significantly associated with poorer overall survival in colon cancer, and could be used as an independent prognostic factor for colon cancer. GSEA analysis found that the processes related to stimulate tumor proliferation and anti-apoptosis were significantly enriched in the hypoxia high risk score group. The expression of immunosuppressive cells and most immune checkpoints in the high risk score group was significantly higher than that in the low risk score group. In vitro cell experiments showed that knockdown the expression of ANKZF1 could inhibit the proliferation, migration and invasion of colon cancer cells.

Conclusion

Hypoxia plays an important role in evaluating the TIME and predicting the prognosis of colon cancer.

Keywords:

• colon cancer
• hypoxia
• prognosis
• risk score
• tumor immune microenvironment

Introduction

Colon cancer is the main malignant tumor of the digestive tract and one of the most malignant tumors in the world, with a high incidence.¹ Colon cancer has a high recurrence rate and mortality.^2,3 Although surgery, chemotherapy, radiotherapy and immunotherapy for colon cancer have improved, survival rates remain low.^4,5 Therefore, it is necessary to identify reliable prognostic models to predict clinical outcomes.

Hypoxia is a typical feature observed in the tumor immune microenvironment (TIME). It is caused by insufficient oxygen supply.⁶ It has previously been reported that many transcription factors in tumor cells are activated under hypoxic conditions, thereby inducing downstream signals to regulate tumor angiogenesis and cell proliferation.⁷ There are more and more studies on the TIME, in which immune cells play a vital role in the development of tumors.^8–10 Interestingly, increasing evidence suggests an association between tumor hypoxia characteristics and tumor immunosuppression and immune escape. Receptors affected by hypoxia play an important role in the recognition and lysis of a range of NK cell targets. Therefore, under hypoxia conditions, NK cells reduce their ability to kill.¹¹ Hypoxia can promote immunosuppressive cells or immunosuppressive cytokines, which in turn block immune effector cells.¹² Colitis-associated colon cancer is induced in wild-type mice, and the hypoxia and T cell immunity in colon cancer tissues were analyzed. It was found that hypoxia enhances immune suppression by inhibiting CD4⁺ effector T cell and promoting regulatory T cells activity.¹³ Tumor hypoxia is an important obstacle to effective tumor treatment. Therefore, studying the relationship between tumor hypoxia and immunity is helpful to the effect of immunotherapy.

Currently, methods of research tumor hypoxia are still very limited. However, the existence of public databases of gene data and clinical data makes it possible to identify the relationship between hypoxia and prognosis. In this study, we downloaded data related to colon cancer from the GEO and TCGA databases. Subsequently, hypoxia-associated signature was identified to predict the prognosis and TIME of patients with colon cancer. In addition, hypoxia-associated gene ANKZF1 was selected for in vitro cell experiments to detect the effects on the proliferation, migration and invasion of colon cancer cells.

Materials and Methods

Data Source

We searched gene expression data and complete clinical annotations of colon cancer in the GEO and TCGA databases (Table 1). For the data in TCGA, RNA sequencing data (FPKM value: fragments per kilobase of exon per million reads mapped value) and clinical information were downloaded from UCSC Xena (https://gdc.xenahubs.net). Then, convert the FPKM value into a transcript per kilobase million (TPM) values. For the GSE17536 dataset in GEO, we download the original “CEL” file. The RMA algorithm in “affy” package was used for background adjustment and quantile normalization. The annotation file was used for annotation. If multiple probes correspond to the same gene, take the average value. Then, the gene expression matrix file was obtained. In addition, samples with an overall survival (OS) of less than 30 days were excluded based on clinical information. In this study, the TCGA-colon cancer dataset was used as training set, and the GSE17536 dataset was used as test set. The list of hypoxia-associated genes was obtained from the hallmark gene sets of the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/index.jsp).¹⁴

Table 1 Basic Information of the Data Set in This Study

Accession Number	Platform	Number of Colon Cancer Patients	Survival Data
GSE17536	Affymetrix Human Genome U133 Plus 2.0 Array	175	OS/DFS/DSS
TCGA	Illumina RNAseq	433	OS

Abbreviations: OS, Overall survival; DFS, Disease Free Survival; DSS, Disease-Specific Survival.

Identification of Hypoxia-Associated Genes Signature

A total of 200 hypoxia-associated genes were in the Molecular Signatures Database. Based on expression data in the TCGA-colon cancer dataset, univariate Cox regression analysis was performed on hypoxia-associated genes to screen out genes associated with prognosis of colon cancer (P <0.05). Subsequently, these genes were screened using the Lasso-cox regression analysis in the glmnet package to reduce the number of genes for the risk score model. In Lasso-cox regression analysis, a lasso penalty was applied, taking into account both shrinkage and variable selection. Through 10-fold cross-validation, the optimal value of the lambda penalty parameter was determined. The purpose of 10-fold cross-validation is to prevent overfitting.^15,16 The risk score for colon cancer patients was calculated using the following formula:

(Note: n, number of prognostic genes; exp_i: expression value of gene i, β_i: regression coefficient of gene i.)

Patients were divided into high and low risk score groups according to the median risk score. The same formula was used to calculate the risk score in the GSE17536 dataset.

Survival Analysis

The OS of the high and low risk groups was compared by Kaplan-Meier analysis using the “survival” and “survminer” packages based on TCGA-colon cancer and GSE17536 datasets. Univariate Cox analysis was used to identify potential prognostic factors based on TCGA-colon cancer dataset. Multivariate Cox analysis was used to determine risk score as an independent risk factor for OS in colon cancer patients based on TCGA-colon cancer dataset. The time-dependent receiver operator characteristic (ROC) curve was generated by the “timeROC” package, which was used to verify the accuracy of the risk model for predicting patient survival based on TCGA-colon cancer and GSE17536 datasets.

Functional Enrichment Analysis

GSEA was performed through Java program (http://software.broadinstitute.org/gsea/index) based on TCGA-colon cancer dataset. The underlying mechanisms were researched in the “Molecular Signatures Database” of c2.cp.kegg.v7.3.symbols.¹⁷ False discovery rate (FDR) ≤0.05 were considered statistically significant.

Estimation of TIME Cell Infiltration

We used the ssGSEA algorithm to quantify the relative abundance of TIME cell infiltration in each colon cancer sample in TCGA. From Charoentong’s research, we obtained a set of genes that mark each TIME infiltrating immune cell type.^18,19 This gene set stored a variety of human immune cell subtypes, including activated CD8 T cell, macrophage, regulatory T cell and so on. Immune score, stromal score, tumor purity and ESTIMATE score of each patient were calculated using the ESTIMATE algorithm through the “ESTIMATE” software packages in R.

Expression of Cancer Immune Cycle Negative Regulatory Genes

The cancer immune cycle has become the intellectual framework of cancer immunotherapy research. It describes the process of the immune system to eradicate cancer: cancer cell antigen release, cancer cell antigen presentation, initiation and activation, transport of T cells to tumor, infiltration of T cells into tumor, recognition of cancer cells by T cells and killing of cancer cells.²⁰ In order to investigate the expression of genes that negatively regulating these processes in low and high risk score groups, the Tracking Tumor Immunophenotype (http://biocc.hrbmu.edu.cn/TIP/index.jsp) was used to download genes signatures.²¹

Immune Checkpoint and Somatic Mutation

To clarify the potential relationship between hypoxia in the TIME and clinical immunotherapy, six immune checkpoints were randomly selected to analyze their expression in the high and low risk score groups based on the TCGA-colon cancer and GSE17536 datasets. Six immune checkpoints include PDCD1 (PD1), CD274 (PDL1), PDCD1LG2 (PDL2), CTLA4, CD80 and CD86. In addition, in order to identify the mutation status in the high and low risk score groups, somatic mutation data of colon cancer patients were downloaded from TCGA for somatic mutation analysis.

Cell Experiment

Colon cancer cell HCT116 was selected to knockdown ANKZF1. Subsequently, the expression level of ANKZF1 in the knockdown cells was detected by real time-PCR to screen out the effective interference target. GAPDH was used as the internal reference gene. The relative gene expression level was calculated using the 2^−ΔΔCt method. The forward and reverse primer sequences of ANKZF1 were 5’-GCAGAAATCCGGCAATCGAC-3’ and 5’-GCCCTTAGAAGACGCACCAA-3’, respectively. The forward and reverse primer sequences of internal reference gene GAPDH were 5’-CTGGGCTACACTGAGCACC-3’ and 5’-AAGTGGTCGTTGAGGGCAATG-3’, respectively. Subsequently, the cells were divided into normal control (NC, normoxia), Si-ANKZF1 (normoxia), NC+1% O₂ (hypoxia) and Si-ANKZF1+1% O₂ (hypoxia) groups. Cell counting kit-8 (CCK-8) kit was used to detect the cell activity of each group at 0 h, 24 h and 48 h. The absorbance was measured at 450 nm using a microplate reader. In addition, flat panel clone proliferation experiments were also performed and the results were photographed using camera photography. The number of cell clones formed was statistically analyzed using Image J. Subsequently, the cell migration and invasion of each group were detected by cell wound scratch assay and transwell assay. Matrix metallopeptidase 2 (MMP2) and matrix metallopeptidase 9 (MMP9) play key roles in cancer immune escape and may also play a role in colon cancer progression.²² Low expression of epithelial cadherin (E-cadherin) affects the progression of colon cancer and is associated with poor prognosis.²³ Reducing the expression of neural cadherin (N-cadherin) may be involved in the regulation of colon cancer progression.²⁴ Hypoxia-inducible factor-alpha (HIF-α) is a transcription factor that regulates hypoxia response and plays an important regulatory role in colon cancer.²⁵ Western blot was used to detect the protein expression levels of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in each group.

Expression of CA-IX and PD-L1

Electronic expression analysis of CA-IX and PD-L1 was performed in the TCGA (colon cancer tissue samples:paracancer tissue samples=433:41) and GSE110224 (colon cancer tissue samples:paracancer tissue samples =17:17) datasets. The Wilcoxon test was used for statistical significance analysis. Subsequently, immunofluorescence was also performed in colon cancer tissue samples and paracancer tissue samples. Paraffin sections of tissue samples were routinely dewaxed, repaired, and washed twice times (5 min/time) with phosphate-buffered saline (PBS). Then, the blocking serum was added dropwise, and the blocking solution was wiped off with filter paper after 60 minutes at 37°C in a wet box. After adding the primary antibodies (CA-IX: Proteintech, 66243-1; PD-L1: Proteintech, 66248-1), the slides were incubated overnight at 4°C, and washed with PBS five times (3 min/time). The secondary antibody (goat anti-rabbit IgG (H&L):MDL, MD912526) was added dropwise and the sections were incubation at room temperature for 60 min, and washed with PBS three times (5 min/time). Wipe off the PBS outside the specimen with filter paper, and then add 4’,6-diamidino-2-phenylindole (DAPI) dropwise and incubate for 10 minutes in the dark. Finally, the slides were sealed with glycerol and immediately observed under fluorescence microscope. Blue coloring was defined as positive expression.

Statistical Analysis

All statistics were performed using R software (version 3.6.3, https://www.R-project.org). Wilcox.test was used to screen for statistically differentially expressed genes and infiltrating immune cells. When plotting Kaplan-Meier curves, log-rank was used to test the significant difference in OS between groups. In addition, univariate and multivariate Cox regression analyses were performed. P <0.05 was considered statistically significant.

Results

Identification of Hypoxia-Associated Genes Signature Prognostic Model

We performed univariate Cox regression analysis on 200 hypoxia-associated genes, and then screened out 19 genes that were significantly related to the prognosis of colon cancer (Figure 1A). Subsequently, Pearson correlation analysis was performed for these genes (Figure 1B). In addition, in order to more fully understand and demonstrate the associations between these hypoxia-associated genes and their relationship to OS, interaction network was mapped using Cytoscape (Figure 1C). Then, Lasso-cox regression analysis was performed on these 19 genes (Figure 1D and E). Ultimately, 14 hypoxia-associated genes (ALDOB, ANGPTL4, ANKZF1, CSRP2, DDIT3, ENO3, GLRX, GPC1, PGF, PGM2, PPARGC1A, SIAH2, STC2 and TKTL1) were selected to construct the risk score model. We used the coefficients calculated by the lasso regression analysis to calculate the risk score. The formula was as follows: Risk Score = (−0.043 ×ALDOB) + (0.200 × ANGPTL4) + (0.404 × ANKZF1) + (0.311 × CSRP2) + (0.092 × DDIT3) + (0.199 × ENO3) - (0.146 × GLRX) + (0.040 × GPC1) + (0.068 × PGF) - (0.441 × PGM2) - (0.328 × PPARGC1A) + (0.278 × SIAH2) + (0.152 × STC2) + (0.115 × TKTL1).

Figure 1 Hypoxia gene signature in colon cancer. (A) Forest map of 19 hypoxia-associated genes by univariate Cox analysis; (B) Pearson correlation analysis of 19 hypoxia-associated genes; (C) Interaction network of 19 hypoxia-associated genes mapped by Cytoscape. The Orange circle represents that hypoxia-associated genes were the risk factors of OS. The blue circle represents that hypoxia–associated genes were the favorable factors for OS. The blue line represents a negative correlation between two genes, and the depth of the color represents the strength of the correlation. The Orange line represents a positive correlation between two genes, and the depth of the color represents the strength of the correlation. (D) The lasso regression analysis; (E) The partial likelihood deviance for the lasso regression.

Abbreviation: OS, overall survival.

Prognostic Value of Hypoxia Risk Signature

Patients in TCGA-colon cancer and GSE17536 datasets were divided into high and low risk score groups according to the median risk score. The compareGroups package was used for clinical information statistics of patients in high and low risk score groups (Tables S1 and S2). In this study, we found that mortality was significantly higher in the high risk group than in the low risk group (Figure 2A and B). Subsequently, the expression of 14 hypoxia-associated genes in the high and low risk score groups was analyzed. Heat map results showed that 10 of the 14 hypoxia-associated genes were lowly expressed in the low risk score group, suggesting that high risk score patients were prone to forming hypoxic microenvironment (Figure 2C). The differential expression of 14 hypoxia-associated genes also suggested different disease progression in patients in the high and low risk score groups. In addition, Kaplan-Meier analysis was performed to evaluate the prognostic value of colon cancer hypoxia signature. The results showed that the low risk score group had a better survival advantage than the high risk score group (Figure 2D). This suggests that hypoxia high risk score was related to the poor OS in the TCGA-colon cancer dataset, which has been further validated by the GSE17536 dataset.

Figure 2 Distributions of risk score (A), OS status (B), heat map of signature genes (C) and survival of patients with high risk score and low risk score (D). The TCGA-colon cancer dataset on the left is training set and the GSE17536 dataset on the right is test set. P < 0.05 was considered statistically significant.

Abbreviations: TCGA, the cancer genome atlas; OS, overall survival.

Hypoxia Signature and Prognosis Assessment

In order to assess the accuracy of hypoxia signature in predicting 1-year, 3-year and 5-year survival rates, time-dependent ROC analysis was performed using data from TCGA-colon cancer and GSE17536 datasets. In the case of AUC >0.5, the greater the AUC value, the higher the accuracy. In the TCGA-colon cancer dataset, the area under curve (AUC) for 1-year, 3-year and 5-year was 0.710, 0.708 and 0.782, respectively (Figure 3A). In the GSE17536 dataset, the AUC for 1-year, 3-year and 5-year was 0.66, 0.7 and 0.81, respectively (Figure 3B). This indicates that the hypoxia signature has high prognostic prediction accuracy. Univariate Cox analysis showed that risk score was significantly correlated with OS of colon cancer patients (P <0.05) (Figure 3C). Multivariate Cox analysis showed that the risk score was associated with poor OS and could be used as an independent prognostic factor for colon cancer patients (P <0.05) (Figure 3D).

Figure 3 Prognostic value of the hypoxia gene signature in colon cancer. (A) ROC curves showing the predictive efficiency of the hypoxia gene signature on the 1-, 3-, and 5-years survival rate in TCGA-colon cancer dataset; (B) ROC curves showing the predictive efficiency of the hypoxia gene signature on the 1-, 3-, and 5-years survival rate in GSE17536 dataset; (C) Univariate Cox analysis of hypoxia gene signature to evaluate the independent prognostic value of OS in colon cancer patients; (D) Multivariate Cox analysis of hypoxia gene signature to evaluate the independent prognostic value of OS in colon cancer patients. **P <0.01; ***P <0.001.

Abbreviations: OS, overall survival; ROC, receiver operator characteristic, AUC, area under curve.

Gsea

In order to further understand the relevant signaling pathways activated in the high risk score group, GSEA was used to compare high and low risk score groups. The results showed that the activities of tumor proliferation and anti-apoptosis related processes were higher in the high risk score group. For example, cell cycle, PI3K-Akt, JAK-STAT, MAPK, mTOR, VEGF and other pathways (Figure 4). Therefore, it is speculated that these pathways may be involved in promoting disease progression in patients in the high risk score group.

Figure 4 Gene set enrichment analysis (GSEA) results in the high risk score group.

Immune Status of Low and High Risk Score Groups

In this study, we investigate the ability of hypoxia risk signature to assess the TIME. The ssGSEA algorithm was used to quantify the infiltration level of immune cells in colon cancer samples. The boxplot of ssGSEA results shows that compared to the low risk score group, the high risk score group had higher infiltration levels of immunosuppressive cells (for example, regulatory T cell, macrophage), while the immune infiltration levels of activated CD4 T cell and neutrophil were lower (Figure 5A). In addition, ESTIMATE algorithm results also confirmed that in patients with high risk score, the immune score, stromal score, and ESTIMATE score were significantly lower than those of patients with low risk score, while the tumor purity was significantly higher than that of patients with low risk score (Figure 5B–E).

Figure 5 Difference in immune cells (A), immune score (B), stromal score (C), ESTIMATE score (D) and tumor purity (E) between high and low risk score groups. *P <0.05; **P <0.01; ***P <0.001;****P <0.0001. P < 0.05 was considered statistically significant.

High Risk Score Indicates an Immunosuppressive Microenvironment

Tracking Tumor Immunophenotype was used to download genes signatures. In the high risk score group, most of the genes related to the negative regulation of the cancer immune cycle were up-regulated (Figure 6). The results suggest that the cancer immune cycle is low active in high risk score patients.

Figure 6 Expression of cancer-immune cycle negative regulatory genes in high and low risk score groups. (A) Heat map of cancer-immune cycle negative regulatory genes in high and low risk score groups in the TCGA-colon cancer dataset. The clustering is constructed using the full chain method together with the Euclidean distance. Each row and column represents a cancer-immune cycle negative regulatory gene and a sample, respectively. Orange indicates above the reference channel (high expression genes). Blue indicates below the reference channel (low expression genes). (B) Heat map of cancer-immune cycle negative regulatory genes in high and low risk score groups in the GSE17536 dataset. The clustering is constructed using the full chain method together with the Euclidean distance. Each row and column represents a cancer-immune cycle negative regulatory gene and a sample, respectively. Orange indicates above the reference channel (high expression genes). Blue indicates below the reference channel (low expression genes). (C) Box plot of cancer-immune cycle negative regulatory genes in the TCGA-colon cancer dataset; (D) Box plot of cancer-immune cycle negative regulatory genes in the GSE17536 dataset. *P <0.05; **P <0.01; ***P <0.001;****P <0.0001.

Abbreviation: ns, no statistical significance.

Immune Checkpoint and Somatic Mutation Analysis

The expression of six immune checkpoints in high and low risk score groups was analyzed based on TCGA-colon cancer (Figure 7A) and GSE17536 (Figure 7B) datasets. The expression of most immune checkpoints was significantly higher in high risk score group. These results suggest that patients with high risk score tend to suppress the TIME by up-regulating immune checkpoints. In addition, we also analyzed the genes with the top 20 mutation frequency between patients with high risk score and patients with low risk score based on TCGA-colon cancer dataset (Figure 7C and D). The results showed that the frequency of genetic mutations was high in both groups. However, the mutated genes in the two groups were different.

Figure 7 Expression of immune checkpoints and genes mutations between high and low risk score groups. (A) Expression of immune checkpoints in high and low risk score groups in the TCGA-colon cancer dataset; (B) Expression of immune checkpoints in high and low risk score groups in the GSE17536 dataset; (C) Waterfall chart of genes mutations in the low risk score group in the TCGA-colon cancer dataset; (D) Waterfall chart of gene mutations in the high risk score group in the GSE17536 dataset. *P <0.05; **P <0.01.

Abbreviation: ns, no statistical significance.

Knockdown ANKZF1 Inhibited the Proliferation, Migration and Invasion of Colon Cancer Cell HCT116

Colon cancer cell HCT116 was selected to knockdown ANKZF1. The expression level of ANKZF1 in the knockdown cells was detected by real time-PCR to screen out the effective interference targets. The results showed that the knockdown effect of target 1 was the best (Figure 8A), so we chose Si-RNA1 for subsequent experiments. CCK8 assay showed that the activity of colon cancer cells significantly increased under hypoxia condition, and the activity of cancer cells in Si-ANKZF1 and Si-ANKZF1+1% O₂ groups was inhibited after ANZF1 knockdown (Figure 8B). Flat panel clone experiment showed that the proliferation of cancer cells significantly increased under hypoxia condition, and the proliferation of cancer cells in Si-ANKZF1 and Si-ANKZF1+1% O₂ groups was inhibited after ANKZF1 knockdown (Figure 8C and D). Cell wound scratch assay showed that the migration of cancer cells significantly increased under hypoxia condition, and the migration of cancer cells in Si-ANKZF1 and Si-ANKZF1+1% O₂ groups was inhibited after ANKZF1 knockdown (Figure 9A and B). Transwell assay showed that the invasion of cancer cells significantly increased under hypoxia condition, and the invasion of cancer cells in Si-ANKZF1 and Si-ANKZF1+1% O₂ groups was inhibited after ANKZF1 knockdown (Figure 9C and D). The above results suggest that knockdown the expression of ANKZF1 can help inhibit the proliferation, migration and invasion of colon cancer cells, and also indicate that ANKZF1 may play an important molecular regulatory role in the progression of colon cancer. Subsequently, protein expression levels of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in each group were detected by Western blot. The results showed that knockdown ANKZF1 expression also affected the expression of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in cancer cells (Figure 10A and B). This further suggests that ANKZF1 may play an important role in colon cancer progression.

Figure 8 Activity and proliferation of colon cancer cells HCT116 after knockdown ANKZF1. (A) Real time-PCR was used to detect the expression of ANKZF1 to screen out effective interference targets in HCT116; (B) CCK8 was used to detect the activity of cancer cells in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups; (C) Flat panel clone experiment was used to detect the proliferation of cancer cells in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups; (D) Histogram of cell proliferation detected by flat panel clone experiment. *P <0.05; **P <0.01; ***P <0.001; ^###P <0.001; ****P <0.0001.

Abbreviation: NC, normal control.

Figure 9 Migration and invasion of colon cancer cells HCT116 after knockdown ANKZF1. (A) Cell wound scratch experiment was used to detect the migration of cancer cells in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups. Observation under 40x microscope. (B) Histogram of cell migration detected by cell wound scratch experiment; (C) Transwell was used to detect the invasion of cancer cells in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups. Observation under 200x microscope. (D) Histogram of cell migration detected by transwell. *P <0.05; ***P <0.001; ****P <0.0001.

Abbreviation: NC, normal control.

Figure 10 Western blot was used to detect the protein expression levels of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups. (A) The expression bands of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups were detected by Western blotting; (B) The relative expression levels of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α proteins in NC, Si-ANKZF1, NC+1% O₂ and Si-ANKZF1+1% O₂ groups. *P <0.05; **P <0.01; ***P <0.001; ****P <0.0001.

Abbreviations: NC, normal control; ns, no statistical significance.

Expression Analysis of CA-IX and PD-L1

It was previously reported that CA-IX is a reliable marker of hypoxia and is highly expressed in hypoxia microenvironment.²⁶ PD-L1 is widely regarded as a potential immune marker and therapeutic target for cancers. Therefore, we verified the expression of CA-IX and PD-L1 in colon cancer tissue samples and paracancer tissue samples. We analyze the mRNAs expression of CA-IX and PD-L1 in TCGA dataset. The results showed that CA-IX expression was significantly increased in colon cancer (Figure S1A) (P <0.05). In addition, the expression of PD-L1 in colon cancer tissues also had an increasing trend but was not significant (Figure S1B). We further verified the mRNAs expression of CA-IX and PD-L1 in the GSE110224 dataset. The results showed that CA-IX and PD-L1 were significantly increased in colon cancer tissues (Figure S1C and D) (P <0.05). Subsequently, the proteins expression of CA-IX and PD-L1 in colon cancer tissue samples and adjacent tissue samples was also verified by immunofluorescence. Compared with normal tissues adjacent to cancer, the number and degree of cells staining positive for CA-IX and PD-L1 proteins in colon cancer tissues had increase trended (Figure S1E and F).

Discussion

Colon cancer is a cancer that develops from the epithelial cells lining the colon or rectum of the gastrointestinal tract and is one of the most common cancers in the world.²⁷ With the improvement of colon cancer detection and treatment methods, morbidity and mortality have been reduced, but 30–50% of patients have recurrence and metastasis within five years after treatment.²⁸ Previous studies have found that immune cell infiltration,²⁹ autophagy-related genes,³⁰ tumour-stroma ratio³¹ and other characteristics can be used to predict the prognosis of patients with colon cancer. However, there are few studies on the relationship between hypoxia and the prognosis of colon cancer. Hypoxia is a typical feature of TIME.⁶ Recent studies have shown that hypoxia plays an important role in tumor prognosis.^32–34 A study has shown that hypoxic-associated gene signature can be used to predict the prognosis of colon cancer patients.³⁵ In order to deepen our understanding of predicting the prognosis of colon cancer using hypoxia risk score model, we also identified a novel hypoxic-associated gene signatures based on mRNA analysis data. The results showed that colon cancer patients in the hypoxia high risk score group had a poor prognosis. Different from the study of Yixin Xu et al, this study was a novel risk score model identified based on 14 hypoxia-associated genes. Subsequently, risk score model gene ANKZF1 was selected for in vitro cell experiments to detect the effects on the proliferation, migration and invasion of colon cancer cells.

In this study, the identified risk model consisted of 14 hypoxia-associated genes (ALDOB, ANGPTL4, ANKZF1, CSRP2, DDIT3, ENO3, GLRX, GPC1, PGF, PGM2, PPARGC1A, SIAH2, STC2 and TKTL1), which were significantly associated with colon cancer prognosis. Aldolase B (ALDOB) has been reported to alter the traction force of colorectal cancer to hypoxia, which has been linked to tumor metastasis.³⁶ The expression level of the angiogenic gene angiopoietin like 4 (ANGPTL4) in colorectal cancer was up-regulated under hypoxia.³⁷ High expression of ankyrin repeat and zinc finger peptidyl tRNA hydrolase 1 (ANKZF1) is associated with lower OS in colon cancer through involvement in angiogenesis and several cancer signaling pathways.³⁸ In this study, knockdown the expression of ANKZF1 in colon cancer cells ANKZF1 can inhibit the proliferation, migration and invasion of cells. Moreover, knockdown ANKZF1 expression also affected the expression of MMP2, MMP9, N-Cadherin, E-Cadherin and HIF-α in cancer cells. This further suggests that ANKZF1 may play an important role in colon cancer progression. Cysteine and glycine rich protein 2 (CSRP2), an aggressive actin binding protein, was up-regulated by hypoxia in breast cancer cell lines.³⁹ In addition, CSRP2 was abnormally expressed in colorectal cancer and was a potential therapeutic target.⁴⁰ Studies in colon cancer have shown that DNA damage inducible transcript 3 (DDIT3) was significantly induced in the hypoxia response.⁴¹ Low expression level of enolase 3 (ENO3) was associated with positive prognosis of colon cancer and can be used as a prognostic biomarker in patients with colon cancer.⁴² Glutaredoxin (GLRX) was abnormally expressed in colon cancer.⁴³ Glypican 1 (GPC1) expression was increased in colon cancer and was significantly related to prognosis.^44,45 Previous studies have found that hypoxia can activate the expression of placental growth factor (PGF).⁴⁶ Furthermore, PGF can induce the formation of blood vessels in colon cancer.⁴⁷ The down-regulation of phosphoglucomutase 2 (PGM2) may be related to poor prognosis of colorectal cancer, and may also be involved in the regulation of the process of colorectal cancer by regulating cell cycle, purine metabolism and spliceosome.⁴⁸ PPARG coactivator 1 alpha (PPARGC1A) was significantly down-regulated in colon cancer stroma and was a potential prognostic biomarker.⁴⁹ Siah E3 ubiquitin protein ligase 2 (SIAH2) is a ring finger ligase, which is an important regulator of hypoxia activation pathway.⁵⁰ In addition, SIAH2 expression plays an important role in the adaptation of colorectal cancer cells to hypoxia.⁵¹ Down-regulation of stanniocalcin 2 (STC2) significantly inhibited the growth and migration of colon cancer cells under hypoxic conditions.⁵² STC2 was also a potential tumor biomarker for colorectal cancer diagnosis and prognosis.⁵³ Transketolase like 1 (TKTL1) expression can be induced in colon cancer cells under hypoxia condition, and plays a vital role in cell proliferation.^54,55 In addition, strong TKTL1 protein expression is associated with poor prognosis in colon cancer.⁵⁶ Therefore, the hypoxia risk model of colon cancer constructed by these 14 hypoxia-associated genes has potential reliability.

Tumor hypoxia is a negative prognostic factor associated with oncogenic signal activation, immune escape and resistance to therapy.⁵⁷ Macrophages stimulate an anti-tumor immune response under hypoxia.⁵⁸ Hypoxia can enhance immunosuppression by inhibiting the function of CD4+ effector T cells and promoting Treg activity, thus affecting the immune detection of colon cancer.¹³ In our study, the proportion of immunosuppressive cells such as macrophages and Treg in the hypoxia high risk score group was significantly higher, while the proportion of immune cells such as activated CD4 T cells and B cells was reduced. In addition, in patients with high risk score of hypoxia, the immune score, stromal score, and ESTIMATE score were significantly lower than those of patients with low risk score, while the tumor purity was significantly higher than that of patients with low risk score. These results imply that the hypoxia risk model may predict TIME of colon cancer. Up-regulation of immune checkpoints and their ligands in the TIME can modulate tumor immune escape.⁵⁹ A study found that immune checkpoint inhibitors can improve survival in patients with clear cell renal cell carcinoma.⁶⁰ In this study, immune checkpoints like PD1, PD-L1, CD80, and PDL2 were significantly up-regulated in the high risk score group. Therefore, identifying the underlying mechanisms of hypoxia in immune escape and immune checkpoints may help identify potential therapeutic targets for colon cancer. It was previously reported that CA-IX is a reliable marker of hypoxia and is highly expressed in hypoxia microenvironment.²⁶ CA-IX is overexpressed in colorectal cancer, and this overexpression leads to a poor prognosis.⁶¹ The expression of CA-IX is associated with the proliferation of colon cancer cells.⁶² In this study, CA-IX expression was up-regulated trend in colon cancer tissues. This further suggests that hypoxia plays an important role in mediating the progression of colon cancer.

Nevertheless, this study also has some limitations. The hypoxia risk score model obtained in this study based on public data lacks clinical validation to demonstrate its effectiveness. Although the risk score model was analyzed in two different databases, there are still potential risks. Therefore, it is necessary to collect samples for further research. The immunofluorescence results showed that compared with the control group, CA-IX was increase in colon cancer, but it was not significant. It may be caused by sample heterogeneity and technical issues with experimental operations. We will continue to study further in the later period.

Conclusion

In this study, a novel hypoxia risk model for colon cancer was identified based on 14 hypoxia-associated genes. It can not only predict the survival of colon cancer patients, but also evaluate TIME. Analysis of the relationship between hypoxia and prognosis and immune cells will help improve the level of cancer treatment and prognosis of colon cancer patients. In addition, in vitro cell experiments showed that knockdown the expression of the risk score model gene ANKZF1 in colon cancer cells can inhibit cell proliferation, migration and invasion. In short, this study lays the foundation for understanding the value of hypoxia-associated genes in colon cancer.

Data Sharing Statement

All data generated or analyzed during this study are included in this published article. We searched for colon cancer public data from TCGA (https://tcga-data.nci.nih.gov/tcga/) and GEO (http://www.ncbi.nlm.nih.gov/geo). The accession numbers are TCGA-colon cancer; GEO: GSE17536 and GEO: GSE110224, respectively.

Ethics Approval and Consent to Participate

This study was approved by the Ethics Committee of The First Affiliated Hospital of Hebei Medical University (2013106). All participants were informed as to the purpose of this study, and that this study complied with the Declaration of Helsinki. The informed consent was obtained from the all participants.

Disclosure

The authors declare that they have no competing interests.

Additional information

Funding

This study was supported by the Key Science and Technology Research Project of Health and Family Planning Commission of Hebei Province (grant no. 20221165), Natural science foundation of Hebei Province (No. H2020206374 & No. H2021206306), Clinical medicine excellent talents project of Hebei Province (No. LS202001).