{'fraction_non_alphanumeric': 0.03493112607325962, 'fraction_numerical': 0.017295694607449504, 'mean_word_length': 5.3851311378426345, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Construction and Verification of a Hypoxia-Stemness-Based Gene Signature for Risk Stratification in Esophageal Cancer

Background Numerous studies have shown that esophageal cancer (ESCA) contains areas of intertumoral hypoxia. It is widely accepted that the association of hypoxia with cancer stemness in the tumor microenvironment of ESCA is of profound clinical significance. However, reliable prognostic signatures based on hypoxia and cancer stemness are still lacking in ESCA. Material/Methods The t-SNE algorithm was used to estimate the hypoxia status based on the transcriptome profiles of the discovery cohort in the TCGA database. Median values of the stemness index were used to group and identify stemness-associated differentially expressed genes (DEGs). The LASSO method and Cox regression model were combined to screen for prognostic genes and to establish a genetic signature based on hypoxia-stemness. The robustness of the prognostic model was then tested in an external independent validation cohort of the GEO database. Results A total of 8 genes – FBLN2, IL17RB, CYP2W1, AMTN, FABP1, FOXA2, GAS1, and CTSF – were identified to construct a gene signature for ESCA risk stratification. Overall survival was significantly lower in the high-risk group than in the low-risk group in both the internal discovery set and the external validation set. The risk score was found to be an independent prognostic factor for ESCA patients. In addition, a higher risk score was significantly associated with the sensitivity of ESCA patients to gefitinib, bexarotene, dasatinib, and imatinib. Conclusions The hypoxia-stemness-based genetic signature established for the first time in our study could be a promising tool for ESCA cancer risk stratification.


Background
Esophageal cancer (ESCA) is one of the most commonly diagnosed cancer types and is the sixth most lethal cancer worldwide. With an increasing annual incidence and an overall 5-year survival rate of approximately 15-20%, the prognosis for patients with esophageal cancer remains poor [1]. To improve the prognosis of patients with esophageal cancer, multimodal combination therapy has been incorporated into clinical applications in recent decades. However, surgery remains the mainstay of treatment for patients with esophageal cancer. The burgeoning of adjuvant treatment has played a central role in improving survival outcomes, especially in low-stage esophageal cancer patients. Neoadjuvant chemotherapy has been shown to be superior to surgery alone, with a yield of 12-15%, thus becoming an essential component of current treatment for patients with curable esophageal cancer [2]. Unfortunately, despite advances in early screening, radical surgical resection, and recent advances in immunotherapy, esophageal cancer remains a silent disease usually diagnosed in the middle to late stages, leading to its high mortality rate. Therefore, there is an urgent need to explore the characteristics of ESCA to develop new therapeutic approaches.
As a hallmark of cancer, hypoxia is also associated with a higher capacity for metastatic progression, tumor recurrence, immune evasion, and resistance to radio-chemotherapy [3]. Since tumor hypoxia cannot be predicted based on tumor size, TNM stage, or pathological grade, there is an urgent need for validated molecular biomarkers capable of assessing ESCA hypoxia status. To the best of our knowledge, only HIF-1a and carbonic anhydrase IX (CAIX) are valid molecular markers capable of evaluating the hypoxic response in ESCA [4,5]. An early study demonstrated a significant relationship between HIF-1a expression and VEGF expression, venous infiltration, and, ultimately, clinical prognosis in ESCA [6]. During cancer development and progression, tumor cells, including CSCs, are often in a hypoxic state, and this hypoxic microenvironment maintains CSCs to acquire a more aggressive phenotype [7].
As the name implies, cancer stem cells (CSCs) are a subpopulation of cancer cells with characteristics similar to those of normal stem cells. The current widely accepted Cancer Stem Cell theory claims that a single stem cell is sufficient to initiate and generate the remaining tumor cells. The heterogeneity of specific tumor cell subpopulations with different potential for self-renewal and differentiation has led to a significant focus on the concept of cancer stem cells (CSCs) [8]. It is known that hypoxia regulates CSCs stemness by targeting the Wnt/b-catenin, Notch, Hedgehog, PI3K/AKT, and STAT3 pathways [9]. Studies have shown that esophageal cancer cells exposed to hypoxia are arrested in the G1/S phase. Since cells in the G2/M phase are most sensitive to irradiation, hypoxia-induced arrest in the G0-G1 phase could be considered as a potential mechanism leading to radiation resistance in esophageal cancer cells [10].
The present work hypothesized that the interaction between hypoxia and stemness has some prognostic significance in ESCA. Through a series of systematic bioinformatics analyses, we developed a novel genetic signature by integrating stemness and hypoxic states followed by incorporation into the current clinicopathological signature and staging system, aimed at improving the prognosis of patients with esophageal cancer and supporting personalized treatment.

Datasets
As the largest publicly available database of cancer genetic information, TCGA (https://portal.gdc.cancer.gov/) covers data including gene expression data, miRNA expression data, copy number variants (CNVs), DNA methylation, and singlenucleotide polymorphisms (SNPs). We downloaded processed esophageal cancer mRNA expression data (level 3 data) from the TCGA esophageal cancer RNA-seq dataset, including a normal group (n=11) and a tumor group (n=161). Similarly, we also downloaded the series matrix file data file of GSE53624 (n=119) from the NCBI GEO public database (http://www.ncbi. nlm.nih.gov/geo/), annotated with the platform GPL18109, and extracted the complete gene expression profiles and clinical survival information of 119 ESCA patients.

Esophageal Cancer Subtype Classification and Differentially Expressed Genes (DEGs) Analysis
A total of 480 hallmark hypoxia-related gene sets were obtained through the GeneCards database (https://www.genecards.org/) (Relevance score>5). Unsupervised nonnegative matrix factorization (NMF) [11] clustering was performed on the expression profile of hypoxia-related genes through the NMF package. We next assessed the correlation of all candidate genes with overall survival (OS) by using the "survival" package in R packages. The best cluster value is defined as the value of k at which the correlation coefficient starts to decrease. Then, based on the t-distributed Stochastic Neighbor Embedding (t-SNE) method [12], we further used the mRNA expression data of the aforementioned hypoxia genes to verify subtype assignment. Finally, we analyzed the differentially expressed genes (DEGs) between esophageal cancer subtypes through the Limma package [13]. The threshold for differential gene screening was: |logFC| >0.58 and P<0.05. In addition, TCGA esophageal cancer samples were further grouped according to the median of stemness index, and stemness-related differential genes were obtained by differential expression e934359-2 analysis between the 2 groups through the Limma package. The intersection of differential genes related to hypoxia and differential genes related to stemness index yielded the differential gene of the hub (candidate DEGs).

Functional Annotation and PPI Network Construction of Candidate DEGs
To fully explore the functional relevance of candidate DEGs, the R package cluster profile (R4.0) [14] was used for functional annotation of candidate DEGs. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis was further performed to evaluate related functional categories. The GO and KEGG enrichment terms with p and FDR q values less than 0.05 were considered significant. The STRING online database 11.0 (https://string-db.org) was used to explore the protein-protein interaction (PPI) network of candidate DEGs, and a confidence score >0.4 was used as the cut-off criterion. The generated PPI network is further visualized by Cytoscape software to explore the regulatory correlation between candidate DEGs.

Constitution of Hypoxia-Stemness Index Risk Model Based on Hub DEGs
We collect clinical information of ESCA patients based on the candidate DEGs, Cox single factor regression, and LASSO regression algorithms [15] were used in combination to screen out essential genes (Hub DEGs) in esophageal cancer to further construct prognosis-related models. Specifically, risk scores were obtained for each patient with esophageal cancer after incorporating values for each hub gene at the transcriptional level, the estimated regression coefficients derived from LASSO regression analysis were used for weighting. According to the risk score formula, the median risk score value was defined as the cut-off point, and esophageal cancer patients were further stratified into low-and high-risk subgroups. The probability of survival curves was estimated by the Kaplan-Meier estimator and any statistical difference in survival between groups was assessed by the log-rank test. The predictive performance of the models was finally evaluated by the ROC curve.

Drug Susceptibility Tests
Based on the largest publicly available pharmacogenomics database, the Genomics of Drug Sensitivity in Cancer (GDSC, www. cancerrxgene.org), we initially utilized the R software package "pRRophetic" to predicted the chemotherapy response of each esophageal cancer patient. The regression method was used to obtain the IC 50 estimated value for each specific chemotherapy drug treatment, and we performed the tenfold cross-validation approach on the GDSC training set. All parameters were run at default values, and batch effects were corrected using ComBat.

Tumor Immune Infiltration Analysis
We combine the xCell algorithm [16] and the CIBERSORT algorithm [17] to evaluate the abundance of tumor-infiltrating immune cells.The correlation between risk score and immune cell abundance was analyzed by the Spearman correlation test. P<0.05 was considered as the threshold value of statistical difference.

Gene Set Variation Analysis (GSVA)
GSVA was used to transform the gene expression data from expression matrices with individual genes as features to expression matrices with specific sets of genes as features and quantified the gene enrichment results [18]. The multi-category gene set was downloaded from the molecular signatures database (MSigDB v7.0) (http://software.broadinstitute.org/ gsea/msigdb/), and the GSVA algorithm converted gene-level quantification into pathway-level quantification, thereby identifying differentially expressed pathways between samples.

Statistical Analysis
The Kaplan-Meier method and log-rank tests were used for survival analysis. Multivariate analysis was performed using the Cox proportional hazard model. All statistical analyses were performed in R (version 4.0). All statistical tests were two-sided, and P values <0.05 were regarded as statistically significant.

Hypoxia-Stemness-Related DEGs in ESCA Cohort
The baseline clinicopathological characteristics of the TCGA-ESCA cohort and GSE53624 dataset from GEO are shown in Table 1. We download the processed ESCA original mRNA expression data (FPKM) in the TCGA database. A total of 480 hypoxia-related genes were obtained through the GeneCards database (https://www.genecards.org/) (Relevance score >5), and we extracted a total of 469 hypoxia-related regulatory factor gene sets from the TCGA data. The NMF consensus clustering method was used to cluster the TCGA data set containing esophageal cancer samples based on the expression profiles of the above 469 candidate genes, and k=2 was selected as the best clustering value ( Figure 1A). When k=2, the dimensionality of the feature was reduced by the t-SNE method, and it was found that the subtypes were mainly consistent with the two-dimensional t-SNE distribution pattern ( Figure 1B). In addition, we observed a significant prognostic difference in the TCGA-ESCA data set. Compared with the C2 cluster, cluster C1 shown a longer median survival time (P=0.036) ( Figure 1C). Finally, a total of 1308 hypoxia-related DEGs were obtained through the difference analysis between the 2 subtypes through the Limma package in R (Figure 1D, 1F).

e934359-3
A total of 1071 DEGs related to stemness index were obtained by analyzing the difference between the 2 groups through the Limma package (Figure 1E, 1G). We used Venn diagrams to further screen the genes that are differentially expressed in both hypoxia subtypes and stemness index to obtain 282 genes as hypoxia-stemness-related DEGs for subsequent analysis ( Figure 1H).

Functional Enrichment of Hypoxia-Stemness-Related DEGs and Protein-Protein Interaction Network Construction
GO and KEGG enrichment analysis was used to annotate the functions of these hypoxia-stemness-related DEGs. As shown in Figure 2A, GO enrichment showed significantly enriched GO biological processes (BP) terms such as extracellular matrix organization, extracellular structure organization, connective tissue development, and cartilage development. The significant GO-cellular component (CC) terms included collagen-containing extracellular matrix and cell-cell junction. For molecular function (MF) analysis, the significant terms included extracellular matrix structural constituent, growth factor binding, and cell adhesion molecule binding. KEGG pathway analysis also identified numerous signaling pathways that were significantly enriched in the candidate DEGs, for example, the AGE-RAGE signaling pathway and PI3K-Akt signaling pathway ( Figure 2B). To analyze the role of candidate DEGs in the occurrence and development of esophageal cancer, we explored the regulatory correlation between candidate DEGs and established a PPI network based on the candidate DEGs. As shown in Figure 2C, according to the PPI network, the top 10 candidate DEGs with higher degrees were screened, including Fibronectin 1 (FN1), Collagen type I alpha 1 chain (COL1A1), Matrix metallopeptidase 2 (MMP2), Collagen type III alpha 1 chain (COL3A1), Collagen type I alpha 2 chain (COL1A2), Secreted protein acidic and cysteine rich (SPARC), Periostin (POSTN), Biglycan (BGN), Lysyl oxidase (LOX), and Collagen type V alpha 1 chain (COL5A1). Among these genes, FN1 had the highest node degree (82).

Construction of Hypoxia-Stemness Risk Signature
To further identify the essential genes (Hub DEGs) in the hypoxia-stemness-related DEGs sets, we extracted the clinical information of TCGA-ESCA patients, combining univariate Cox regression analysis and LASSO regression analysis to screen out essential genes (Hub DEGs) (Figure 3A, 3B). The results showed that a total of 14 prognosis-related genes were screened by Cox univariate regression ( Table 2). We randomly divide TCGA-ESCA patients into a training set and a validation set at a ratio of 4: 1. As shown in Figure 3C, the optimal risk score value corresponding to each esophageal tumor sample was obtained by LASSO regression analysis for follow-up analyses (Risk score=FBLN2 x (-0.1731926)+IL17RB×0.01043405+CYP2W1× 0.03672903+AMTN× 0.09540487+FABP1× 0.20776812+FOXA2× 0.21556933+GAS1× 0.2433375+CTSF× 0.2849325). Two sets of TCGA-ESCA patients are further divided into high-risk groups and low-risk groups according to the risk score individually. Kaplan-Meier curves were used for survival analysis. As shown in Figure 3D and 3E, the overall survival of high-risk esophageal cancer patients in the training set (P=0.005) and validation set (P=0.041) was significantly lower than that of the lowrisk group. In addition, the ROC curve analysis of the training

DATABASE ANALYSIS
set and the validation set indicates that the model has effective verification performance (Figure 3G, 3H).

External Dataset Validated the Robustness of the Hypoxia-Stemness Risk Classifier
The hypoxia-stemness risk model was further validated by an external independent cohort (GSE53624). Using the risk score classifier, Kaplan-Meier analysis was used to evaluate the survival difference between groups to verify the stability of the hypoxia-stemness risk model. As shown in Figure 3F, similar to the findings in the discovery cohort, the external validation cohort also showed high-risk groups with poor survival and low-risk groups with better outcomes (log-rank test, P=0.021). We analyzed the ROC curve to further verify the accuracy of the risk score model (Figure 3I). The results showed that the risk score model had strong predictive power for the prognosis of patients with esophageal cancer.

Hypoxia-Stemness Risk Classifier and Clinicopathologic Characteristics in Esophageal Cancer
Univariate Cox regression analysis showed that the clinical features of tumor stage (HR=2.671, P<0.001), N stage (HR=1.819, P<0.001), M stage (HR=4.847, P<0.001) and risk score (HR=3.346, P<0.001) were significantly related to survival rate in esophageal cancer patients, while multivariate Cox regression analysis indicated risk score (HR=4.172, P<0.001) to be an independent risk factor of overall survival in patients with esophageal cancer ( Table 3). Next, we divided esophageal cancer patients into different groups according to the clinicopathologic characteristics, and the corresponding risk score values were included in different groups of each sample. The results of each clinical feature subgroup are displayed in the form of box plots. The result of the Kruskal-Wallis rank sum test indicates that the distribution of risk score was significant between the 2 clinical characteristics of N stage (P=1.939e-02) ( Figure 4C) and fustat (P=3.3939e-03) (Figure 4E). There were no significant relations between risk score and stage (Figure 4A), T stage ( Figure 4B) or M stage ( Figure 4D). The comprehensive results show that the risk score derived from the hypoxia-stemness model has great applicability and stability for the clinicopathologic classification of esophageal cancer patients.

Multi-Omics Analysis Validated the Clinical Predictive Value of the Hypoxia-Stemness Risk Classifier
First, based on the transcription data, we used CIBERSORT and xCell algorithms to evaluate the distribution of tumor-infiltrating immune cells. Utilizing the xCell algorithm, the risk score was significantly positively correlated with hepatocytes, NKT, preadipocytes, plasma cells, and macrophages M2. In contrast, a significantly negative correlation was found between the risk score and keratinocytes, sebocytes, cDC, mast cells, monocytes, epithelial cells, CD4+ T cells, and myocytes ( Figure 5A). Meanwhile, the result of CIBERSORT was employed to generate a boxplot of immune cells in all tumor samples based on the risk score group. The risk score was significantly positively correlated with T cells regulatory (Tregs), B cells naïve, and T cells CD4 memory resting. On the other hand, there was   Figure 1A). In addition, we also observed increased CD8 expression ( Figure 5B) and decreased CD274(PD-L1) (Figure 5C) expression in high-risk groups compared to low-risk groups. However, we failed to find any significant difference between the high-and low-risk groups for tumor mutation burden (TMB) (P=0.083) (Supplementary Figure 1B), microsatellite instable (MSI) (P=0.18) (Supplementary Figure 1C) or T cell dysfunction and exclusion (TIDE) (P=0.65) (Supplementary Figure 1D).
In this study, drug response data were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) dataset, and we used the R package "pRRophetic" to predict the chemotherapy sensitivity of each tumor sample. The results show that the level of risk score is significantly related to the sensitivity e934359-11  of patients to Gefitinib (P=2.2e-07) (Figure 5F), Bexarotene (P=1.7e-06) (Figure 5G), Dasatinib (P=2e-04) (Figure 5H), and Imatinib (P=0.004) (Figure 5I). While there was no significant relationship between risk score and Metformin (Figure 5D) or Mitomycin C (Figure 5E).
We further explored the mutation profiles of patients in the highrisk and low-risk groups. The results showed that the proportion of TP53 and other gene mutations in the high-risk group was significantly higher than that in the low-risk group (Figure 5J).

GSVA Analysis Revealed a Close Relationship Between the Hypoxia-Stemness Risk Classifier and Tumor Proliferation
Since the constructed hypoxia-stemness risk model is significantly related to the prognosis of patients with esophageal cancer, it is valuable to explore its potential molecular mechanisms. GSVA is an improved GSE method, which is essentially a non-parametric and unsupervised method. Its value lies in integrating traditional analysis methods such as clustering, survival analysis, and correlation analysis into a pathwaycentric analysis [18]. As shown in Figure 6, different pathways of esophageal cancer patients in high-risk and low-risk groups are mainly enriched in angiogenesis, G2M checkpoint, glycolysis androgen response, protein secretion, xenogamy metabolism, and P53 pathway. It is suggested that the disturbance of these signaling pathways in the high-risk and low-risk groups can activate the proliferative process and ultimately affects the prognosis of patients with esophageal cancer.

Discussion
Hypoxia occurs in tumors when oxygen consumption exceeds oxygen supply in the tumor tissue. Intra-tumor hypoxia can promote the malignancy and aggressiveness of tumors and further adversely affect the clinical prognosis of cancer patients, especially in solid tumors [19]. However, uncontrolled proliferation of tumors can stimulate the growth of microvasculature. However, unlike the precise distribution of blood vessels in normal tissues, neovascularization in tumor tissues is characterized by disorderly distribution, and excessive misdistribution of tumor capillaries induced by hypoxic microenvironment can increase the distance between capillaries, thus exceeding the upper limit of oxygen exchange and diffusion capacity between capillaries [20]. Hypoxic areas within the tumor have been identified as one of the independent prognostic factors for cancer. Differences in the severity and duration of hypoxia result in oxygen gradients within the tumor, leading to heterogeneity and plasticity and promoting a more metastatic and aggressive phenotype [21]. More importantly, hypoxia mediates resistance to conventional chemotherapy and radiotherapy, which is ultimately associated with poor prognosis in cancer patients [22]. A growing number of studies have shown that hypoxia-mediated drug delivery contributes significantly to the development of chemoresistance in many carcinomas, including esophageal cancer [23]. As a critical mediator of the hypoxic response in the tumor microenvironment, the crucial role of hypoxia-inducible factors (HIFs) in resistance to conventional anti-cancer therapy for esophageal cancer has been demonstrated [24]. A recent study showed that high expression levels of the long noncoding RNA (lncRNA) E2F1 messenger RNA (mRNA) stabilizing factor (EMS), Wilms' tumor 1-associating protein (WTAP), and low expression levels of miR-758-3p were associated with poorer survival in patients with esophageal cancer. In addition, they demonstrated that the interaction between lncRNA-EMS/miR-758-3p/WTAP might be a potential mechanism for hypoxia-mediated cisplatin resistance in patients with esophageal cancer [25].
Cancer stem cells (CSCs) are defined as a small subpopulation of malignant cells at the top of the cell pyramid, characterized by the ability to self-renew and to give rise to differentiated cells [26]. Several clinical trials have shown that current surgical procedures alone or in combination with neoadjuvant therapy do not significantly improve the prognosis of patients with esophageal cancer [27][28][29]. The presence of CSCs is one of the leading causes of cancer recurrence and is mainly responsible for treatment resistance, which explains the inadequate response to conventional clinical treatment for esophageal cancer [30]. CSCs in esophageal cancer are directly involved in the regulation of carcinogenesis, tumor proliferation, distant metastasis, immune evasion, and drug resistance [31]. Notch, Wnt/ b-catenin, Hedgehog, Hippo pathways, and immune system interactions have been shown to be central mechanisms involved in regulating the stemness of CSCs [32]. Intratumoral hypoxic regions are thought to be niches where treatment-resistant and more aggressive cell subtype, cancer stem cells (CSCs), localize and accumulate [33]. CSCs-associated genes, including Lgr5, ALDH1A1, NANOG, SOX2, and p75NTR, were reported to be upregulated in ESCC under Cocl2 mimic hypoxic conditions. In addition, they demonstrated that targeting hypoxia by specific knockdown of HIF-1a inhibits cell proliferation and spheroid formation, suppresses the expression of CSC-related genes, and inhibits the activity of the Wnt/b-catenin pathway [34].
Considering the widely varying prognosis of patients with esophageal cancer, it would be valuable to establish an effective classification system to stratify patients with different levels of risk. Currently, due to the complexity of the tumor microenvironment (TME), few transcriptional expression profile-based biomarkers have been developed and assessed for their prognostic value in ESCC. Solid tumors are highly disorganized normal organs with many cell types, including fibroblasts, endothelial cells, immune cells, and malignant tumor cells. Complex metabolite exchange mechanisms within the tumor synergistically manage the metabolism of the tumor [35]. Therefore, it is difficult for a single biomarker to effectively assess the state of hypoxia within a tumor. Since hypoxia and stemness in the tumor microenvironment are symbiotic, there are no published biomarkers based on transcriptional expression patterns to jointly assess their status in esophageal cancer. In our study, mRNA expression data of hypoxic genes were used to verify esophageal cancer subtype assignment based on the T-SNE approach, and differentially expressed genes (DEGs) among esophageal cancer subtypes were obtained by the Limma method. Regarding stemness status, the idea of a stemness index was adopted for grouping patients with esophageal cancer and further to obtain differentially expressed genes (DEGs) associated with stemness index. In the present study, a total of 242 hypoxia-stemness index overlap DEGs were selected. Prognostic prediction models were further constructed by integrating differential genes associated with hypoxia-stemness index using LASSO regression. Finally, 8 prognosis-related genes (FBLN2, IL17RB, CYP2W1, AMTN, FABP1, FOXA2, GAS1, and CTSF) were selected to construct a risk scoring system. A risk score formula was built for each esophageal cancer patient, and the patients were divided into low-and high-risk groups. The stability of the prediction model was finally validated by assessing the survival differences between the 2 subgroups in both internal (TCGA-ESCC) and external datasets (GSEA53624), and univariate and multivariate Cox analyses showed that only the risk score was an independent prognostic variable in the TCGA-ESCC cohort.
The vital role of risk signature genes identified in this study has been previously reported in several types of cancer. The

DATABASE ANALYSIS
human FBLN2 gene encodes a protein that is a member of the fibronectin family. Fibulin-2 was initially found to be expressed at the site of the conversion of polarized epithelial cells into mesenchyme during endocardial cushion tissue development. In addition, fibulin-2 has been reported to be involved in cell migration and adhesion activities [36]. Notably, an isomer of FBLN2, FBLN2S, has previously shown tumor-suppressive and anti-angiogenic effects in nasopharyngeal carcinoma [37]. These results are consistent with the results of the present study, in which FBLN2 was identified as a gene associated with good prognosis of esophageal cancer. In contrast, overexpression of FBLN2 is related to unfavorable pathological features and poor prognosis in urethral epithelial cancer patients [36]. Another risk gene, IL17RB, has been reported to enhance cancer stemness and resistance to the chemotherapeutic drug gemcitabine in pancreatic cancer [38]. Moreover, FoxA2, a forkhead box transcription factor, is a member of the FoxA subfamily, and the expression level of FoxA2 was significantly increased under hypoxia environment in prostate adenocarcinomas. More importantly, the interaction between FpxA2 and HIF-1a stimulates a specific subset of HIF target genes, including Sox9, Jmjd1a, and Plod2, in the presence of hypoxia [39]. However, the above-mentioned signature genes have rarely been studied in the context of combined hypoxia and cancer stemness. The risk signature genes identified in this study may provide targets for elucidating the molecular mechanisms of esophageal cancer.
In this study, we also conducted a multi-omics study to investigate the clinical assessment validity of the risk model. We observed a positive correlation between high-risk scores and M2 macrophages, Tregs. As a drug delivery vehicle, macrophages play an essential role in hypoxia-mediated immune evasion due to their targeting tropism to the intra-tumor hypoxic microenvironment. M2 macrophages, also known as reparative macrophages, are involved in constructive tissue processes, including wound healing and tissue repair. In addition, M2 macrophages evade the immune surveillance system in the organism by secreting various immunosuppressive cytokines [40]. We further explored the mutation profiles of patients in the highand low-risk groups, and the results showed that the proportion of TP53 mutations was significantly higher in the highrisk group than in the low-risk group. Early studies have shown that tumors with TP53 mutations have higher hypoxia scores in various tumor types, including breast and lung adenocarcinoma, thus supporting the idea that TP53 mutations are a genomic consequence of tumor hypoxia [41]. In addition, multiple enriched biological processes and signaling pathways, including EMT, significantly associated with risk scores were identified. Esophageal cancer cells undergo EMT through activation of the Wnt/beta-catenin, Notch, and Hedgehog pathways in the tumor hypoxic microenvironment, acquiring characteristics of CSCs such as enhanced treatment resistance and malignancy. This makes EMT a vital bridge between hypoxia and cancer stem cells [42]. It is known that intratumorally hypoxia is associated with treatment resistance. The presence of cancer stem cells (CSCs) is also critical for tumor progression and chemotherapy resistance. More importantly, hypoxia affects patient prognosis not only by inducing tumor regeneration after treatment, but also by promoting the dissemination of stem cells [43]. In our study, patients with esophageal cancer in the high-risk group were significantly less sensitive to gefitinib, bexarotene, dasatinib, and imatinib, compared with the low-risk group. A randomized trial of second-line treatment for esophageal cancer reported that patients with epidermal growth factor receptor (EGFR) amplification had a better prognosis after treatment with gefitinib [44]. The relationship between EGFR and hypoxia has been demonstrated in various tumors [45,46].

Conclusions
In conclusion, 282 overlapping DEGs based on the hypoxiastemness index were screened for the first time in ESCA, and 8 prognosis-related genes (FBLN2, IL17RB, CYP2W1, AMTN, FABP1, FOXA2, GAS1, and CTSF) were selected to construct a risk scoring system. Furthermore, this novel risk scoring system effectively predicted the clinical prognosis of ESCA patients, including TCGA-ESCA and GSEA53624. However, a limitation of this study is that the robustness of our risk scoring model still needs to be further validated in an external independent ESCA cohort.

Data Availability Statements
Some or all data, models, or code generated or used during this study are available from the corresponding author by reasonable request.