Genome-Wide Analysis Reveals Hypoxic Microenvironment Is Associated With Immunosuppression in Poor Survival of Stage II/III Colorectal Cancer Patients

Background: Hypoxia is associated with a poorer clinical outcome and resistance to chemotherapy in solid tumors; identifying hypoxic-related colorectal cancer (CRC) and revealing its mechanism are important. The aim of this study was to assess hypoxia signature for predicting prognosis and analyze relevant mechanism. Methods: Patients without chemotherapy were selected for the identification of hypoxia-related genes (HRGs). A total of six independent datasets that included 1,877 CRC patients were divided into a training cohort and two validation cohorts. Functional annotation and analysis were performed to reveal relevant mechanism. Results: A 12-gene signature was derived, which was prognostic for stage II/III CRC patients in two validation cohorts [TCGA, n = 509, hazard ratio (HR) = 2.14, 95% confidence interval (CI) = 1.18 – 3.89, P = 0.01; metavalidation, n = 590, HR = 2.46, 95% CI = 1.59 – 3.81, P < 0.001]. High hypoxic risk was correlated with worse prognosis in CRC patients without adjuvant chemotherapy (HR = 5.1, 95% CI = 2.51 – 10.35, P < 0.001). After integration with clinical characteristics, hypoxia-related gene signature (HRGS) remained as an independent prognostic factor in multivariate analysis. Furthermore, enrichment analysis found that antitumor immune response was suppressed in the high hypoxic group. Conclusions: HRGS is a promising system for estimating disease-free survival of stage II/III CRC patients. Hypoxia tumor microenvironment may be via inhibiting immune response to promote chemoresistance in stage II/III CRC patients.


