Microdissecting the Hypoxia Landscape in Colon Cancer Reveals Three Distinct Subtypes and Their Potential Mechanism to Facilitate the Development of Cancer

Background Hypoxia contributes to tumor progression and confers drug resistance. We attempted to microdissect the hypoxia landscape in colon cancer (CC) and explore its correlation with immunotherapy response. Materials and Methods The hypoxia landscape in CC patients was microdissected through unsupervised clustering. The “xCell” algorithms were applied to decipher the tumor immune infiltration characteristics. A hypoxia-related index signature was developed via the LASSO (least absolute shrinkage and selection operator) Cox regression in The Cancer Genome Atlas (TCGA)-colon adenocarcinoma (COAD) cohort and validated in an independent dataset from the Gene Expression Omnibus (GEO) database. The tumor immune dysfunction and exclusion (TIDE) algorithm was utilized to evaluate the correlation between the hypoxia-related index (HRI) signature and immunotherapy response. Quantitative reverse transcription polymerase chain reaction (qRT-PCR) and western blotting were performed to verify the mRNA expression levels of five key genes. The Cell Counting Kit-8 (CCK-8) assay and flow cytometry were performed to examine the cell viability and cell apoptosis. Results Patients were classified into hypoxia-high, hypoxia-median, and hypoxia-low clusters in TCGA-COAD and verified in the GSE 17538 dataset. Compared with the hypoxia-low cluster, the hypoxia-high cluster consistently presented an unfavorable prognosis, higher immune scores, and stromal scores and elevated infiltration levels of several critical immune and stromal cells. Otherwise, we also found 600 hypoxia-related differentially expressed genes (HRDEGs) between the hypoxia-high cluster and the hypoxia-low cluster. Based on the 600 HRDEGs, we constructed the HRI signature which consists of 11 genes and shows a good prognostic value in both TCGA-COAD and GSE 17538 (AUC of 6-year survival prediction >0.75). Patients with low HRI scores were consistently predicted to be more responsive to immunotherapy. Of the 11 HRI signature genes, RGS16, SNAI1, CDR2L, FRMD5, and FSTL3 were differently expressed between tumors and adjacent tissues. Low expression of SNAI1, CDR2L, FRMD5, and FSTL3 could induce cell viability and promote tumor cell apoptosis. Conclusion In our study, we discovered three hypoxia clusters which correlate with the clinical outcome and the tumor immune microenvironment in CC. Based on the hypoxia cluster and HRDEGs, we constructed a reliable HRI signature that could accurately predict the prognosis and immunotherapeutic responsiveness in CC patients and discovered four key genes that could affect tumor cell viability and apoptosis.


