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)(7)(8)(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)(11)(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)(17)(18). Multiple genes transcribed by HIFs, including Glut1, KLF8, VEGFA, ITGb1, 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 m atrix 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.

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 [log 2 (x+1)] and normalized by the "sva" and "limma" packages of R software. The baseline characteristics of the screening and validation cohorts are shown (Supplementary Table S2).

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 = (∑coefficient x * expression of signature gene x ) (gene x 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.
(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% CO 2 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 Table S1.

Statistical Analyses
All analyses were performed by R version 4.0.2 (http://www.Rproject.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.

Patient Characteristics
The baseline characteristics of the screening and validation cohorts are shown in

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 (Table S3 and Figures 1A-C). To better visualize the interactions between these hypoxia genes, the STRING online data base was used to analyze the protein-protein interaction network ( Figure 1D). 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, Table S4). 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 (l) ( Figures 1E, F). Finally, the Cox proportional hazards analysis confirmed that two genes (SERPINE1, EFNA3) ( Figure 1G) 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.

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 ( Figure 2A). The mortality rate in the low-risk group was significantly lower than that in the highrisk group ( Figures 2B, C). 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) ( Figure 2D). Similar results were found in the validation cohort ( Figures 2E-H).

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-b signaling pathway were activated ( Figure 3A). 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 3B).

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 ( Figures S1A-H). 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, Figures 3C-F).

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, Figures 4A-H). 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.
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 (Figures S2A, B). 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, Figure S2A). However, the prognostic protective gene EFNA3 showed the opposite trend for most aspects except B cells ( Figure S2B, 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 ( Figures 5A, B). In the screening cohort, the expression level of five key immune checkpoint molecules (PD-1, PD-L1, CTLA-4, HAVCR2 and TGF-b) in the high-risk group was significantly higher than those in the low-risk group, and significantly positively correlated with risk score (Figures 5C-G). Similar results were obtained in the validation cohort, except that CTLA4 and PD-L1 showed no significant difference between the high-risk and lowrisk groups (Figures 5H-L).

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; Figure 6A). 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) ( Figure 6B). 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, Figure 6C). 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; Figure 6D). 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. 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 (Figures 6E, F). 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; Table S5).

The Correlation Between Risk Score and Clinicopathological Characteristics
We conducted correlation analyses between clinicopathological factors and risk score in the screening cohort ( Figures S3C), and tumor grade and T stage were significantly associated with risk score (Figures S3C, D). In the validation cohort, age, T stage, and N stage were significantly associated with risk score ( Figures  S3H, J, K).

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 ( Figure 8A). 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 ( Figure 8B) and were used to construct nomogram 1 ( Figure 8D). 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 ( Figure 8C) and were used to construct nomogram 2 ( Figure 8E).

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, Figures 10A-D). 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 ( Figures 10E, F, except for SERPINE1 in HGC-27, P > 0.05) and gastric cancer tissues (P < 0.05) (Figures 10G, H).

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-b 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-b signaling pathway and NOTCH signaling pathway (40)(41)(42)(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)(45)(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-b (TGF-b) 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-b-receptor 2 and TGF-b, thereby stabilizing the phenotype and functions of Treg (50). In our research, the immunosuppressive cytokines, such as IL-10 and TGF-b, 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-b, and HAVCR2). In this study, immune checkpoint molecules of PD-1, TGF-b, 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)(53)(54)(55)(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.

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 cocorresponding 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).