Effects of Hypoxia in Intestinal Tumors on Immune Cell Behavior in the Tumor Microenvironment

Background Imbalanced nutritional supply and demand in the tumor microenvironment often leads to hypoxia. The subtle interaction between hypoxia and immune cell behavior plays an important role in tumor occurrence and development. However, the functional relationship between hypoxia and the tumor microenvironment remains unclear. Therefore, we aimed to investigate the effect of hypoxia on the intestinal tumor microenvironment. Method We extracted the names of hypoxia-related genes from the Gene Set Enrichment Analysis (GSEA) database and screened them for those associated with colorectal cancer prognosis, with the final list including ALDOB, GPC1, ALDOC, and SLC2A3. Using the sum of the expression levels of these four genes, provided by The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, and the expression coefficients, we developed a hypoxia risk score model. Using the median risk score value, we divided the patients in the two databases into high- and low-risk groups. GSEA was used to compare the enrichment differences between the two groups. We used the CIBERSORT computational method to analyze immune cell infiltration. Finally, the correlation between these five genes and hypoxia was analyzed. Result The prognosis of the two groups differed significantly, with a higher survival rate in the low-risk group than in the high-risk group. We found that the different risk groups were enriched by immune-related and inflammatory pathways. We identified activated M0 macrophages in TCGA and GEO databases and found that CCL2/4/5, and CSF1 contributed toward the increased infiltration rate of this immune cell type. Finally, we observed a positive correlation between the five candidate genes’ expression and the risk of hypoxia, with significant differences in the level of expression of each of these genes between patient risk groups. Conclusion Overall, our data suggest that hypoxia is associated with the prognosis and rate of immune cell infiltration in patients with colorectal cancer. This finding may improve immunotherapy for colorectal cancer.


INTRODUCTION
Since Stephen Paget proposed the "seed and soil" theory of cancer development and tumor metastasis (1), the understanding of the tumor microenvironment has gradually deepened. Several components of the tumor microenvironment contribute toward tumor occurrence and development (2,3). An imbalance between nutrient supply and demand within the tumor often leads to hypoxia, glucose deficiency, and consequently an acidic tumor microenvironment (4). Specifically, tumor cells can use immune escape mechanisms to drive metastasis and invasion in anoxic environments (5). This has increased research interest into the relationship between hypoxia and the tumor immune microenvironment.
Hypoxia can reduce the activity of various immune cells in the tumor microenvironment and the production of corresponding immune stimulators to increase the release of suppressors and the expression of immune checkpoint inhibitors (5,6), suggesting a close relationship between hypoxia and immune cells.
This study aimed to investigate the relationship between hypoxia-related genes and the immune microenvironment of colorectal cancer by (1) quantifying the influence of these genes on the tumor immune microenvironment, (2) screening hypoxiarelated genes associated with intestinal tumor prognosis, (3) developing a hypoxia risk score model, (4) validating the model's ability to predict patient prognosis and risk using data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, (5) identifying the enrichment of model-related pathways, (6) using the hypoxia risk score as the entry point to explore differences in the infiltration rate of immune cells, and (7) identifying genes that affect immune cell infiltration to evaluate the relationship between hypoxia risk score and the tumor immune microenvironment.