Introduction
Colon cancer (CC) is the ffth most frequent malignant disease with 1,148,515 new cases diagnosed in 2020 and accounting for 576,858 cancer-associated deaths around the world [1]. Te 5-year survival probability for colorectal cancer ranges from 90% in early-stage patients to 14% in distant-stage patients [2]. Te American Joint Committee on Cancer (AJCC) staging is a critical assessment system for the treatment management of CC [3], and patients with stage III or high-risk stage II may need to undergo a combination treatment of curative resection and adjuvant therapy [4]. However, most of the distant-stage patients miss the radical surgical opportunity and die due to metastasis or recurrence.
Owing to the tumor heterogeneity and diverse molecular pattern, patients with the same AJCC stage exhibit tremendous survival diferences. Tus, it is imperative to conduct in-depth microdissection and develop new prognostic biomarkers for patients with CC.
Hypoxia is a specifc hallmark of solid tumors, owing to the unrestricted growth and abnormal vascularization during the tumor progression [5]. Hypoxia promotes the tumor metastatic cascade, including invasion, migration, and distant metastasis [6]. Te hypoxia-inducible factor (HIF)-1α pathway contributes greatly to the modulation of hypoxia-related downstream gene expression and pathway activity in cancer cells under hypoxic conditions [7]. Hypoxia also promotes the epithelial-to-mesenchymal transition (EMT) process and facilitates the invasion of CC cells by activation of HIF1A, whereas treatment with HIF1A-specifc small interfering RNAs (siRNAs) suppresses these processes [8]. Our previous study [9] has constructed a hypoxia-related long noncoding RNAs signature that is tightly associated with the prognosis and drug sensitivity in patients with hepatocellular carcinoma. Te HIF-1α signaling pathway also confers drug resistance under hypoxic stress in colorectal carcinoma [10,11]. Hence, we speculate that the hypoxic exposure level in tumor tissues probably has a critical impact on the prognosis and treatment efectiveness of CC.
Over the past decade, immunotherapeutic treatment based on immune checkpoint inhibitors (ICIs) has resulted in revolutionary long-term benefts in the therapy of several cancer types [12]. ICIs such as anti-PD-1 (programmed cell death 1) and anti-PD-L1 (programmed cell death 1 ligand 1) have achieved a durable response in a subset of microsatellite instability-high (MSI-H) patients [12], whereas the MSI-L/MSS (MSI-low/microsatellite stability) patients who constitute the majority of CC patients have not obtained satisfactory benefts from ICI treatment. Interestingly, hypoxia has been reported to afect tumor plasticity, heterogeneity, and the immune resistance phenotype [13]. Hypoxia not only recruits myeloid-derived suppressive cells (MDSCs), cancer-associated fbroblasts (CAFs), and regulatory T cells (Tregs) to induce tumor immunosuppression [14] but also augments the expression level of immune checkpoints such as PD-L1 to promote tumor immune evasion [15]. Hence, targeting the hypoxic microenvironment may improve the efcacy of cancer immunotherapy [16]. Nevertheless, there is still a defciency in comprehensive delineation of the interplay among hypoxia, tumor immune infltrating patterns, and immunotherapy response in patients with CC.
In the current study, we discovered the hypoxia cluster in CC patients using unsupervised clustering based on two publicly available datasets (TCGA-COAD and GSE17538) and investigated the intrinsic correlation between hypoxia and the tumor immune microenvironment by the xCell algorithm and TIDE. Additionally, we developed a reliable hypoxia-related index (HRI) prognostic signature that exhibited good performance in predicting clinical prognosis and immunotherapy response in two independent datasets by the LASSO cox regression model. Finally, in vitro experiments were supplied to explore the results at the cell level. Our fndings may deepen the understanding of the hypoxia role in the tumor microenvironment and provide benefcial information for immunotherapy in CC.
Another publicly available, independent microarray dataset, GSE17538, was downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih. gov/geo/). Te TCGA-COAD cohort consisted of 402 primary CC samples and 39 adjacent normal tissues. Only 348 patients with complete clinical data and overall survival (OS) time of ≥1 month and 39 normal samples were used as the discovery cohort. Te GSE17538 dataset was composed of two subsets, GSE17536 (177 CC patients) and GSE17537 (55 CC patients), and the nonbiological batch was corrected using the "ComBat" function via the R "sva" package. In total, 210 patients with CC in GSE17538 with complete clinical and histopathological grade information were enrolled as the independent validation cohort. Detailed information on all enrolled patients in the previous two datasets is listed in Supplemental Table 1  (Table S1).

Microdissecting the Hypoxia-Specifc Cluster of CC.
Te "HALLMARK_HYPOXIA" gene set ("h.all.v7.2.symbols.gmt") includes 200 hypoxia-specifc genes (Table S2), which have been demonstrated to typically represent the biological process under hypoxia conditions and was gathered from the molecular signatures database (MsigDB) [18]. TCGA-COAD cohort (348 patients) and GSE17538 dataset (210 patients) were assigned into diferent groups by the unsupervised clustering method according to the expression of the previous 200 hypoxia-specifc genes, respectively, via the "km" method in the R "ConsensusClusterPlus" package. Survival analysis for hypoxia-specifc clusters was performed by the R "survival" package, and the survival diference among these clusters was determined by the log-rank test.

Identifying Hypoxia-Related Diferentially Expressed
Genes (HRDEGs). Analysis of diferentially expressed genes (DEGs) between the 348 COAD cancer samples and 39 adjacent normal samples in TCGA-COAD cohort was carried out by the "limma" package, according to the standard of the absolute value of log2 (fold change) greater than 1 and an adjusted p value less than 0.05. With the same method and criteria, DEGs between the hypoxia-high and hypoxia-low clusters were further examined. HRDEGs were defned as the intersection of the previous two gene lists.

Development of the HRI Signature.
A univariable Cox regression model was applied to screen the prognostic HRDEGs in TCGA-COAD cohort. Te LASSO penalty Cox regression model, which can avoid overftting and select the most contributive variables through tuning the penalty parameter, was employed to develop the optimal HRI signature using the "glmnet" package [20]. Te fnal HRI score formula is defned as follows: risk score � n k�1 expk * coefk, where k is the sequence number of the prognostic gene in the HRI signature, expk represents the corresponding gene expression of each patient, and coefk represents the corresponding LASSO coefcient.
2.6. Evaluation and Validation of the Prognostic Capability of the HRI Signature. HRI scores of CC patients in TCGA--COAD cohort (discovery dataset) and GSE17538 (validation dataset) were calculated using the previous formula. Patients in each dataset were assigned to the HRI high-or low-risk group according to their respective median HRI scores. Survival analysis for each dataset was carried out by the "survival" package, and the survival diferences were determined by the log-rank test. Time-dependent receiver operating characteristic (ROC) curves were drawn to evaluate the performance for prognosis prediction using the "timeROC" package. Multivariable Cox regression was conducted to determine whether the HRI signature was independent of other clinical parameters (age, sex, AJCC stage, and histopathological grade) in prognostic prediction.

