Screening and Validation of the Hypoxia-Related Signature of Evaluating Tumor Immune Microenvironment and Predicting Prognosis in Gastric Cancer.

Background: Hypoxia is one driving factor of gastric cancer. It causes a series of immunosuppressive processes and malignant cell responses, leading to a poor prognosis. It is clinically important to identify the molecular markers related to hypoxia.
Methods: We screened the prognostic markers related to hypoxia in The Cancer Genome Atlas database, and a risk score model was developed based on these markers. The relationships between the risk score and tumor immune microenvironment were investigated. An independent validation cohort from Gene Expression Omnibus was applied to validate the results. A nomogram of risk score model and clinicopathological factor was developed to individually predict the prognosis.
Results: We developed a hypoxia risk score model based on SERPINE1 and EFNA3. Quantified real-time PCR was further applied to verified gene expressions of SERPINE1 and EFNA3 in gastric cancer patients and cell lines. A high-risk score is associated with a poor prognosis through the immunosuppressive microenvironment and immune escape mechanisms, including infiltration of immunosuppressive cells, expression of immune checkpoint molecules, and enrichment of signal pathways related to cancer and immunosuppression. The nomogram basing on the hypoxia-related risk score model showed a good ability to predict prognosis and high clinical net benefits. Conclusions: The hypoxia risk score model revealed a close relationship between hypoxia and tumor immune microenvironment. The current study potentially provides new insights of how hypoxia affects the prognosis, and may provide a new therapeutic target for patients with gastric cancer.

Introduction

Gastric cancer is a public health burden, ranking fifth in global incidence and fourth in mortality among all cancers (1). Therapeutic strategies are still based on the American Joint Committee on Cancer (AJCC) tumor/node/metastasis (TNM) staging system (2, 3). However, due to the high heterogeneity of gastric cancer, patients with similar clinicopathological characteristics could have different prognosis, suggesting the current TNM staging system are inadequate for predicting prognosis and risk stratification (4, 5). Therefore, it is clinically important to develop a novel biomarker to better guide clinical treatment and improve prognosis.

Tumor cells always grow faster than their blood vessels. Owing to the inadequate blood supply, the supply of oxygen and nutrients to the tumor cells is unbalanced, thereby forming a hypoxic microenvironment (6–9). Hypoxia is one of the characteristics of tumor microenvironment (TME) that can lead directly to the malignant characteristics, including tumor proliferation, migration and invasion, resulting in a poor prognosis (10–12). Previous studies have shown a significant relationship between hypoxia and poor prognosis of GC (13, 14), and hypoxia plays a key role in metastasis (15). In the hypoxic microenvironment, hypoxia-inducible factors (HIFs) are key transcription factors that allow cancer cells to survive under hypoxic conditions and promote tumor progression (16–18). Multiple genes transcribed by HIFs, including Glut1, KLF8, VEGFA, ITGβ1, etc., can promote GC metastasis and lead to poor prognosis (19).

TME is the internal environment in which tumor cells are produced and survive. It is composed of immune cells, endothelial cells, mesenchymal cells, inflammatory mediators and extracellular matrix molecules (20, 21). The immunological components of the TME can inhibit or promote tumor development (22). Recently, the significance of hypoxia in promoting tumor immunosuppression and immune escape has received increasing attentions (23, 24). It is important to understand the potential mechanisms that are involved between hypoxia and the tumor immune microenvironment. Therefore, the establishment of a hypoxia-based signature may help to identify the potential prognostic value of hypoxia, and improve the comprehension of the immunogenomic profile of gastric cancer.

Here, we established a hypoxia-related signature related to prognosis by The Cancer Genome Atlas (TCGA) data base, which was validated by the Gene Expression Omnibus (GEO) data base. Potential mechanisms of the hypoxia-related signature were further investigated.

Materials and Methods

Patients