Raw Data
The intestinal cancer-related RNA-sequencing and clinical data used in this study were obtained from TCGA (https://portal.gdc. cancer.gov/) and GEO (gse39582) databases.

Construction and Grouping of the Hypoxia Model
The screening of hypoxia-related genes identified the genes independently related to the prognosis of intestinal cancer. We then multiplied the expression levels of each gene in TCGA and GEO databases by their respective expression coefficients, with the resulting sum defined as the risk score for each patient. Based on the median risk score value, the patients in the two databases were divided into high-and low-risk groups for the follow-up evaluations.

Survival Analysis
The survival and survminer R packages (The R Foundation for Statistical Computing, Vienna, Austria) were used to analyze the prognosis of 445 and 579 patients in TCGA and GEO databases, respectively. Patients in TCGA database were followed up for 12 years, while those in the GEO database were followed up for 16 years. The Kaplan-Meier method was used to plot the survival curves, and the log-rank test was used to evaluate the statistical significance between them, with p-values <0.05 considered statistically significant.

ROC Curve Analysis
ROC curve analysis was conducted using the survival, survminer, and timeROC packages in R. The 1-, 3-, and 5-year survival rates of patients in TCGA and GEO databases were evaluated. The area under the ROC curve for the 1-, 3-, and 5-year survival rates increased gradually and exceeded 0.5, which was defined as the threshold for the accurate prediction of survival by the model. Survival for each group is shown as risk columns and risk curves in Figures 2A, B.

Heat Maps
In this study, heat maps of gene expression were drawn using the pheatmap package in R.

PPI Network Analysis
A network of anoxic genes was constructed using the STRING database. R software was used to select the 50 genes with the largest numbers of adjacent nodes for subsequent analysis.
provides the names of all hypoxia-related genes. We then used the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (http://string-db.org/cgi/input.pl), a protein-protein interaction (PPI) network database, to construct a model of PPIs between hypoxia-related genes ( Figure 1A). By counting the number of adjacent nodes for each protein, we identified the core genes. The top 50 core genes with adjacent nodes are presented in Figure 1B. We then downloaded the intestinal tumor gene expression data and associated clinical information available from TCGA database. We extracted the data on hypoxiarelated gene expression and screened genes associated with prognosis using univariate Cox regression analysis ( Figure 1C). Among the identified genes, ALDOB was associated with a low risk, while five genes (GPC1, ALDOC, ENO2, SERPINE1, and SLC2A3) were associated with a high risk of developing tumor malignancy. The multivariate Cox regression analysis of genes associated with patient prognosis revealed four genes ( Figure 1D) that were then used to construct prognostic models. These four genes (ALDOB, GPC1, ALDOC, and SLC2A3) had model coefficients of −0.1574, 0.2994, 0.2647, and 0.2074, respectively.

Effects of Hypoxia-Related Genes on Prognosis
Multivariate Cox regression analysis identified four hypoxiarelated genes associated with intestinal tumor prognosis, which were used for modeling. We multiplied the expression levels of these genes (as reported in TCGA and GEO databases) by the corresponding coefficients to obtain the risk score for each patient. The median risk score value was then used to divide the patients in TCGA database into high-and low-risk groups. The subsequent survival analysis revealed significant differences between the high-and low-risk groups (p < 0.05; Figures 2A, B). included in TCGA database ( Figure 2C). The area under the ROC curve for the GEO database was > 0.05 ( Figure 2D), indicating that the model accurately predicted the survival rate. To value the survival rate of patients more intuitively in the highand low-risk groups, we used risk histograms to show the differences in survival status between the two databases ( Figures  3A, B). We observed a higher proportion of surviving patients in the low-risk group than in the high-risk group. These results further demonstrated that our model effectively distinguished between high-and low-risk patients. We also analyzed the interactions between the hypoxia-related genes that were identified as affecting patient prognosis in the model ( Figures  3C, D). We further demonstrated the relationship between patient risk and survival using risk curves, in which the risk scores for both groups of patients were plotted ( Figures 3E, F). We observed a longer survival time in the low-risk group than in the high-risk group. Moreover, the number of deaths in the low-risk group decreased over time ( Figures 3G, H). Finally, we compared the expression level of each gene included in the model between the high-and low-risk groups using thermography ( Figures 3I, J).

Effects of Different Clinical Characteristics on Intestinal Tumor Prognosis
Clinical characteristics differ in the impact they have on patient prognosis. Thus, we analyzed the impact of clinical characteristics on the prognosis of patients included in TCGA and GEO databases. We first used univariate Cox regression analysis to evaluate the impact of clinical characteristics on the survival time and prognosis of patients included in the two databases ( Figures  4A, B). We found that patient sex affected neither survival nor prognosis, while other factors affected survival and prognosis and were associated with increased risk. However, the p-value of the risk score in our model was <0.05 for patients in TCGA database, indicating that the risk score also affected patient prognosis and survival. In contrast, the p-value was >0.05 for patients in the GEO database. Multivariate analysis of these factors showed that the pvalues for age and tumor-node-metastasis (TNM) stage were both <0.05, indicating that these variables were independent prognostic factors ( Figures 4C, D). We also observed differences in the expression levels of hypoxia-related genes between patients with different T stages included in the two databases ( Figures 4E-H). The expression level of SLC2A3 in different T stages (T1, T2, T3, and T4) differed significantly between patients in TCGA and GEO databases (p < 0.05).

Enrichment of Pathways in Hypoxia-Related Risk Groups
We observed differences in the enrichment level of hypoxiarelated genes and associated pathways between the high-and low-risk groups. To understand the level of pathway enrichment, we used GSEA software (UC San Diego, San Diego, CA, USA; Broad Institute, Cambridge, MA, USA) to compare the pathways between different risk groups. The high-risk group in TCGA database showed many enriched pathways related to apoptotic and immune functions ( Figure 5A) compared with the low-risk group. Among these were apoptosis, hypoxia, interleukin (IL) 2/ signal transducer and activator of transcription (STAT) 5 signaling, IL6-STAT3 signaling, and the inflammatory response pathways. The other related pathways are shown in Table 1. In contrast, the pathways enriched in the low-risk group mainly included those associated with oxidative phosphorylation and lipid metabolism, among others. The pathways with false discovery rate (FDR) q values >0.05 are shown in Table 2. The low-risk group of patients in the GEO database also showed the enrichment of apoptosis and IL2-STAT5 signaling pathways, in addition to, the p53 and peroxisome and phosphatidylinositol 3kinase-mammalian target of rapamycin serine/threonine protein kinase B signaling pathways ( Figure 5B). The enrichment of other pathways is shown in Table 3. The high-risk group also showed the enrichment of the Hedgehog and Wnt/beta-catenin signaling pathways. Pathways with FDR values >0.05 are shown in Table 4.

Immune Cell Infiltration
The results of the pathway enrichment analysis of hypoxia-related genes for both risk groups in TCGA and GEO databases showed the enrichment of pathways related to inflammation, immunity, and other factors. Based on these findings, we assessed the rate of immune cell infiltration in each risk group according to the constructed hypoxia-related gene model. Figures 5C, D shows the infiltration rate of immune cells in the high-and low-risk groups in TCGA and GEO databases, respectively. TCGA database showed differences in the infiltration rate of two types of immune cells in the high-and lowrisk groups ( Figure 5E; p < 0.05). In contrast, the rate of infiltration of 10 types of immune cells differed between the high-and low-risk groups in the GEO database ( Figure 5F; p < 0.05). Among them, the infiltration rate of activated M0 macrophages differed between the high-and low-risk groups in both databases. Using the Tracking Tumor Immunophenotype online platform (http://biocc.hrbmu.edu. cn/TIP/index.jsp), we screened the immune-related genes for those that played important roles in the regulation of this immune cell type. Heat maps were then drawn to visualize the expression of these genes relative to that of hypoxia-related genes in the high-and low-risk groups of patients in TCGA and GEO databases ( Figures 6A, B).
The expression levels of CCL2/4/5, CSF1, and CX3CL1 were significantly different between the high-and low-risk groups in both databases (p < 0.05). Next, we plotted the correlation curves between the expression levels of these four genes and the risk scores, which showed that the expression levels of these four genes were positively correlated with the patients' risk score. We also observed differences in the expression levels of these four genes between the high-and low-risk groups (Figures 6C-L). (F) Immune cells whose infiltration is significantly associated with the risk of hypoxia in the GEO database (p < 0.05).

DISCUSSION
Hypoxia is a feature of tumor physiology specifically that of mechanisms associated with the acquisition of some malignant attributes, such as metastasis, invasion (7-9), and drug resistance (10). In these processes, hypoxia-related genes act on the corresponding pathways or on the regulating immune cells. Hypoxia-inducible factor (HIF) is activated during hypoxia. Some immunosuppressive factors, such as vascular endothelial growth factor, are HIF target genes, which affect both angiogenesis and immunosuppression (11). Hypoxia also results in upregulated EGFR expression, which promotes ligand-independent epidermal growth factor receptor signaling (12,13). This process increases the rate of tumor glycolysis, resulting in metabolic competition (14). In our study, we screened for hypoxia-related genes in the gut and found that the core genes (ALDOB, GPC1, ALDOC, and SLC2A3) were closely related to patient prognosis. The rate of aldolase-B and fructose-bisphosphate B-driven fructose metabolism is significantly increased in patients with colon cancer and liver metastasis, and in those with colorectal villous polyps (15,16). GPC1 has also been shown to be overexpressed in various malignancies (17,18). A recent study has reported GPC1 enrichment in tumor-derived exosomes (17). In melanoma, NME1 has been shown to inhibit metastasis by activating ALDOC transcription (19). Solute carrier family 2, member 3 can increase glucose uptake in anoxic cells and, thus, increase the rate of glycolysis. These findings suggest that hypoxia-related gene expression levels are closely related to tumor development and metabolism (20); therefore, they were included in our hypoxia-related gene model. With a deepening understanding of the mechanisms of hypoxia, its influence on tumor prognosis is also increasingly being understood (21,22). To further investigate the relationship between the expression of hypoxia-related genes and patient prognosis, we evaluated patient prognosis by taking the product sum of the expression levels and the coefficients of ALDOB, GPC1, ALDOC, and SLC2A3 in TCGA and GEO databases as the risk score. The prognosis of patients in the high-and low-risk groups differed significantly; the survival rate of patients in the low-risk group was significantly higher than that of their counterparts. However, the ROC curve analysis, which aimed to evaluate the accuracy of the survival estimates, showed that the curves obtained from the GEO database corresponded poorly to the observed prognosis. Thus, we also analyzed the influence of other factors, including age, sex, TNM staging, and our proposed risk score, on patient prognosis. We concluded that our proposed risk score showed a good correspondence with patient prognosis. However, the results of the multivariate analysis in the GEO database were not statistically significant. This finding suggests that the risk score alone cannot be used as an independent prognostic factor. The samples in the GEO database were all colorectal adenocarcinoma samples that had been obtained in France. Because the tumor type and region are very specific, these samples may not accurately reflect the relationship between hypoxia and colon cancer prognosis. Moreover, the prognosis of colon cancer is related to disease stage and clinical, histological, genetic, and molecular factors, among others. These factors should be considered in future studies. We next applied GSEA to identify pathways enriched in the high-and low-risk groups in the two databases. We found that most of the enriched pathways were related to inflammation, immune response, and apoptosis. Hypoxia and cell death in tumor tissue produce large amounts of cell debris and trigger the release of inflammatory factors, which can attract macrophages and monocytes, and can induce macrophage polarization. After polarization, macrophages secrete inflammatory factors (23). These findings suggest a close relationship between hypoxia, inflammation, and the immune response. In addition, some apoptosis-related genes, such as p53, which is a tumor suppressor gene, are closely related to tumor apoptosis. Nilay et al. (24) found that mutations in p53 in gastric and esophageal cancer cells can induce hypoxia signaling. This finding was confirmed in a study involving nude mice.
The results of many previous studies and those of the present study have shown that hypoxia can recruit immune cells into the tumor microenvironment. Hypoxia-induced tumor-derived cytokines, such as IL-10 and transforming growth factor-beta, can induce tumor-associated macrophages to differentiate into M2 macrophages with immunosuppressive effects (25). When monocytes are stimulated by inflammatory factors, such as interferon-gamma and lipopolysaccharide, they activate M1 macrophages, which can secrete inflammatory factors, such as IL-6 and tumor necrosis factor-alpha, and can phagocytize invasive pathogens and tumor cells (23,26). Hypoxia can also lead to immune escape through the role of immune cells (5); for example, hypoxia can reduce T-cell activity (27). Hypoxia is closely related to immune cell function. In our study, we observed significant differences in the level of activated M0 macrophages between the high-and low-risk groups in TCGA database. In the GEO database, the level of infiltration of activated CD4 memory T cells and M0 macrophages, activated and resting mast cells, and neutrophils differed significantly between the high-and low-risk groups. Both databases showed significant differences in the infiltration rate of activated M0 macrophages between the high-and low-risk groups. When we screened the genes that regulated this immune cell type, we found different levels of CCL2/4/5 and CSF1expression between the high-and low-risk groups. Further analyses confirmed the significant correlation between the expression level of each of these genes and the risk of hypoxia. Three of these four genes encode chemokines, which play a chemotactic role in immune cells, such as natural killer cells and monocytes, which are closely related to tumor development. CCL2 and CCL5 play important roles in prostate cancer metastasis and drug resistance (28,29), while CCL4 is associated with the clinical characteristics of breast cancer (30). CSF1 is expressed in almost all tumors (31,32). It recruits macrophages other than alveolar macrophages through CSF1 receptors and regulates their differentiation (31,33). Macrophages, a major component of the tumor microenvironment, contribute to tumorigenesis by promoting angiogenesis, immunosuppression, invasion, and metastasis (34,35).
In conclusion, hypoxia plays an important role in the tumor microenvironment. Screening for genes that can affect the rate of immune cell infiltration revealed a correlation between these genes and the hypoxia risk score. Our findings show that hypoxia-related genes can affect the prognosis of intestinal cancer and may play a role in immune infiltration in intestinal cancer. Analysis of the relationship between hypoxia and immune cells may improve immunotherapy and tumor treatment.