Single Sample Gene Set Enrichment Analysis (GSEA).
Te gene list of critical immune function pathways (Table S4) was collected from the previous studies [21]. Single-sample GSEA (ssGSEA) [22], a particular kind of GSEA that can calculate the relative score for a predefned gene list at a single sample level, was utilized to calculate the relative scores of the previous immune function pathways using the "GSVA" package in R.

Analyzing the Immune Landscape of Hypoxia-Specifc
Clusters. Te "xCell" algorithm, which can efectively infer immune and stromal cell abundance from the mixture transcriptomic profles [23], was applied to comprehensively delineate the tumor immune microenvironment (TIME).

Evaluating the HRI Predictive Ability of Immunotherapy
Response. Te tumor immune dysfunction and exclusion (TIDE) algorithm, which can calculate the TIDE scores representing the dysregulation of tumor immune escape for tumor samples and function as a representative biomarker to predict responsiveness to immune checkpoint blockade [24], was employed to examine the HRI predictive capability of immunotherapy response in CC patients.

Quantitative Reverse Transcription Polymerase Chain
Reaction (qRT-PCR). Tirty pairs of clinical samples (including tumors and corresponding adjacent normal samples) of patients diagnosed with CC were gathered at Nanfang Hospital of Southern Medical University. Te samples were immediately preserved at −80°C postcollection after surgical resection until RNA extraction. All patients gave informed consent for sample collection and usage. Te present research was supported by the Institutional Ethical Committee Board of Nanfang Hospital (NFEC-201809-K3). Total RNA from 30 pairs of clinical tissues was isolated using an RNAex Pro Reagent (Accurate Biology, China). qRT-PCR reactions were performed using the Evo M-MLV RT Premix for qPCR (Accurate Biology, China) and SYBR ® Green Premix Pro Taq HS qPCR Kit (Accurate Biology, China). GAPDH was utilized as the internal standard, and each sample was analyzed in triplicate. All PCR primer sequences are presented in Table S5. Relative quantifcation of mRNA expression levels of RGS16, SNAI1, CDR2L, FRMD5, and FSTL3 was analyzed via the 2 −ΔΔCt method.
2.11. Cell Culture and Cell Transfection. Human colon cell line HCT116 was obtained from ATCC. Ten, the cells were cultured in DMEM with 10% FBS at 37°C in 5% CO 2 .
Te plasmid and scramble were purchased from Biosystems (General Biosystems, Anhui, China). siRNA and siRNA scramble were obtained from the GenePharma Corporation (Shanghai, China). According to the introduction, all siRNA and vectors were transfected using a lipofectamine 3000 transfection kit (Invitrogen, USA). qRT-PCR was performed to test the transfection efciency.

Western
Blot. Proteins were extracted using RIPA (CWBIO, China), subjected to SDS-PAGE gel electrophoresis, and then transferred to a nitrocellulose membrane, incubated with primary antibodies, and incubated overnight at 4°C. Te secondary antibody was then incubated for 1 h at room temperature. Immobilon ECL substrate was used for signal detection and image acquisition. transfected with siRNA or plasmid were placed on 96-well plates and cultured for 24 h, 48 h, 72 h, and 96 h. Ten, the OD450 value was detected using the Termo Scientifc Varioskan Flash spectrophotometer (Termo Scientifc, Finland).

Statistical Analysis.
Numerical variable diferences with normal distribution were determined using Student's ttest or analysis of variance for two or more groups, respectively. Te Wilcoxon rank-sum test or Kruskal-Wallis test were performed to determine the numerical variable diferences with nonnormal distribution for two or more groups, respectively. Categorical variable diferences were examined via the chi-square test. Spearman correlation analysis was conducted to investigate the correlation between the continual variables. Univariate Cox and LASSO penalty Cox regression analyses were utilized to perform survival analyses. Survival diferences were examined by the Kaplan-Meier curve and log-rank test. A two-tailed p value of <0.05 was set to indicate statistical signifcance. For multiple testing, the p value was corrected by the Benjamini-Hochberg method. We utilized R software (version 3.6.3) to perform all the statistical analyses.

