Identification of prognostic hypoxia-related genes signature on the tumor microenvironment in esophageal cancer

Abbreviations: EC: Esophageal cancer; EAC: Esophageal adenocarcinoma; ESCC: Esophageal squamous cell carcinoma; PD-L1: Programmed death-ligand 1; TME: Tumor microenvironment; TCGA: The cancer genome atlas; PPI: Protein-protein interaction network analysis; TMN: Tumour node metastases; HR: Hazard ratio; GSEA: Gene set enrichment analysis; ROC: Receiver operating characteristic; AUC: Area under curves; HIF: Hypoxia-inducible factor; TAMs: Tumor-associated macrophages; SLAM: Signaling lymphocyte activation molecule; CD274: Cluster of differentiation 274; B7-H1: B7 homolog 1

1.   Introduction

Esophageal cancer (EC) is the eighth most common cancer in the world. Non-industrialized countries account for nearly four out of every five instances, with Asia and Africa holds the greatest rates ^[1].Esophageal adenocarcinoma(EAC) and esophageal squamous cell carcinoma(ESCC) account for over 95% of esophageal cancers ^[2]. Esophageal cancer is related to a poor prognosis. Despite advances in diagnosis and treatment, esophageal cancer patients' five-year overall survival rate is 15 to 20 percent globally and in the United States ^[3].

Hypoxia is an essential component of the tumor microenvironment and one of the most common features of solid tumors ^[3]. Tumor cells are more vulnerable to proliferation, metastasis, drug resistance, and poor prognosis when exposed to hypoxia ^[4,5].

Hypoxia stress can cause tumor resistance, resulting in immunosuppression ^[6,7]. Several studies have shown that the expression of programmed death-ligand 1 (PD-L1) in tumors is up-regulated under hypoxia ^[8]. More and more studies began to focus on analyzing the tumor microenvironment (TME) in recent years. Immune cells in tumor microenvironment play an essential role in tumorigenesis ^[9,10,11]. These related immune cells may have anti-tumor or tumor-promoting functions, and their specific mechanism is very complex ^[12].

Moreover, there are several (non-)aggressive hypoxia biomarkers that can be identified in esophageal cancer. In esophageal cancer, the molecular responses associated with hypoxia are an important clinical phenomenon, and elevated levels of hypoxia-related biomarkers are often associated with poorer treatment outcomes such as overall survival, disease-free survival, complete remission, and local control ^[13]. Various myeloid cells make up the immune portion of the TME. Hypoxia can disrupt lymphocyte differentiation and function by altering the expression of co-stimulatory receptors and the cytokines released by these cells ^[14].

The high incidence and s severe adverse complications of gastrointestinal tumors seriously affect the long-term survival rate and quality of Life In the patient population. It is well known that early detection of the disease will lead to longer survival rates and better quality of life ^[15]. At present, a variety of molecules have been reported to be used as prognostic markers and targets for the development of new therapies, such as Autophagy-related microRNAs ^[16], Circular RNAs ^[15], curcumin ^[17,18]

Multiple models have been developed to forecast the prognosis of esophageal cancer. However, there is no risk model based on hypoxia-related genes to predict the prognosis of esophageal cancer. Therefore, we generated a hypoxia-related gene-based signal in this study to look into the possible value of hypoxia-related genes in the prognosis of esophageal cancer patients.

2.   Materials and methods

2.1.   Data collection

GEO publicly provides genomic data to support MIAME-compliant data submissions ^[14]. Human lncRNA, miRNA and mRNA expression were freely accessible from NCBI-GEO (GSE124401 and GSE136547). To be specific, the threshold of up-regulated or down-regulated lncRNA/mRNA was |log2 (fold-change) | > 1 with adjust P < 0.05.