INTRODUCTION
Colorectal cancer (CRC) is the third most common cancer type in men (10.9% of all cancers) and the second in women (9.5% of all cancers) worldwide (1). Currently, the therapeutic treatment for individual patients is mainly based on the TNM staging (2). However, for stages II and III CRC patients, the current staging method cannot provide indications of how and when to use chemotherapy, neither is the possible effectiveness of chemotherapy. As a result, chemotherapy frequently leads to diverse unpredictable responses due to risk factors such as T stage, tumor differentiation, and microsatellite instability (MSI) status (3)(4)(5). A previous study has shown that some stage II CRC patients have worse prognosis, whereas some stage III patients are overtreated (6). Therefore, investigating new reliable biomarkers can help to choose the appropriate therapy at the right time for CRC patients.
Hypoxia is a common microenvironmental condition found in most solid tumors and is associated with poor prognosis (7)(8)(9). Hypoxia-inducible factors can be activated under the hypoxia condition in the cellular environment. Evidence has shown that these factors can promote the tumor invasion and metastasis (10)(11)(12). Furthermore, recent studies have found that the hypoxic tumor microenvironment could affect the chemotherapy efficiency to CRC patients (13,14). Therefore, it is important to identify prognostic hypoxia-inducible factors that could help to choose the appropriate therapy against the CRC.
Currently, hypoxic gene signatures are found to have prognostic and predictive effects for diverse cancers including head and neck, bladder, soft tissue sarcoma, and cervical cancers, which allow clinicians to coordinate the therapeutic agent use for patients with most benefit (7,8,(15)(16)(17)(18). In this study, we analyzed hypoxia-related genes from large amounts of CRC transcriptional data and created a hypoxia-related gene signature (HRGS) for CRC prognosis. The prognostic prediction value of the HRGS was systematically validated in a metavalidation cohort. In addition, we demonstrated a relevant mechanism underlying the poor survival outcome of the high hypoxic risk CRC patients at stage II/III, and this finding may help to detect out new therapeutic target for CRC patients. This would help improve therapeutic strategy for CRC patients.

Patients (Data Source)
This meta-analysis study used the gene expression data of frozen CRC tumor tissue samples from six public cohorts. To be included in the study, patients needed to meet the following inclusion criterion: patients with pathologic diagnosed as CRC. The exclusion criterion was patients without survival information. In total, data of 1,877 patients, including 309 CRC patients without adjuvant chemotherapy, were used. The two largest individual datasets, CIT/GSE39582 and The Cancer Genome Atlas (TCGA) CRC cohort, were used for training and independent validation. The remaining four microarray datasets (GSE14333, GSE17536, GSE37892, and GSE33113) obtained from the Gene Expression Omnibus (GEO) database were merged as a metavalidation cohort. Data of the TCGA CRC cohort were downloaded from Broad GDAC Firehose (http://gdac.broadinstitute.org/), and transcripts per million of level 3 RNA-Seq data in log2 scale were used for the analysis. Other datasets were obtained in processed format from GEO database through R using the Bioconductor package "GEOquery." The batch effects were corrected using "combat" algorithm implemented in the R package "sva, " and z scores of each gene were used for the following analyses. The staging classification in each dataset was based on pathologic stage. Data were collected from September 27 to December 26, 2018.

Construction and Validation of HRGS
The Molecular Signatures Database (MSigDB) is one of the most widely used and comprehensive databases of gene sets for performing gene set enrichment analysis (GSEA) (19)(20)(21). We created a list of hypoxia-related genes (HRGs) including all gene sets found in MSigDB (version 6.2) with the keyword "hypoxia." In order to construct a prognostic HRGS, we assessed the association between all HRGs found in this meta-analysis and the patients' disease-free survival (DFS) in GSE39582 dataset using the log-rank test with 1,000 times randomization (80% proportion of samples each time). HRGs, which have frequently been significantly associated with the patients' DFS, were selected as candidates for the construction of HRGS. To minimize the risk of overfitting, we applied a Cox proportional hazards regression model combined with the least absolute shrinkage and selection operator (LASSO) (glmnet, version 2.0-16). The penalty parameter was estimated by 10-fold cross-validation in the training dataset at 1 SE beyond the minimum partial likelihood deviance.
In order to separate patients into low-or high-risk groups, the optimal HRGS cutoff was determined by a time-dependent receiver operating characteristic (ROC) curve (survival ROC, version 1.0.3) at 5 years in the training dataset. The Kaplan-Meier estimation method was used to estimate the ROC curve. The predictive value of HRGS that corresponded to the shortest distance between the ROC curve and the point representing both the 100% true positive rate and 0% false-positive rate was selected as the cutoff value.
The prognostic value of the HRGS was assessed in stage II/III CRC patients and patients with all stages in the training and independent validation cohorts in univariate analyses, respectively. The prognostic value of HRGS was also examined with available clinical and pathologic variables in multivariate analyses.

Functional Annotation and Analysis
To investigate the biological characteristic of the HRGS, we conducted enrichment analysis for differentially expressed genes between hypoxic risk groups in TCGA CRC dataset using the R package "gProfileR." The biological pathways of interest were further analyzed by GSEA in R using the Bioconductor package "HTSanalyzeR." (22) In addition, we estimated the proportion of infiltrated immune cells and stromal cells within the tumor tissues using the ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) (23). The proportion of different types of immune cells such as lymphocytes, monocytes, and neutrophils and necrosis percentage were calculated using CIBERSORT (24).

Statistical Analysis
Statistical analysis was conducted using SPSS (version 22.0.0, IBM SPSS Statistics, IBM Corp., Armonk, NY) and R software (version 3.5.1; http://www.Rproject.org). Means with standard deviations or medians with interquartile ranges were calculated

Discovery and Training of the HRGS
According to the CIT gene microarray dataset (GSE39582), a total of 309 eligible CRC patients without chemotherapy were enrolled in the analysis as a discovery cohort. Among 3,444 HRGs selected from MSigDB, 3,184 HRGs were detected in the meta-analysis in this study. Filtering based on conditions that the median absolute deviation is >0.5 and expression level is greater than the total median expression level, 1,636  (Figure 1). We used satisfactory RFS cutoff at 5 years in a time-dependent ROC curve analysis to train HRGS for stratifying high and low hypoxic risk groups (Supplementary Figure 1A). The prognosis correlation coefficient of each gene in HRGS was determined ( Table 1), and the related risk score calculation model was defined (Supplementary Figures 1B,C). As expected, GSEA showed that the hypoxia pathways were enriched in the high hypoxic risk patients from GSE39582, confirming that the HRGS was hypoxicrelated signature (Supplementary Figure 1D). Using stage II/III CRC patients in the CIT dataset (n = 469) as training dataset, more recurrence cases were arranged in the hypoxia high-risk group than in the low-risk group (

Validating the HRGS Prognostic Capability in TCGA and Metavalidation Cohort
To validate the prognostic power of HRGS, the TCGA dataset (n = 624) and the metavalidation cohort (n = 687) composed of GSE17536, GSE33113, GSE37892, and GSE14333 were used. All datasets include the transcriptional data and prognostic information. No significant difference was observed among the three cohorts regarding the clinical and pathologic factors ( Table 2, Supplementary Tables 1, 2). To evaluate the HRGS signature, we compared the nomogram with or without HRGS in R and analyzed the area under the curve (AUC) of the ROC curve. The AUC of nomogram with HRGS was 0.61, which was better than the nomogram without HRGS 0.56 (Supplementary Figures 2A,B). Decision curve analysis (DCA) was conducted to determine the clinical usefulness of the predicted nomogram by quantifying the net benefits at different threshold. As expected, addition of HRGS improved the DCA in patients of GSE39582, TCGA, and metavalidation cohorts (Supplementary Figure 2C).
To examine whether HRGS could predict the effect of adjuvant chemotherapy for CRC patients, we separated the nonchemotherapy-treated and the chemotherapy-treated group into high and low hypoxia groups using HRGS test and analyzed the prognosis data of these patient groups. In CRC patients without adjuvant chemotherapy, the DFS of the high hypoxia group was worse than that of the low hypoxia group in both the training (HR = 5.10, 95% CI = 2.51-10.35, P < 0.001, Figure 3A) and TCGA cohorts (HR = 2.54, 95% CI = 1.03-6.30, P = 0.037, Figure 3D). In patients who received adjuvant chemotherapy, no significant difference in DFS was found between the two groups (Figures 3B,E; P > 0.05). Moreover, the DFS of the high hypoxia group without chemotherapy treatment had a similar outcome as compared with the chemotherapy-treated patients (Figures 3C,F; P > 0.05).

Functional Annotation of the HRGS
To explore the possible underlying mechanisms responsible for the worse outcome of DFS in the high hypoxic risk patients, enrichment analysis of differentially expressed genes was conducted by identifying several overrepresented biological processes in the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Interestingly, the most valuable biological processes were found to be associated with immune system response ( Figure 4A). To further evaluate the role of genes corresponding to the HRGS, GSEA was performed in TCGA CRC cohort and found that the hypoxia environment was statistically related to interferon α (IFN-α), IFN-γ, interleukin 6 (IL-6), and IL-2 ( Figure 4B), which are important components of the immune response network.
ESTIMATE algorithm was used to further validate the immune system response between different hypoxic risk patients. A lower stromal score (P < 0.001), immune score (P < 0.001), and ESTIMATE score (P < 0.001) were found in the high hypoxic risk group (Supplementary Figures 4A,B). Further analysis of the specific immune cell types were identified using CIBERSORT, lower percentages of CD4 T cells, and M1 macrophages, whereas higher percentages of regulatory T cells and M2 macrophages were detected in the high hypoxia group (Figures 4C,D). To investigate whether the HRGS can predict the prognosis of immunotherapy, we evaluated HRGS in a cohort of patients with programmed cell death (PD-1) blockade treatment (25). The results show that high hypoxic risk patients have a tendency for worse prognosis after immunotherapy but cannot reach statistical significance (Supplementary Figure 5A).
In addition, high-risk group patients have a lower score of T-effector cell infiltration score (IMmotion150 Teff signature), immune infiltration (Javelin signature), and Merck18, which are consistent with our results (Supplementary Figure 5B).
To investigate the therapeutic strategy for high hypoxic risk patients, we analyzed the half maximal inhibitory concentration (IC 50 ) of different drugs in the Pharmacogenomics database Genomics of Drug Sensitivity in Cancer (GDSC) with HRGS. Using the same risk score calculation model, all cell lines in the GDSC were divided into low and high hypoxic risk groups and analyzed the IC 50 of different drugs. A total of five drugs have statistical significance: AZ6102, fulvestrant, irinotecan, temozolomide, and topotecan (Supplementary Figure 5C). As expected, high hypoxic risk group cells were more resistant to the traditional chemotherapeutic drugs, such as fulvestrant, irinotecan, temozolomide, and topotecan. Interestingly, high hypoxic risk group cells were more sensitive to AZ6102, a selective TNKS1/2 inhibitor.

DISCUSSION
Current cancer-related clinical trials have not included the hypoxia status as a variable factor despite large variability of tumor microenvironment due to the hypoxia status (26)(27)(28)(29). Although hypoxia tumor microenvironment has been shown to FIGURE 4 | Functional annotation of the HRGS. Enrichment analysis of the differentially expressed genes between risk groups in GO and KEGG (A). GSEA showed IFN-α, IFN-γ, IL-6, and IL-2 were depressed in the high hypoxic risk patients (B). Immune cells were assessed based on data from TCGA (C). CD4 T cells and M1 macrophages were depressed while T cells and M2 macrophages were enriched in the high hypoxia group (D,E). P values comparing hypoxia high-and low-risk groups were calculated with t tests.
affect the chemotherapy outcome in various cancer types (14,(30)(31)(32), there is no tool that could distinguish the high/low hypoxia risk and predict prognostic response to chemotherapy in CRC. In this study, we selected various HRGs to create an HRGS for CRC patients. The data suggested that the HRGS can stratify stage II/III CRC patients into subgroups with different DFS at a 5year follow-up duration. Moreover, HRGS test showed a better C-index outcome as compared with the existing prognostic tool, Oncotype DX. These results indicate that the HRGS test could be an effective prognostic tool to distinguish the hypoxia status among CRC patients.
Current TNM stage classification could not efficiently separate the patient group for whom chemotherapy was effective within stage II CRC patients. Studies have shown that chemotherapytreated patients have >5% improvement in the 5-year survival rate (2,33), and the treatment did not reduce the recurrence risk (34)(35)(36). A more efficient selection method is required to determine the patient group for whom chemotherapy was effective within stage II CRC patients. A hypoxia gene signaturebased test may be a potential candidate, as it is known that the hypoxia microenvironment could promote tumor invasion and metastasis (10,11,37), and HRGS could effectively select patients for individual-based treatment in other cancer types (7,8). However, the predictive potential of HRGS for the disease prognosis and the chemotherapy outcome in CRC patients has not been examined yet. In this study, we identified high and low hypoxic risk groups within CRC patients using HRGS.
The high hypoxic risk group within non-chemotherapy-treated CRC patients had a similar disease prognosis as compared to the chemotherapy-treated CRC patients and had a significantly worse prognosis as compared with low hypoxic risk group. These results implied that high hypoxic risk related to tumor recurrence and HRGS could help to determine the non-chemotherapytreated patients who could benefit from chemotherapy or other adjuvant treatment.
To find a novel therapeutic strategy for high hypoxia risk patients, it is important to understand the intrinsic mechanism among hypoxia-related poor prognosis. Our previous study showed that the hypoxia tumor microenvironment was related to dysregulated cell cycle machinery and PI3K-AKT-mTOR pathways (38). To date, a genome-wide mechanistic analysis of how hypoxia tumor microenvironment induces poor chemotherapy response is still lacking. Analysis of an HRGS test identified that the hypoxia high-risk group has significantly lower scores of stromal and immune cell infiltrations than those observed in the hypoxia low-risk group. Previous studies have postulated that tumor hypoxia reduces the antitumor effect by suppressing the microenvironment immune response (37,39). A large number of immunosuppressive cells, such as myeloid-derived suppressor cells, tumor-associated macrophages, and T-regulatory cells, were found in the hypoxic zones of the tumor (37,40,41). In line with these data, we also found that the enriched T-regulatory cells and M2 macrophages, but less amount of M1 macrophages in the high-hypoxic risk patients. These data suggest that hypoxia microenvironment was probably associated with immune suppression and enhanced tumor progression. As hypoxic microenvironment could provide tumor resistant to immunotherapy (42,43), which has emerged as an effective therapy for CRC (44,45), HRGS may have the potential to predict the outcome of the immunotherapy. Further studies will be required to verify this hypothesis.
Despite the exciting finding in this study, some limitations to the study should be noted. First, although we included as many datasets as possible for validation of our HRGS, the results of this study should be carefully evaluated because of the shortage of the retrospective design. Second, gene expression signature is subject to sampling bias caused by intratumor genetic heterogeneity (46). Although we have reduced as many as possible cross-study batch effects by constant ordering and excluding HRGs, their complex nature implies that not all batch effects can be addressed. To be more objective, further validation using prospective data from multiple centers would be ideal and necessary before its application in clinical practice.
In conclusion, the proposed prognostic HRGS is a promising test system to estimate the DFS of stage II/III CRC patients and to predict the possible beneficial effect from the chemotherapy. Hypoxia tumor microenvironment may promote chemoresistance in stage II/III CRC patients by inhibiting the local immune response.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board (IRB) of The Sixth Affiliated Hospital of Sun Yat-sen University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
Y-fC, Z-lY, M-yL, PL, X-jW, and FG contributed to study concept and design, acquisition, analysis, interpretation of data, and drafting of the manuscript. Y-fC, Z-lY, BZ, X-hL, Z-rC, M-yL, and Y-xT contributed to data collections and manuscript review. Y-fC, Z-lY, Y-fZ, and JK contributed to study concept and design, analysis and interpretation of data, and critical revision of the manuscript for important intellectual content. PL, X-jW and FG supervised the study. All authors read and approved the final manuscript.