Te Discovery of Hypoxia-Related Cluster
Using Unsupervised Clustering. In total, 348 patients with complete clinical information in TCGA-COAD cohort were categorized into three diferent clusters by unsupervised clustering (Figure 1(a) and Figures S1A-S1B). Clusters 1, 2, and 3 consisted of 82, 182, and 84 patients, respectively. Te detailed lists are shown in Table S6. Te principal component analysis confrmed a clear distinction among the three clusters ( Figure 1(c)). To clarify the relationship between the clusters and hypoxia, the HIF1A messenger RNA (mRNA) expression, which represents the mRNA level of the master regulator HIF-1α under hypoxic conditions, was compared among the three clusters. Notably, cluster 3 possessed the highest HIF1A mRNA level, while cluster 1 exhibited the lowest HIF1A mRNA level (Figure 1(d)). GSVA further showed that cluster 3 had the highest activity in the "HALLMARK_HYPOXIA" pathway, whereas cluster 1 displayed the lowest pathway activity (Figure 1(e)). Tese results demonstrated that the previous three clusters were strongly correlated with hypoxia exposure in CC tissues. Henceforth, we defned clusters 3, 2, and 1 as the hypoxiahigh, hypoxia-median, and hypoxia-low subtypes, respectively. Survival analysis revealed a signifcant OS difference among the three hypoxia-specifc clusters (global p value � 0.045, Figure 1(f )). Te hypoxia-high subtype had the poorest OS outcome compared with the hypoxia-low (p � 0.031) and hypoxia-median (p � 0.046) subtypes.
To further verify the hypoxic landscape in CC, the independent microarray dataset GSE17538 was explored using the same unsupervised clustering method. Notably, 210 patients with complete clinical characteristics in GSE17538 were likewise classifed into three diferent clusters (Figures 1(b), S1C-S1D, and 1(g)), with detailed lists shown in Table S7), namely, cluster 1 (77 patients), cluster 2 (68 patients), and cluster 3 (65 patients). Similarly, cluster 3 had the highest level of HIF1A mRNA expression and the activity of the "HALLMARK_HYPOXIA" pathway, while cluster 1 exhibited the lowest level for the previous two indices (Figures 1(h) and 1(i)). Tus, we also defned clusters 3, 2, and 1 in GSE17538 as the hypoxia-high, hypoxiamedian, and hypoxia-low subtypes, respectively. In addition, there was a signifcant OS diference among the three clusters (global p value � 1.35e − 04, Figure 1(j)). Te hypoxia-high cluster showed the poorest OS outcome compared with the hypoxia-low (p value 3.22e − 05) and hypoxia-median (p value 0.019) clusters. Te previous results confrmed that the hypoxia exposure landscape is closely correlated with the clinical outcomes in patients with CC.

Distinct Molecular Patterns among the Hypoxia-Specifc
Clusters. Owing to the close relationship between hypoxiaspecifc clusters and clinical outcomes, we continued to explore the underlying molecular mechanisms. GSVA results for the hallmark gene sets showed that the relative activities of several tumor aggression-associated pathways, including "EPITHELIAL_MESENCHYMAL_TRANSITI ON," "ANGIOGENESIS," "MYOGENESIS," "API-CAL_JUNCTION," "APICAL_SURFACE," "HYPOXIA," and "IL6_JAK_STAT3_SIGNALING," were elevated in the hypoxia-high group compared with those in the hypoxialow group in both TCGA-COAD and GSE17538 datasets (Figures 2(a) and 2(b)).

Identifcation of HRDEGs.
In total, 1756 DEGs (1748 upregulated and 8 downregulated genes, Figure 3(a)) between the hypoxia-high and hypoxia-low clusters (|log2FC| greater than 1 and adjusted p value less than 0.05) were identifed. Using the same criteria, we acquired 2745 DEGs (1442 upregulated and 1303 downregulated genes, Figure 3(b)) between the tumor tissues and adjacent normal samples. Furthermore, 600 overlapping genes for the previous two gene lists (Figure 3(c), detailed lists shown in Table S8) were categorized as HRDEGs. Gene Ontology (GO) function enrichment analysis demonstrated that these HRDEGs were predominantly enriched in several biological process (BP) terms, including "extracellular matrix organization," "positive regulation of cell adhesion," and "cellsubstrate adhesion" (Figure 3(d)). KEGG pathway analysis further showed a strong linkage between the HRDEGs and the following pathways: "cytokine-cytokine receptor interaction," "PI3K−Akt signaling pathway," and "focal adhesion" (Figure 3(e)). Tese enriched terms were closely  3.4. Development of HRI Signature. TCGA-COAD cohort with 348 patients was utilized as the discovery cohort to construct an HRI signature. Te univariate Cox regression model yielded 22 prognostic genes out of the aforementioned 600 HRDEGs (Figure 4(a)). LASSO penalty Cox regression selected the optimal HRI signature according to the "lambda. min" standard, which represents the lambda (tuning parameter) with minimal cross-validation error. Ultimately, 11 selected optimal prognostic HRDEGs were incorporated to develop the HRI score signature ( Figures S2A-S2B, detailed gene list shown in Table S9). Te fnal HRI calculation formula was as follows (detailed formula development is described in the "Methods" section): HRI = (−0.140) * CD177 expression + 0.045 * CP expression + 0.006 * RGS16 expression + 0.013 * PGM5 expression + 0.206 * SNAI1 expression + 0.010 * CALB2 expression + 0.041 * OSBPL1A expression + 0.043 * CDR2L expression + 0.012 * FRMD5 expression + 0.096 * FSTL3 expression + 0.069 * TUBB2B expression. Te HRI scores for CC patients in TCGA-COAD cohort were calculated using the previous HRI calculation formula (Table S10). Patients were assigned into the HRI high-or low-risk groups based on the median HRI score. Survival analysis uncovered that the HRI high-risk group exhibited a signifcantly poorer OS outcome than the low-risk group (p = 6.321e − 06, Figure 4(b)). Time-dependent ROC curves showed that the areas under the curve (AUCs) of 1-, 3-, 5-, and 6-year survival predictions were 0.682, 0.699, 0.768, and 0.753, respectively (Figure 4(c)), indicating good prognostic prediction. To further validate the reliability of the signature, the HRI scores for 210 patients in the validation cohort GSE17538 were calculated using the same formula (Table S11). Because the data type of dataset GSE17538 (microarray data) was diferent from that of the TCGA-COAD cohort (sequencing data), we classifed all patients in GSE17538 into HRI high-risk or low-risk groups according to the median HRI score of the dataset GSE17538. Similarly, the HRI high-risk group possessed a poorer OS prognosis than the low-risk counterpart (p = 7.956e − 06, Figure 4(d)). Te AUCs for 1-, 3-, 5-, and 6-year survival predictions were 0.647, 0.645, 0.716, and 0.754, respectively (Figure 4(e)). Tese results verifed the robustness and reliability of the HRI signature in diferent platform datasets.
Subsequently, we investigated the relationship between HRI scores and HIF1A mRNA expression. Notably, the HRI high-risk group consistently showed higher HIF1A expression than the low-risk counterpart in both TCGA--COAD (Figure 4(f )) and GSE17538 (Figure 4(g)), demonstrating that the HRI scores indeed refected the hypoxic exposure level in CC tissues.