The RNA-seq transcriptome data and corresponding clinicopathological information from The Cancer Genome Atlas (TCGA) are obtained for 82 ESCC patients and 78 EAC patients (http://cancergenome.nih.gov/). The MSigDB database provided the hypoxia-associated gene set (Hallmark-hypoxia) (https://www.gsea-msigdb.org/gsea/index.jsp). Hypoxia-related genes have been shown to be upregulated in the lab in reaction to low oxygen levels. In this study, we used an illustration package in R to map the heatmaps of gene expression.

The data in the TCGA database is freely accessible to the public and this study follows a strict access policy and publication guidelines, thereby ethical review and approval by the ethics committee are not required for this study.

2.2.   Protein-protein interaction (PPI) network analysis

The Search Tool to Retrieve Interacting Genes (STRING) database (http://string-db.org/cgi/input.pl) is used to construct the hypoxia gene network. The 50 genes with the most adjacent nodes are chosen for further investigation using R software.

2.3.   The hypoxia model's construction and grouping

Testing for oxygen deprivation genes and discovering genes that are separately associated with esophageal cancer prognosis. In a multivariate Cox regression, correlation coefficients are obtained after finding statistically relevant hypoxic genes in a univariate Cox regression. Moreover, the risk score of each patient is determined by multiplying the expression levels of each gene in the TCGA databases by the expression coefficients for each gene. Based on median risk score values, patients in both databases are categorized into high-risk and low-risk categories.

The risk score equation is defined as follows: Risk score = ∑Ni = 1(Expi × Coei)

For example: In EAC, N = 5, Exp_iis the expression of each of the five hypoxic genes and Coe_i is the multivariate Cox regression coefficient that corresponds.

2.4.   Independence of prognostic genetic features from other clinical features.

We used R's "survival" kit and conducted univariate and multivariate Cox regression analysis to identify hypoxia-related genes linked to prognosis. Then, to determine whether hypoxia risk score is an independent factor affecting prognosis, age, gender, tumor node metastases (TNM) level and the proposing risk score are used. Statistical significance is identified as a bilateral P value of less than 0.05. The hazard ratio (HR) is estimated, as well as the 95 percent confidence intervals.

2.5.   Gene set enrichment analysis (GSEA)

During enrichment between high and low-risk groups in two groups of esophageal cancer, using GSEA (https://www.gsea-msigdb.org/gsea/index.jsp) to identify significant variations between the groups of genes expressed in the MSigDB pool (h.all.v7.0.cymbols.gmt; c5.bp.v7.0.symbols.gmt). For each analysis, 1000 permutations of gene sets are performed. A risk index is calculated using the phenotype mark.

2.6.   Distinguish of immune cell types fractions and gene expression

CIBERSORT ^[25] is used to calculate immune cell types in low-and high-risk classes. In each sample, the number of all predicted immune cell type scores equals 1. The online platform(http://biocc.hrbmu.edu.cn/TIP/index.jsp) for tracking tumor immunophenotypes to screen for immune-related genes that play an important role in the control of this immune cell type. In R serve, the ggExtra packages, ggplot2 and ggpubr to explore the relationship between the hypoxic risk and the gene expression, as well as the expression levels that differ between the high and low risk of hypoxia, in patients looking for genes that perform a critical role in the regulation of immune cells.

3.   Results

3.1.   Hypoxia-related genes extraction and screening

The gene expression data and clinical data of esophageal cancer and hypoxia related genes were downloaded from the TCGA database and GSEA website respectively. We later used STRING to model PPIs among hypoxia-related genes (Figure 1A, B). The number of nodes in the network for every protein is used to identify the central genes. Figures 1C, D showed the top 50 essential genes. To find genes linked to prognosis, we used univariate Cox regression analysis and multivariate Cox analysis.

PGM1 is correlated with low risk in the EAC, while two genes (PGK1 and SLC2A3) are associated with a high risk of malignancy development by using univariate Cox regression analysis (Figure 1E). Finally, multivariate Cox regression analysis identified three hub genes related to patient prognosis (Figure 1F). The model coefficients for these three genes (PGK1, PGM1 and SLC2A3) are -2.8966, 0.5736 and 1.6160, respectively.

Furthermore, univariate Cox and multivariate Cox showed that EGFR relates to a low risk, while ATF3 was correlated with a high risk of malignancy development in ESCC (Figure 1G, H). The model coefficients for these two genes (EGFR and TPF3) are -0.482 and 1.394, respectively.

3.2.   Association between hypoxia related genes and prognosis

We used the genes from the Cox regression study to model and used the median risk score to classify patients into high-risk and low-risk types. Subsequently, the survival rate was significantly increased in the low-risk categories in both ESCC and EAC groups (p < 0.05; Figure 2A, B).

The accuracy of this model in predicting survival is tested using receiver operating characteristic (ROC) curve analysis. The EAC three-gene prognostic model has AUCs of 0.728, 0.778 and 0.921 after one, three, and five years respectively (Figure 2C). The ESCC three-gene predictive model has Area Under Curves (AUCs) of 0.654, 0.871 and 0.982 after one, three, and five years respectively (Figure 2D). This means that the model gets more accurate over time in predicting survival rates for EAC and ESCC patients at 1, 3 and 5 years, implying that the model correctly predicted survival rates in EAC and ESCC.

Figures 3A, B showed the interaction between hub hypoxia-related genes. The risk curve showed that the survival time of the low-risk group was higher than that of the high-risk group Figure 3C, D). In addition, the high-risk category showed that the number of deaths increased over time (Figure 3E, F). Heat Maps showed the expression level of the hub gene in different groups. The expression of these genes was significantly different in high and low-risk groups (Figure 3G, H).

3.3.   Effects of different clinical characteristics on esophageal cancer 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 EAC and ESCC 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 TCGA-EAC and TCGA-ESCC database, we found that patient gender and tumour stage affected neither survival nor prognosis, while node stage affected EAC patients, metastases stage affected ESCC patients (Figure 4A, B). The p-value of the risk score in our model was < 0.05 for patients in both databases, indicating that the risk score affected patient prognosis and survival. Multivariate Cox regression analysis indicating that risk score about hypoxia were independent prognostic factors (Figure 4C, D).

3.4.   GSEA identifies hypoxia-related signaling pathways

We use GSEA to find hypoxia-related signaling pathways. When compared to the TCGA-EAC database's low-risk category, the high-risk group had several enriched pathways linked to apoptotic, hypoxia, and inflammatory response pathways (Supplemental Table 1). On the other hand, the pathways associated with metabolism are enriched in the low-risk community (Supplemental Table 2), and similar pathways are shown in the TCGA-ESCC database in metabolism (Supplemental Tables 3 and 4)

3.5.   Immune landscape between low and high hypoxia risk

According to GSEA, we investigated the ability to characterize hypoxia risk in the evaluation of the immune microenvironment. Using the CIBERSORT method combined with the LM22 feature matrix, we estimated the differences in immune infiltration of 22 immune cell types in patients with low- and high-risk esophageal cancer. The outcomes of immune infiltration for patients with EAC and ESCC were summarized in the Figure 5A, B.

In EAC patients, the proportion of M2 macrophages is significantly higher in high risk (Figure 5C). In ESCC patients, the proportion of activated mast cell and activated CD4 memory T cell is significantly higher in in high risk (Figure 5D, E). Visualization of the expression of immune-related genes from the online platform mentioned in the method were by Heatmap. Several immune-related genes differed significantly between the high- and low-risk groups (p < 0.05) (Figure 6A, D).

In two types of esophageal cancer, CD244 and CD274 are found to be up regulated in high-risk categories and their expression levels are in positive correlation with the risk scores of the patients, according to the correlation curves between their risk scores and their expression levels (Figure 6B, C, E, F).

4.   Discussion

The poor prognosis and increasing occurrence of esophageal cancer illustrate the importance of better diagnosis and prediction methods before treatment ^[19]. Numerous reports were linking hypoxia-related genes to esophageal cancer ^{[20,21,22,23]}. In addition, for patients with esophageal cancer, mRNA gene markers dependent on certain tumor microenvironment features have high predictive ability.

Hypoxia may cause multiple transcription factors to react by inducing or repressing genes, resulting in the initiation of an adaptive transcriptional response ^[24]. Hypoxia-inducible factor (HIF) is the most significant of these, dubbed a "chief regulator" of the "metazoan adaptive response to hypoxia" ^[25,26,27].

To learn more about the connection between hypoxia-related gene expression and patient prognosis, the researcher constructed a hypoxia analysis model and divided the hypoxia analysis model into two risk groups. Besides, the difference in prognosis between patients in the high-risk and low-risk groups is statistically significant. In both ESCC and EAC, the survival rate of patients in the low-risk group is substantially higher than that of the control group.

We looked at the influence of other variables on the prognosis of patients, which included age, gender, TNM stage, and our proposed risk score on the prognosis of the patients. We conclude that our risk score is an independent predictor of esophageal cancer and is valid in both esophageal squamous cell carcinoma and esophageal adenocarcinoma.

Next, the enriched pathways in the high-risk and low-risk groups are identified using GSEA. Most of the enriched pathways in esophageal cancer is linked to angiogenesis, apoptosis, hypoxia, the inflammatory response, and metabolic pathways, according to our findings. In tumor tissue, hypoxia can lead to the release of many cell debris and inflammatory factors, macrophages and monocytes are attracted and cause macrophage polarization. After polarization, macrophages secrete inflammatory factors ^[28]. These results suggest that hypoxia, inflammation, and the immune system's response are all related. Therefore, we studied the differences in immune function between the high and low risk groups in the hypoxia model.

We discovered that patients with a high risk of hypoxia had a substantially higher proportion of M2 Macrophage cells. Macrophages are a significant component of the immune infiltration of solid tumors, as they differentiate into tumor-associated macrophages (TAMs), which are preferentially present in hypoxic zones of the tumor ^[29]. Macrophages can be stimulated by a multitude of stimuli and polarized to functionally distinct phenotypes, including classically activated (M1) and alternatively activated (M2) macrophages ^[30]. Infiltrating macrophages are classified as the M2 phenotype in most tumors, and they provide an immunosuppressive microenvironment for tumor development. In addition, many cytokines, chemokines, and proteases are secreted by tumor-associated macrophages, which promote tumor angiogenesis, development, metastasis, and immunosuppression ^[31,32]. Cancer patients typically receive intensive surgery of the tumor tissue, as well as radiochemical or hormone therapy. However, none of these drugs are unique for tumors, and they all have significant side effects. As a result, using macrophages to mediate tumor resistance is being considered as a possible therapy. The discovery of mechanisms underlying macrophage functional polarization into M1 or M2 cells could lead to new insights into developing macrophage-centered diagnostic and therapeutic strategies for a variety of diseases.

We also looked for genes that could influence immune cell infiltration rates and looked at the relationship between the risk of hypoxia and these genes. Though EAC and ESCC have been treated as distinct entities due to their widely differing and overlapping etiologies, and targeted treatments for each have been investigated. Interestingly, regardless of the type of esophageal cancer, PD-L1 and CD244 levels increased in the high-risk population.CD244 belongs to the Signaling Lymphocyte Activation Molecule (SLAM)class of Ig Superfamily receptors. It is a transmembrane receptor made up of an extracellular segment with two immunoglobulin (Ig) -like domains, a transmembrane region, and a cytoplasmic domain with tyrosine-based motifs. The rest of the SLAM family. Unlike other SLAM family receptors, it does not function as a self-ligand; instead, it binds CD48, a transmembrane receptor widely expressed in hematopoietic cells ^[34,35]. In a recent review on squamous cell carcinoma of the head and neck, researchers evaluated the therapeutic potential of CD244 and found that CD244 plays a role in the overall immunosuppressive environment and therefore has the potential to be a new immunotherapeutic target for malignancies ^[33].

Another gene was Cluster of differentiation 274 (CD274), also known as B7 homolog 1 (B7-H1) or PD-L1, which was the first functionally characterized ligand of the co-inhibitory PD-1 ^[34]. PD-L1 is an immunosuppressive molecule that can inhibit the activation of T cells, leading to tumor progression ^[35]. A large number of studies have confirmed that the expression of PD-L1 on tumor cells is related to the poor prognosis of tumor patients, such as lung cancer ^[36,37], gastric cancer ^[38], and colorectal cancer ^[39]. Two studies have reported a positive correlation between PD-L1 expression and poor Esophageal cancer prognosis ^[40,41]. In addition, numerous studies suggest a crucial role for PD-L1 in cancer immune escape ^[42,43,44]. Some studies reported that blocking PD-L1 could enhance anti-tumor immunity ^{[45,46,47,48]}. PD-L1 upregulation occurs in approximately 40% of esophageal cancers. Preliminary clinical data show remission rates of 22–27% in metastatic gastroesophageal cancer treated with PD-1/PD-L1 inhibitors and 10–17% in unselected patients ^[49].

Our research provides a more theoretical basis for PD-L1 to become a targeted treatment for esophageal cancer. In the future, we will pay more attention to the therapeutic effect of PD-L1 in esophageal cancer. However, due to the low accuracy of PD-L1 immunohistochemical staining, PD-L1 as the only predictive biomarker for cancer immunotherapy is problematic ^[50]. Therefore, it is necessary to understand the tumor microenvironment better and use other biomarkers, for example, in combination with other biomarkers to predict the prognosis of esophageal cancer, better to identify the benefits of PD-1/PD-L1 checkpoint blockade Therapy patients.

Our result was the first hypoxia prediction model based on esophageal cancer, and at the same time identified crucial genes linked to the immune microenvironment in esophageal cancer, pointing to possible molecular targets for esophageal cancer immunotherapy, It will enrich the research on the hypoxia and immunity therapy in esophageal cancer, which deserves in-depth exploration by follow-up researchers.

However, there are some limitations to our study. First, the development and progression of esophageal cancer are influenced by many other complex mechanisms, and the use of a single feature to construct a predictive model may be inherently flawed. Secondly, our evidence is based on a bioinformatics methodology that needs to be validated by more cell experiment and clinical trials before it can be considered initially. Furthermore, our inability to find molecularly targeted drugs in a variety of situations, including normoxia and hypoxia (acute and chronic).

Acknowledgements

Thanks for the support of "Yinzhou district science and technology project in the field of agricultural and social development" with NO.20201YZQ010102 and NO.20211YZQ010088. The datasets generated during the current study are available in the TCGA database (https://portal.gdc.cancer.gov/).

Conflict of interest

The authors declare that there are no conflicts of interest related to this paper.