The Clinical data (375 cancer and 72 non-cancerous samples) and FPKM RNA-seq data from TCGA data base (https://www.cancer.gov/tcga) was applied as a screening cohort. The data of 433 cancer samples (GSE84437) from GEO data base (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi) was applied as a validation cohort (25). RNA-seq and microarray data included were transformed [log2(x+1)] and normalized by the “sva” and “limma” packages of R software. The baseline characteristics of the screening and validation cohorts are shown ( ).

Development of a Risk Score Model

Univariate analysis was firstly applied to identify potentially hypoxia-related genes that have a statistically significant difference of prognosis in gastric cancer patients from TCGA data base. Least absolute shrinkage and selection operator (LASSO) method was then applied to shrink the scope of gene screening (26). Finally, Cox proportional hazards analysis was used to identify highly hypoxia-related genes. The risk score formula was constructed as: Risk score = (∑coefficientx * expression of signature genex) (genex indicated the identified genes). The regression coefficient was obtained from Cox proportional hazards analysis. The patients of gastric cancer were divided into a high-risk and a low-risk groups by the cut-off value of the median risk score.

Tumor Immune Microenvironment

To investigate the relationships between risk score and TME, the ESTIMATE algorithm was applied to determine immune score, stromal score, ESTIMATE score, and tumor purity of individual patient in the screening and validation cohorts (27). Wilcoxon test was applied to compare the differences between the high-risk and low-risk groups in terms of immune score, stromal score, ESTIMATE score, and tumor purity. The TIMER web server (http://timer.cistrome.org/) was applied to analyze the correlations between signature genes and immune cells. The TIMER algorithm was used to assess the abundances of six immune infiltration cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) and tumor purity (28).

The ssGSEA method was used to the transcriptome to assess immune cell infiltrations (29). We obtained a set of marker genes including immune cell types, immune-related pathways and functions (30). We used the R package called “GSVA” to perform ssGSEA to obtain the normalized enrichment score (NES) of each immune-related item.

Development and Assessment of a Predictive Nomogram

Univariate analysis and Cox proportional hazards analysis were conducted on risk score of target genes and patient clinicopathological characteristics to determine independent prognostic factors related to prognosis. The predictive nomograms were developed by including all independently prognostic factors.

GC Cell Lines and Tissue Samples

The human gastric epithelial cell line GES-1 and gastric cancer cell lines AGS, SGC-7901, HGC-27, MKN-45 and MGC-803 were purchased from the Chinese Academy of Sciences (Shanghai, China). Cells were cultured in RPMI 1640 medium (HyClone, Logan, UT, USA) with 10% fetal bovine serum (FBS, Invitrogen) and 1% penicillin/streptomycin in a humidified atmosphere of 5% CO2 at 37°C. Totally, 39 pairs of gastric cancer together with their adjacent non-cancerous tissues (> 5 cm away from cancer tissue) were collected. This study was approved by the Ethics Committee of the Fourth Affiliated Hospital of China Medical University (EC-2021-KS-043). All patients included in this study provided written informed consent in accordance with the Declaration of Helsinki.

Quantitative Real-Time PCR Analysis

Total RNA was extracted using Trizol reagent (Invitrogen, Eugene, OR) and was used to synthesize cDNA using Prime-Script RT Master Mix (TaKaRa, Shiga, Japan), and quantitative real-time PCR (qRT-PCR) was performed by TaKaRa SYBR® Premix Ex Taq™ (TaKaRa, Shiga, Japan). All primers of qRT-PCR were listed in .

Statistical Analyses

All analyses were performed by R version 4.0.2 (http://www.R-project.org). The calibration curve and area under the curve (AUC) were used to evaluate the predictive performance of the predictive nomogram. The clinical benefit was further evaluated by the decision curve analysis (DCA) (31, 32). An independent validation cohort was applied to validate these findings. All tests were two-sided, and a P value less than 0.05 was considered as statistically significant.

Results

Patient Characteristics

The baseline characteristics of the screening and validation cohorts are shown in . In the screening cohort, a total of 234 (63.1%) patients were male and 137 (36.9%) were female. Among them, 17 (4.6%) patients were T1, 74 (19.9%) were T2, 177 (47.7%) were T3, and 103 (27.8%) were T4 cases. 117 (31.5%) patients were N0, 97 (26.1%) were N1, 79 (21.3%) were N2, and 78 (21.1%) were N3 cases. Accordingly, 8 (2.2%) patients were Grade I, 126 (34.0%) were Grade II, and 237 (63.8%) were Grade III. Considering the TNM staging system, 46 (12.4%) patients were stage I, 119 (32.1%) were stage II, 165 (44.4%) were stage III, 41 (11.1%) were stage IV cases.

In the validation cohort, a total of 296 (68.4%) patients were male and 137 (31.6%) were female. Among them, 11 (2.5%) patients were T1, 38 (8.8%) were T2, 92 (21.2%) were T3, and 292 (67.5%) were T4 cases. Accordingly, 80 (18.5%) patients were N0, 188 (43.4%) were N1, 132 (30.5%) were N2, and 33 (7.6%) were N3 cases.

Screening of Hypoxia-Related Risk Signature in Gastric Cancer

The hallmark hypoxia-related 200 genes, was obtained from the Molecular Signatures data base (MSigDB version 6.0). Among them, the TCGA data base contains 197 hypoxia-related genes, and 41 differentially expressed genes (DEGs) have been identified ( and ). To better visualize the interactions between these hypoxia genes, the STRING online data base was used to analyze the protein-protein interaction network ( ). We evaluated the hypoxia-related genes in the screening cohort, and identified 14 of the 41 genes that were significantly associated with prognosis (all P < 0.05, ). The LASSO method was further used to analyze these 14 genes, which minimized the potential over-fitting problem and established the minimum standard. Five of the 14 genes in the model are under the optimal adjustment parameter (λ) ( ). Finally, the Cox proportional hazards analysis confirmed that two genes (SERPINE1, EFNA3) ( ) met the proportional hazard hypothesis and were finally used to establish the following risk score model: Risk score = (0.223 * expression level of SERPINE1) + (–0.165 * expression level of EFNA3). Of the two signature genes, SERPINE1 was a risk DEG, and EFNA3 was protective. The risk score of individual patient was calculated and all patients were classified into a high-risk and a low-risk groups based on the median risk score.

Figure 1

Identification of the hypoxia risk signature. (A) The Venn diagram shows the hypoxia-related genes in TCGA. (B) The Volcano plot for differentially expressed genes (DEGs) in cancer and non-cancer tissues. (C) The heatmap plot for DEGs in cancer and non-cancer tissues. (D) The PPI network visualizes the interaction between these DEGs. (E, F) The LASSO method identified five genes associated with prognosis. (G) The Cox proportional hazards analysis identified the hypoxia risk signature.

Prognostic Ability of Hypoxia-Related Risk Score in Gastric Cancer

Hypoxia usually promotes an aggressive tumor phenotype, so the prognostic ability of the hypoxia-related risk score was explored. In the high-risk group of the screening cohort, the heatmap showed that the expression of SERPINE1 was up-regulated and EFNA3 was down-regulated ( ). The mortality rate in the low-risk group was significantly lower than that in the high-risk group ( ). Kaplan–Meier analysis indicated that the prognosis of the low-risk group was significantly superior than that of the high-risk group (log-rank test, P < 0.001) ( ). Similar results were found in the validation cohort ( ).

Figure 2

Prognostic value of the hypoxia risk signature in gastric cancer. (A, E) Heatmaps of the prognostic signature in the screening (TCGA) and validation (GEO) cohorts. (B, F) Patient risk score in the screening and validation cohorts. (C, G) The status distribution of patients in the high-risk and low-risk groups in the screening and validation cohort. (D, H) Kaplan-Meier analysis of patients in the high-risk and low-risk groups in the screening and validation cohorts.

Hypoxia-Related Signaling Pathways

In the screening cohort, we used GSEA to analyze the signaling pathways activated in the hypoxia-related high-risk group. In the high-risk group, the JAK-STAT signaling pathway, NOTCH signaling pathway, pathway in cancer, and TGF-β signaling pathway were activated ( ). These signaling pathways are related to the stimulation of tumor proliferation, migration, invasion, anti-apoptosis, Epithelial-Mesenchymal Transition (EMT), immune escape and drug resistance. These results have been confirmed in the independent validation cohort ( ).

Figure 3

Enrichment of pathways related to hypoxia and analysis of tumor immune microenvironment. (A, B) The enrichment plots show the signaling pathways related to hypoxia in the screening and validation cohorts. (C, D) The heatmaps show 29 immune-related gene sets, immune score, stromal score, ESTIMATE score and tumor purity in the screening and validation cohorts. (E, F) The relationship between risk score and immune score, stromal score, ESTIMATE score, and tumor purity in the screening and validation cohorts.

The Correlation Between Risk Score and TME

The ESTIMATE analysis showed that the immune score, stromal score and ESTIMATE score were significantly positively correlated with the risk score in both the screening and validation cohorts, while tumor purity was significantly negatively correlated with the risk score ( ). It also indicated that the immune score, stromal score and ESTIMATE score of the high-risk group were significantly higher than those of the low-risk group (P < 0.001), while the tumor purity of the high-risk group was significantly lower than that of the low-risk group (P < 0.001, ).

The Correlation Between Risk Score and Immune Cell Subtypes

As the tumors of the high-risk group were proved to be infiltrated with a large number of immune cells, we further analyzed the subtypes of infiltrating immune cells. It indicated that the levels of immune cell infiltration in the high-risk group, including regulatory T cells, macrophages, neutrophils, and mast cells, were higher than those in the low-risk group (P < 0.05, ). Accordingly, the high-risk group reflects the immunosuppressive tumor microenvironment, full of immunosuppressive cells, which is consistent with the poor prognosis of the high-risk group.

Figure 4

Correlation of the risk score with immune cell subtypes in the screening and validation cohorts. (A, E) Regulatory T cells; (B, F) Macrophages; (C, G) Neutrophils; (D, H) Mast cells.

Then, TIMER was applied to evaluate the correlation between the expression levels of EFNA3 and SERPINE1 with tumor purity and infiltrating levels of immune cells ( ). It showed a correlation between immune cell infiltration and the expression levels of EFNA3 and SERPINE1. It showed that the risk gene SPERPINE1 had a significant positive correlation with the infiltration of macrophages, neutrophils and dendritic cells, and a significant negative correlation with tumor purity and B cells (P < 0.05, ). However, the prognostic protective gene EFNA3 showed the opposite trend for most aspects except B cells ( , all P < 0.05).

The Correlation Between Risk Score and Immune Checkpoint Molecules

We compared the immune checkpoint molecules between the high-risk and low-risk groups. The expression levels of many immune checkpoint molecules were higher in the high-risk group than those in the low-risk group ( ). In the screening cohort, the expression level of five key immune checkpoint molecules (PD-1, PD-L1, CTLA-4, HAVCR2 and TGF-β) in the high-risk group was significantly higher than those in the low-risk group, and significantly positively correlated with risk score ( ). Similar results were obtained in the validation cohort, except that CTLA4 and PD-L1 showed no significant difference between the high-risk and low-risk groups ( ).

Figure 5

Relationships between hypoxia risk score and immune checkpoint molecules. (A, G) Heatmaps show the expression level of immune checkpoint molecules in high-risk and low-risk groups in the screening and validation cohorts (*P < 0.05; **P < 0.01; ***P < 0.001). Scatter plots and box plots show the relationship between the risk score and the expression level of (B, H) PD-1, (C, I) HAVCR2, (D, J) TGF-β, (E, K) PD-L1, and (F, L) CTLA4 in the screening and validation cohorts.

The Correlation Between Risk Score and Tumor Mutation Burden and Somatic Mutation

It showed that the risk score was significantly negatively correlated with tumor mutation burden (TMB) (R = –0.36, P < 0.001; ). We further compared the TMB of patients in the low-risk and high-risk groups. It showed that the TMB of the low-risk group was significantly higher than that of the high-risk group (Wilcoxon test P < 0.001) ( ). We determined the optimal cutoff value of TMB (cutoff value = 0.68) by using the minimum P-value method, and divided the patients into a high TMB group (n = 320) and a low TMB group (n = 42). It showed that patients in the high TMB group had a better survival prognosis than those in the low TMB group (log-rank test, P < 0.001, ). We further evaluated the synergistic effect of the TMB grouping and the risk score grouping in the prognostic stratification. It showed that TMB status did not affect the survival prognosis prediction based on the risk score group. The risk score subgroup indicated significant survival differences in both the low and high TMB subgroups (log rank test, high TMB & high-risk vs. high TMB & low-risk, P < 0.001; low TMB & low-risk vs. low TMB & low-risk, P < 0.001; ). Moreover, the high TMB & low-risk group had the best overall survival rate, and the low TMB & high-risk group had the worst overall survival rate.

Figure 6

The correlation between the risk score and somatic variants. (A) The scatter plot depicts the negative correlation between risk score and tumor mutation burden (TMB) in the screening cohort. (B) TMB difference in the high-risk and low-risk groups. (C) Kaplan-Meier curves for high-risk and low-TMB groups of the screening cohort. (D) Kaplan-Meier curves for patients in the screening cohort stratified by both risk score and TMB. (E, F) Waterfall plots display the frequently mutated genes in low-risk and high-risk groups in the screening cohort. The left panel shows the genes ordered by their mutation frequencies. The right panel presents different mutation types.

Furthermore, we estimated somatic variations in gastric cancer driver genes between the low-risk and high-risk subgroups. We used Maftools to access gastric cancer driver genes and further analyzed the top 20 ones with the highest mutation frequency ( ). The results showed that there were significant differences in the mutation frequency of PCLO, TTN, FLG, LRP1B, KMT2D, SYNE1, RYR2, OBSCN, CSMD1, FAT3, ARID1A, ZFHX4, FAT4 and SPTA1 in the high-risk and low-risk groups (Chi-square test, all P < 0.05; ).

The Correlation Between Risk Score and Chemotherapeutic Drugs

We further analyzed the association between the risk score and the efficacy of chemotherapy in the treatment of gastric cancer. It showed that the high-risk group was associated with lower half inhibitory centration (IC50) of chemotherapeutic drugs, such as axitinib (P = 0.0053), bexarotene (P < 0.001), bortezomib (P < 0.001), bryostatin.1 (P = 0.0067), dasatinib (P < 0.001), imatinib (P < 0.001), midostaurin (P < 0.001), nilotinib (P = 0.04), pazopanib (P = 0.0024), sunitinib (P < 0.001), temsirolimus (P < 0.001), and vinblastine (P = 0.031), while the IC50 of methotrexate (P = 0.019) and mitomycin.C (P = 0.0035) was higher, indicating that the risk scores can be used as a potential predictor of chemical sensitivity ( ).

Figure 7

The correlation between low-risk and high-risk groups and chemotherapeutics. Sensitivity to chemotherapeutic drugs is expressed by the half inhibitory centration (IC50) of chemotherapeutic drugs. (A) Axitinib; (B) Bexarotene; (C) Bortezomib; (D) Bryostatin.1; (E) Dasatinib; (F) Imatinib; (G) Midostaurin; (H) Nilotinib; (I) Pazopanib; (J) Rapamycin; (K) Sunitinib; (L) Temsirolimus; (M) Vinblastine; (N) Methotrexate; (O) Mitomycin.C.

The Correlation Between Risk Score and Clinicopathological Characteristics

We conducted correlation analyses between clinicopathological factors and risk score in the screening cohort ( ), and tumor grade and T stage were significantly associated with risk score ( ). In the validation cohort, age, T stage, and N stage were significantly associated with risk score ( ).

Development of Nomograms to Predict Individual Survival Outcomes

We developed nomograms based on the screening cohort and further verified their predictive ability in the validation cohort. It showed that age, T stage, N stage, M stage, and risk score are significant prognostic factors ( ). In the first step Cox proportional hazards analysis, we incorporated age, T stage, N stage, and M stage. It showed that age, T stage, and N stage were independent prognostic factors ( ) and were used to construct nomogram 1 ( ). In the second step Cox proportional hazards analysis, we incorporated age, T stage, N stage, M stage and risk score. It showed that age, T stage, N stage and risk score were independent prognostic factors ( ) and were used to construct nomogram 2 ( ).

Figure 8

Construction of nomograms. (A) Univariate analysis included risk score, age, gender, grade, M stage, T stage and N stage in the screening cohort. (B) Cox proportional hazards analysis included age, M stage, T stage and N stage in the screening cohort. (C) Cox proportional hazards analysis included risk score, age, M stage, T stage and N stage in the screening cohort. (D) Nomogram 1 based on the clinicopathological characteristics. (E) Nomogram 2 based on the risk score and clinicopathological characteristics.

Comparison of Prognostic Performance and Clinical Usefulness Between Nomogram 1 and Nomogram 2

In the screening cohort, nomogram 2 showed superior prognostic ability [AUC 0.684, 95% confidence interval (CI), 0.630–0.735] compared with nomogram 1 (AUC 0.639, 95% CI, 0.584–0.692) ( ). The calibration curves of nomogram 2 at 3 years also showed better consistency between the predicted and observed survivals than that of nomogram 1 ( ). Nomogram 2 showed higher net benefit than nomogram 1 between the threshold probabilities of around 37–60% in predicting 3-year overall survival ( ). Similar results were found in the independent validation cohort ( ).

Figure 9

The areas under the curve (AUC), calibration curve and decision curve analysis (DCA) for predicting patient survival. (A, D) The AUCs assess the accuracy of the nomograms in the screening and validation cohorts. (B, E) The calibration curves assess the consistency of the nomograms in the screening and validation cohorts. (C, F) DCAs assess the clinical usefulness of nomograms in the screening and validation cohorts.

Expression Levels of SERPINE1 and EFNA3 in GC Cell Lines and Tissues

In the screening cohort, the expression of SERPINE1 and EFNA3 in tumor tissues was up-regulated when compared with adjacent non-cancerous tissues and normal tissues (P < 0.05, ). To confirm the expression levels of SERPINE1 and EFNA3 in gastric cancer, we subsequently verified it in gastric cancer cell lines and patient tissues by qRT-PCR experiments. The results showed that, when compared with gastric normal epithelial mucosae cell line GES-1 and adjacent non-cancerous tissues, the expression of SERPINE1 and EFNA3 were significantly higher in gastric cancer cell lines ( , except for SERPINE1 in HGC-27, P > 0.05) and gastric cancer tissues (P < 0.05) ( ).

Figure 10

EFNA3 and SERPINE1 are upregulated in gastric cancer cell lines and tissues. (A, B) Bioinformatics analysis of the expression of EFNA3 and SERPINE1 in cancer and non-cancerous tissues in TCGA. (C, D) Bioinformatics analysis of the expression of EFNA3 and SERPINE1 in 27 pairs of gastric cancer and adjacent non-cancerous tissues in TCGA. (E, F) qRT-PCR results of EFNA3 and SERPINE1 expression level in GES-1 and gastric cancer cell lines. (Data are presented as mean ± SD. NS: P ≥ 0.05, *P < 0.05, **P < 0.01, ***P < 0.001). (G, H) qRT-PCR results of EFNA3 and SERPINE1 expression level in 39 pairs of gastric cancer and adjacent non-cancerous tissues. (Data are shown as –ΔΔCT values).

Discussion

Hypoxia is caused by an imbalance between insufficient oxygen supply and increased oxygen demand (21, 33). It is also one significant characteristic of tumor microenvironment. Tumor cells adapt to and rely on tumor microenvironment, contributing to instability and diversity of gene mutations, and activating a variety of signaling pathways and cytokines, contributing to the angiogenesis, invasion, metastasis, epithelial-mesenchymal transition, cancer stem cell maintenance, immune escape and resistance to radiotherapy and chemotherapy (34, 35). Therefore, understanding the molecular mechanism of hypoxia is critical to improving the survival of cancer therapy.

In this study, we identified two prognosis-related hypoxia genes, SERPINE1 and EFNA3, and establish a hypoxia risk score model based on the two genes. Subsequent survival analysis indicated that the high-risk group was associated with poorer prognosis, which was verified by an independent GEO cohort. GSEA analysis showed that the high-risk group was significantly enriched in pathways for tumor progression, such as the JAK-STAT signaling pathway (36), cancer in pathway, TGF-β signaling pathway (37, 38), and NOTCH signaling pathway (39), leading to poor prognosis. In the hypoxic microenvironment, HIFs are the main regulators of hypoxic response (18, 35). HIFs can cause the malignant phenotype of tumors by activating or enhancing JAK-STAT signaling pathway, TGF-β signaling pathway and NOTCH signaling pathway (40–43). Besides, the TCGA data base and qRT-PCR analysis confirmed the overexpression of these two hypoxia genes in tumor tissues and gastric cancer cell lines when compared with normal tissues and gastric normal epithelial cell line. Finally, the risk score model, age, T stage, and N stage were identified as independent risk factors related to OS and included in the nomogram. It showed that the nomogram was an effective tool for predicting the prognosis. The two-gene signature has a powerful ability to predict the prognosis of patients with gastric cancer, and may be helpful to guide clinical treatment decisions.

Tumor purity can reflect the characteristics of the tumor microenvironment. The risk score was significantly positively correlated with infiltrating immune cells and stromal cells, but negatively correlated with tumor purity. Previous studies showed that low tumor purity is associated with poor prognosis of multiple tumor types (44–46). We speculated that low-purity tumors may recruit more tumor immunosuppressive cells than high-purity tumors, and further studied the relationship between the risk score and the subtypes of infiltrating immune cells. We found that the tumors in the high-risk group contained more infiltrating immunosuppressive cells such as Tregs, macrophages, neutrophils, para-inflammatory and mast cells than the low-risk group. A previous study found that Tregs suppressed the anti-tumor immune response by weakening the cell-mediated immune response to tumors, thereby promoting disease progression (47). Hypoxia can protect tumors from the intrinsic anti-tumor immune response by forming an immunosuppressive microenvironment, which may explain a poor prognosis of the high-risk group.

Cytokine are important factors in regulating tumor immunity. Among them, tumor immunosuppressive cytokines are important factors inhibiting immune cell activity. Transforming growth factor-β (TGF-β) suppresses the immune system by inhibiting the maturation of dendritic cells, inhibiting the activity of NK cells, and reducing the cytotoxicity of T cells (47, 48). Interleukin 10 (IL-10) is an immunosuppressive cytokine secreted by T-helper 2 (Th2) cells, Tregs, and M2 macrophages. It has been shown to impair the proliferation, cytokine production and migration capabilities of effector T cells (49). IL-10 also promotes the stable expression of Foxp3, TGF-β-receptor 2 and TGF-β, thereby stabilizing the phenotype and functions of Treg (50). In our research, the immunosuppressive cytokines, such as IL-10 and TGF-β, were up-regulated in the high-risk group, thereby further promoting immunosuppression.

The correlation between the intrinsic escape mechanism and risk score is clinically important. The inherent immune escape of tumors demonstrates that tumor cells can mediate their own immune escape directly. Previous study has illustrated that the expression of immune check-point molecules and tumor immunogenicity are two important aspects of intrinsic immune escape (51). Immune checkpoint molecules play a key role in tumor progression and carcinogenesis by promoting tumor immunosuppression. Malignant tumors can evade immune killing by stimulating immune checkpoint target genes (such as PD-1, PD-L1, CTLA-4, TGF-β, and HAVCR2). In this study, immune checkpoint molecules of PD-1, TGF-β, and HAVCR2 were up-regulated in the high-risk group. This result indicates that tumor cells in the high-risk group express immune checkpoint molecules to protect themselves from attack.

Another potentially significant intrinsic immune escape mechanism is immunogenicity. Some somatic mutations in tumor DNA produce neoantigens, and the antigens from this mutation are recognized and targeted by the immune system, especially after treatment with drugs that activate T cells (52–56). The more somatic mutations are present in a tumor, the more neoantigens it may form. TMB may represent a better estimate of tumor neoantigen burden (57). Here, we found that the high-risk group had a lower proportion of somatic mutations, and the hypoxia risk score was significantly negatively correlated with TMB. Tumor cells in the hyperoxia group produced fewer neoantigens, thus avoiding being recognized and killed by T cells.

To investigate the role of hypoxia risk in drug treatment, our research showed that, the tumors in the high-risk group were not sensitive to most chemotherapy drugs, such as axitinib, bexarotene, bortezomib, and imatinib. However, the tumors in the high-risk group were more sensitive to methotrexate and mitomycin.c and may benefit from these two chemotherapy drugs. Those in the high-risk group express higher levels of PD-1, HAVCR2 and other immune checkpoint molecules to avoid the attack of anti-tumor immune cells. The high-risk group may benefit from immunotherapy, such as the use of PD-1 and HAVCR2 inhibitors.

Nomograms are commonly used to assess the prognosis of tumors (58, 59). In this study, we constructed two prognostic nomograms. Nomogram 1 is based on clinical characteristics, and nomogram 2 is developed by the combination of clinical characteristics and the hypoxia risk score model. It showed that the prognostic nomogram based on the combination of clinical characteristics and hypoxia risk score model has better predictive ability and higher clinical usefulness. However, owing to the lack of in vitro or in vivo experiments, the reliability of our molecular mechanism analysis may be limited.

Conclusions

In summary, we developed and validated a hypoxia risk score model based on a novel hypoxia-related gene signature revealing the relationship between hypoxia and tumor immune microenvironment. The current study may provide new insights into how hypoxia affects the prognosis, and may be helpful in guiding targeted hypoxia therapy for gastric cancer.

Data Availability Statement

The datasets analyzed during this current study are available in TCGA (http://cancergenome.nih.gov/) and GEO (https://www.ncbi.nlm.nih.gov/geo, GSE84437) data bases. The data are also available from the corresponding author (zhangchundong2007@126.com or daidq63@163.com) on reasonable request.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics Committee of the Medical Association of the Fourth Affiliated Hospital of China Medical University (EC-2021-KS-043). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

C-DZ and D-QD have contributed equally to this work as co-corresponding authors. J-PP and C-DZ wrote the main text and performed data analysis. J-PP, C-DZ, and D-QD designed the study and collected the data. All authors contributed to the article and approved the submitted version.

Funding

This research was funded in part by the China Scholarship Council (201908050148).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.