Correlation between the HRI Signature and Clinical
Parameters. Owing to the remarkable impact of the HRI scores on the patient's clinical outcomes, we investigated the     Journal of Oncology correlation between HRI scores and several critical clinical parameters. Results for TCGA-COAD cohort indicated that patients with T3-4, M1, N2, stage III-IV, and "Vascular Invasion" possessed higher HRI scores than patients with T1-2, M0, N0, stage I-II, and "nonvascular invasion," respectively ( Figures 5(a)-5(e)). Furthermore, patients with stage III-IV, histopathological grade 3, and recurrence in GSE17538 had elevated HRI scores compared with patients with stage I-II, grade 1 or 2, and nonrecurrence, respectively ( Figures S3A-S3C). Additionally, patients with high HRI scores had poorer disease-free survival outcomes in TCGA-COAD (p = 9.12e − 05, Figure 5(f )) and poorer recurrence-free survival outcomes in GSE17538 (p = 0.006, Figure S3D) than patients with low HRI scores.
To identify the independent predictive ability of the HRI signature, a multivariable Cox regression model was further performed on the TCGA-COAD and GSE17538 datasets. Te results indicated that age, AJCC stage, and HRI scores were independent prognostic predictors after adjusting for other clinical parameters such as sex in TCGA-COAD cohort ( Figure 5(g)). Similarly, stage and HRI risk scores were consistently independent of age, sex, and histopathological grade in GSE17538 ( Figure S3E). Te previous evidence demonstrated that the HRI signature can act as an independent indicator of prognosis in CC.

Diferent Molecular Patterns, TIME, and Immunotherapy
Response between the High-and Low-Risk Groups. To further explore the underlying molecular mechanism, we investigated the diferent molecular patterns and TIME between the two HRI risk groups. GSEA results displayed that several critical hallmark pathways, including "API-CAL_JUNCTION," "APICAL_SURFACE," "ANGIOGEN-ESIS," "HYPOXIA," "EPITHELIAL_MESENCHYMAL_ TRANSITION," and "P53_PATHWAY," were substantially enriched in the high-risk group in both TCGA-COAD ( Figure S4A) and GSE17538 ( Figure S4C) datasets. Furthermore, KEGG pathways such as "ADHERENS_JUNCTION," "FOCAL_ADHESION," and "PATHWAYS_IN_CANCER" were signifcantly enriched in the group with high HRI scores in both TCGA-COAD ( Figure S4B) and GSE17538 ( Figure S4D) datasets. Tese results suggest that hypoxia contributes to tumor aggression through the abovementioned oncogenic pathways. Te "xCell" algorithm revealed that the high-risk group holds a higher abundance of macrophages, fbroblasts, and endothelial cells and higher stroma scores and microenvironment scores than the low-risk group in TCGA-COAD ( Figure 6(a)). Te high-risk group in GSE17538 possessed a higher infltrating level of macrophages and higher immune scores and microenvironment scores than the low-risk counterpart ( Figure S5A). Te ssGSEA results displayed that the high-risk group consistently possessed higher scores in several critical immune pathways such as "check−point" and "T_cell_co−inhibition" than the low-risk group in both TCGA-COAD ( Figure 6(b)) and GSE17538 ( Figure S5B) cohorts. Moreover, the mRNA expression level of PD-L1 (CD274) was signifcantly elevated in the HRI high-risk group compared with the low counterpart in both TCGA-COAD (Figure 6(c)) and GSE17538 ( Figure S5C) cohorts, suggesting distinct immune infltration characteristics between the two groups.In addition, compared with the hypoxia-low cluster, the hypoxia-high clusterconsistently presented higher immune scores, stromal scores, and elevatedinfltration levels of several critical immune and stromal cells (endothelialcells, fbroblasts, macrophages, dendritic cells, CD8+ T cells, CD4+ memory Tcells, B cells, and monocytes) in both TCGA-COAD and GSE17538 ( Figure S7A-B). Teabove evidence demonstrated that elevated hypoxia exposure levels in CC tissuescorrelated with higher stromal and immune cell infltration.
Using the TIDE algorithm, we estimated the TIDE scores for CC patients in the TCGA-COAD (Table S12) and GSE17538 (Table S13), respectively. Patients in the HRI high-risk group possessed higher TIDE scores than the corresponding low-risk patients in both TCGA-COAD ( Figure 6(d)) and GSE17538 ( Figure S5D). Moreover, HRI       scores consistently displayed a positive correlation with the TIDE scores in TCGA-COAD (Figure 6(e)) and GSE17538 ( Figure S5E), indicating that higher HRI scores represent greater immune evasion and immunotherapeutic resistance. Accordingly, the low-risk group was predicted to have a signifcantly higher ratio of immunotherapeutic responders than the high-risk group in both TCGA-COAD ( Figure 6(f )) and GSE17538 ( Figure S5F). Te previous results demonstrated that HRI scores representing hypoxia levels in CC tissues have the potential to predict the immunotherapy response.

Correlation Analysis of MSI Status with HRI Signature.
Te MSI status information for CC patients in the TCGA-COAD cohort was curated from TCIA (Te Cancer Immunome Atlas) database (https://tcia.at/home) [25].
Tere were 335 CC patients with complete MSI status in our TCGA-COAD dataset, including 46 MSI-H, 64 MSI-L, and 210 MSS (microsatellite stability), and 15 indeterminate cases, respectively. Te chi-square test revealed a statistically signifcant diference in the constitutive proportion of MSI status between the two HRI risk groups (p � 0.019, Figure S6A). Te HRI high-risk group presented an elevated ratio of MSI-H (17%) and MSI-L (25%) cases compared with the low-risk counterpart (11% and 15% of MSI-H and MSI-L cases, respectively). Subsequently, we stratifed CC patients into diferent subgroups according to their MSI status and performed a subgroup survival analysis. Notably, patients with high HRI scores consistently exhibited poorer OS prognosis than those in the HRI low-risk group, irrespective of MSI status (Figures S6B-S6D).

Validating the mRNA Expression of Five Key Genes by qRT-PCR.
Te HRI signature consisted of 11 HRDEGs, namely, CD177, CP, RGS16, PGM5, SNAI1, CALB2, OSBPL1A, CDR2L, FRMD5, FSTL3, and TUBB2B. Among the HRI prognostic signatures, CD177 was the only protective factor, and the other 10 genes were all risk factors for prognostic prediction in CC. Te mRNA expression levels of RGS16, SNAI1, CDR2L, FRMD5, and FSTL3 were higher in tumor samples than that in adjacent normal tissues in TCGA-COAD cohort (Figure 7(a)), suggesting that these fve key genes participate in the progression of CC. Tus, we experimentally investigated their mRNA expression levels in 30 pairs of clinical samples by qRT-PCR. Te results demonstrated that RGS16 (Figure 7(b)), SNAI1 (Figure 7(c)), CDR2L (Figure 7(d)), FRMD5 (Figure 7(e)), and FSTL3 (Figure 7(f )) consistently exhibited signifcantly higher relative mRNA expression levels in CC tumor samples than in paired adjacent normal tissues.

Te Validation Experiment In Vitro.
To further determine the infuence of the previously selected diferentially expressed genes (RGS16, SNAI1, CDR2L, FRMD5, and FSTL3) on cell proliferation and apoptosis, we interfered with the expression of fve diferential genes and detected the cell proliferation and apoptosis levels. First, we tested the transfection level of the disruptor or plasmid. Te results showed that the expression of mRNA (Figure 8(a)) and protein levels (Figure 8(b)) in the siRNA group was induced, while the pLenti group could signifcantly upregulate the expression of mRNA and protein levels of genes. Subsequently, we detected the cell activity by CCK-8 experiment. Te results showed that the expression of RGS16 had no signifcant efect on cell proliferation and apoptosis (Figures 9(a) and 10(a)), while the high expression of SNAI1, CDR2L, FRMD5, and FSTL3 could promote the proliferation of cancer cells and inhibit the apoptosis of cancer cells, but inhibiting their expression could inhibit the proliferation of cancer cells and promote the apoptosis of cancer cells (Figures 9(b)-9(e) and 10(b)-10(e)).

Discussion
CC ranks the ffth most frequent malignant disease worldwide, and advanced-stage cases are associated with high mortality [2]. Tus, it is urgent to identify novel prognostic predictors and targeted biomarkers. Hypoxia in the tumor microenvironment is a specifc hallmark of solid tumors [5] and contributes to the tumor metastatic cascade [6]. Several studies have constructed diferent hypoxia-related gene signatures for predicting the clinical outcomes of colorectal cancer [26][27][28]. However, these studies mainly aimed to establish a prognostic signature for CC patients and lacked comprehensive microdissection of the hypoxia landscape and its correlation with immunotherapy in CC. Compared with the previously published literature, we identifed three hypoxia-specifc clusters and developed a novel HRI prognostic signature. As far as we know, this is the frst comprehensive investigation of the correlation of the hypoxia landscape with metabolic reprogramming, TIME, and immunotherapeutic response prediction in CC.
Te hypoxic tumor microenvironment in solid tumors maintains a selective pressure for tumor cells to adapt to the hypoxia response and promotes their invasion, migration, and dissemination [6]. Moreover, the HIF-1α pathway plays a pivotal role in the modulation of hypoxia-related downstream gene expression and biological processes in cancer cells under hypoxic conditions [7]. Hence, we classifed the CC patients into three diferent clusters based on the expression levels of the 200 genes in the "HALLMAR-K_HYPOXIA" gene set and verifed the relationship between the clusters and the HIF1A mRNA expression level. Hypoxia stress can decrease the expression of DUSP2 and increase cancer stemness and tumor growth in CC cells [29]. Hypoxia may promote EMT, invasion, and migration of CC cells by activation of HIF1A, whereas treatment with HIF1A-specifc siRNAs suppresses these processes [8]. In agreement with the fndings of the abovementioned studies, the hypoxia-high cluster in our study possessed a higher HIF1A mRNA expression and elevated relative scores in tumor aggression-associated pathways including "EPI-THELIAL_MESENCHYMAL_TRANSITION" and "AN-GIOGENESIS." Accordingly, the hypoxia-high cluster had more CC patients with "Vascular_Invasion" and "Recurred/ Progressed." Cancer cells undergo metabolic reprogramming to reconcile themselves to hypoxic stress [30]. Te HIF-1α pathway contributes greatly to metabolism alteration via glycolysis stimulation and oxidative phosphorylation (OXPHOS) suppression under hypoxic conditions during tumor development [14]. Upon hypoxic conditions, pancreatic ductal adenocarcinoma cells exhibited elevated HIF1A and HIF2A expression levels, increased expression of carbonic anhydrase 9, and activated glycolysis [31]. Our study also showed similar results in both TCGA-COAD and GSE17538 datasets. Tis phenomenon may be because of the cancer cells' metabolic plasticity and metabolic heterogeneity, depending on the complex tumor microenvironment [32].
Hypoxia stress also impacts TIME by inducing an immune suppression or immune evasion phenotype [33]. Te local hypoxic microenvironment recruits immunosuppressive cells, such as MDSCs, tumor-associated macrophages, and CAFs, and upregulates immune checkpoint expression to induce antitumor resistance [15]. CAFs at the invasive front of tumor tissues boost tumor progression and metastasis in CC [34]. In our study, the hypoxia-high cluster consistently had a higher immune score, stromal score, and estimate score than the hypoxia-low cluster. Furthermore, macrophages and fbroblasts showed elevated infltrating levels in the hypoxia-high group, supporting a positive correlation between hypoxia and tumor immune dysfunction. Accordingly, we speculate that the poor clinical outcomes of the patients in the hypoxia-high group partly depend on the immune suppression or evasion mechanism.
To further examine the clinical applicability, we developed a reliable HRI prognostic signature that is strongly correlated with critical clinical characteristics (T, N, M, AJCC stage, and tumor histological grade). Hypoxia-treated  Figure 9: Te CCK-8 assay. Te efect of expression of RGS16 (a), SNAI1 (b), CDR2L (c), FRMD5 (d), and FSTL3 (e) on cell viability was examined using the CCK-8 assay. * * * p < 0.001; * * p < 0.01; * p < 0.05. 16 Journal of Oncology CC cells have been reported to strengthen the metastatic ability of normoxic cancer cells [35]. HIF-1α is a master regulator of the hypoxia-response process of tumor cells under hypoxic conditions [7]. Hypoxia can promote EMT, invasion, and migration of CC cells by the activation of HIF1A [8]. In our study, the high HRI score group consistently possessed a higher HIF1A expression level than the low-risk counterpart, indicating the efectiveness of the HRI score to refect hypoxia exposure in CC tumor tissues. Additionally, patients with M1, N2, T3-4, stage 3-4, and tumor grade 3 had higher HRI risk scores than those with M0, N0, T1-2, stage 1-2, and grade 1, respectively. Tis further demonstrated that the HRI signature was strongly related to tumor progression and metastasis in CC. Immunotherapy involving anti-PD-1/PD-L1 has resulted in revolutionary therapeutic benefts for various cancer types, and MSI-H status has proven to be an efective predictor of immunotherapeutic efcacy [36]. However, MSI-L/MSS patients, who represent most CC patients, have not acquired a satisfactory response from immunotherapy. High tumor mutational burden (TMB) in tumors is linked to favorable clinical outcomes; however, the TMB varies markedly among diferent cancer types, and there is a lack of a well-defned standard of high TMB [37]. Our study demonstrated that the HRI signature is positively correlated with the TIDE score, which represents the immune dysfunction and exclusion of tumor samples, and patients with a low HRI score are predicted to be more responsive to immunotherapy. Te previous evidence suggests that the HRI score has the potential to be a complementary measure to MSI-H status and TMB in the personalized management of immunotherapy. Hypoxia gene sets were reported to be enriched in nonresponding pre-anti-PD-1 tumor samples with melanoma [38]. In agreement with these studies, CC patients with a high HRI score representing severe hypoxia exposure are predicted to have a lower response to immunotherapy. As targeting the hypoxic microenvironment may ameliorate the efects of cancer immunotherapy [16], we speculate that these patients with high HRI scores may acquire greater efcacy of immunotherapy in combination with antihypoxia drugs.
Te HRI signature consists of 11 HRDEGs, and we focused on the 5 key genes (RGS16, SNAI1, CDR2L, FRMD5, and FSTL3), which exhibited elevated expression levels in tumor tissues and are prognostic risk factors for CC. RGS16 has already been reported to possess a higher expression level in colorectal cancer tissue than in the corresponding normal tissue and serves as an unfavorable prognostic marker [39]. Overexpression of SNAI1 (also known as SNAIL) is linked to increased stemness and decreased radiation sensitivity in CC cells [40]. A previously published study [41] reported the FRMD5 is a novel downstream gene targeted by the β-catenin/TCF7L2 complex in CC cells. CDR2L is widely present in ovarian cancer tissues and is abundantly expressed in testicular and prostate cancer tissues [42]. Knockdown of FSTL3 remarkably inhibited the aggression phenotype of lung cancer cells [43]. In the subsequent cell activity and apoptosis experiments, we found that the low expression of SNAI1, CDR2L, FRMD5, and FSTL3 could reduce the activity of cancer cells and increase the apoptosis rate of cancer cells. But RGS16 does not exhibit similar functions. According to previous literature reports, high expression of SNAI1 can promote the invasion ability of cancer cells [44], low expression of FRMD5 can weaken the metastatic ability of cancer cells [45,46], and low expression of FSTL3 also has similar functions [47,48]. Te reason why RGS16 has no similar function may be that its mechanism of afecting prognosis is diferent from other genes. According to previous reports, the population with low expression of RGS16 presents a better prognosis than the population with high expression [39]. Terefore, RGS16 may afect the prognosis of patients by regulating the activity of immune cells and has no direct impact on the activity and apoptosis rate of cancer cells [49]. Collectively, these fve key genes may act as oncogenic genes that contribute to the progression of CC, and their molecular mechanism is worth further studying to explore new therapeutic targets.
Nevertheless, there are still several limitations to our study. Te HRI signature was identifed in TCGA-COAD cohort and validated in another independent dataset, but these public datasets are mostly attributed to retrospective studies and may induce indispensable biases to some extent. Tus, prospective research will be required at a future date. Furthermore, although the HRI score is demonstrated to have a reliable predictive capability of immunotherapy response in CC by bioinformatical analysis, well-designed clinical trials should be performed to further prove its clinical efectiveness.

Conclusion
In conclusion, we discover three hypoxia clusters (hypoxia-H, hypoxia-L, and hypoxia-M) which correlate with the clinical outcome and the tumor immune microenvironment in CC. Furthermore, we found 600 HRDEs. Based on the 600 HRDEGs, we constructed a reliable HRI signature that could accurately predict the prognosis and immunotherapeutic responsiveness in CC patients. Finally, we discover fve key genes which are diferently expressed between tumors and adjacent tissues. Of them, four genes could afect tumor cell viability and apoptosis.

Ethical Approval
Tis research work was approved by the Institutional Ethical Committee Board of Nanfang Hospital (NFEC-201809-K3).

Conflicts of Interest
Te authors declare that they have no conficts of interest.