Development and Validation of a Ferroptosis-Related Gene Signature and Nomogram for Predicting the Prognosis of Esophageal Squamous Cell Carcinoma

Esophageal squamous cell carcinoma (ESCC) is a common malignant tumor with high mortality and poor prognosis. Ferroptosis is a newly discovered form of cell death induced by iron-catalyzed excessive peroxidation of polyunsaturated fatty acids (PUFAs). However, the prognostic value of ferroptosis-related genes (FRGs) for ESCC remains unclear. Based on the ESCC dataset from the Gene Expression Omnibus (GEO) database, we identified 39 prognostic FRGs through univariate Cox regression analysis. After LASSO regression and multivariate Cox regression analyses, a multigene signature based on 10 prognostic FRGs was constructed and successfully divided ESCC patients into two risk groups. Patients in the low-risk group showed a significantly better prognosis than patients in the high-risk group. In addition, we combined the risk score with clinical predictors to construct a nomogram for ESCC. The predictive ability of the nomogram was further verified by ROC curves and calibration plots in both the training and validation sets. The predictive power of the nomogram was demonstrated to be better than that of either the risk score or clinical variable alone. Furthermore, functional analysis revealed that the 10-FRG signature was mainly associated with ferroptosis, differentiation and immune response. Connectivity map analysis identified potential compounds capable of targeting FRGs in ESCC. Finally, we demonstrated the prognostic value of SRC gene in ESCC using the clinical samples and found that SRC inhibition sensitized ESCC cells to ferroptosis inducers by in vitro experiments. In conclusion, we identified and verified a 10-FRG prognostic signature and a nomogram, which provide individualized prognosis prediction and provide insight into potential therapeutic targets for ESCC.

Esophageal squamous cell carcinoma (ESCC) is a common malignant tumor with high mortality and poor prognosis. Ferroptosis is a newly discovered form of cell death induced by iron-catalyzed excessive peroxidation of polyunsaturated fatty acids (PUFAs). However, the prognostic value of ferroptosis-related genes (FRGs) for ESCC remains unclear. Based on the ESCC dataset from the Gene Expression Omnibus (GEO) database, we identified 39 prognostic FRGs through univariate Cox regression analysis. After LASSO regression and multivariate Cox regression analyses, a multigene signature based on 10 prognostic FRGs was constructed and successfully divided ESCC patients into two risk groups. Patients in the low-risk group showed a significantly better prognosis than patients in the high-risk group. In addition, we combined the risk score with clinical predictors to construct a nomogram for ESCC. The predictive ability of the nomogram was further verified by ROC curves and calibration plots in both the training and validation sets. The predictive power of the nomogram was demonstrated to be better than that of either the risk score or clinical variable alone. Furthermore, functional analysis revealed that the 10-FRG signature was mainly associated with ferroptosis, differentiation and immune response. Connectivity map analysis identified potential compounds capable of targeting FRGs in ESCC. Finally, we demonstrated the prognostic value of SRC gene in ESCC using the clinical samples and found that SRC inhibition sensitized ESCC cells to ferroptosis inducers by in vitro experiments. In conclusion, we identified and verified a 10-FRG prognostic signature and a nomogram, which provide individualized prognosis prediction and provide insight into potential therapeutic targets for ESCC.

INTRODUCTION
Esophageal cancer is the seventh most common cancer worldwide and ranks sixth in terms of cancer mortality (Bray et al., 2018). There are two main histological subtypes, esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC), which have almost completely distinct etiologies, geographic patterns, and biological characteristics (Bray et al., 2018). Despite the rapid increase in the incidence of EAC in Western countries, ESCC remains the predominant histological type of esophageal cancer in Eastern Asia, accounting for over 90% of all new esophageal cancer cases each year (Bray et al., 2018). ESCC has a poor 5-years overall survival (OS) rate and a high incidence of recurrence and metastasis (Pennathur et al., 2013). Although the tumor-node-metastasis (TNM) staging system serves as the standard method for predicting the OS of cancer patients, there still exist differences in the survival of patients with the same TNM stage (Rice et al., 2017). Therefore, ESCC prognosis prediction based on the TNM staging system needs further improvement. In recent decades, due to advancements in high-throughput technologies, such as microarray and RNA sequencing, gene expression profiling has been widely used to discover molecular biomarkers associated with the phenotype or prognosis of cancer (Zhan et al., 2016). Recently, multigene signatures, such as Oncotype DX for breast cancer or ColoPrint chips for colon cancer, have proven to be of great prognostic value for cancers, and the described signatures can be employed to guide the prognostic evaluation, treatment and management of cancers (Birnbaum et al., 2017). Many studies on the gene expression profiles of ESCC have been reported, but no gene signature has been applied for survival prediction of this disease thus far. Hence, there is an urgent need to identify key molecular biomarkers related to the prognosis of ESCC.
Ferroptosis is an iron-catalyzed form of cell death that occurs through excessive peroxidation of polyunsaturated fatty acids (PUFAs) (Dixon et al., 2012). In recent years, several small molecules and FDA-approved clinical drugs have been identified to induce ferroptosis in cancer cells (Hassannia et al., 2019). The efficacy of cancer suppression by ferroptosis inducers in various studies highlights the potential of ferroptosis as a novel anticancer strategy (Hassannia et al., 2019). A previous study reported that upregulation of glutathione peroxidase 4 (GPX4), a key negative regulator of ferroptosis, and downregulation of heme oxygenase-1 (HMOX1), which promotes ferroptosis by increasing the labile iron pool (LIP), were poor prognostic factors in ESCC (Shishido et al., 2020). 5-Aminolevulinic acid (5-ALA) induces ferroptosis by modulating GPX4 and HMOX1 in ESCC (Shishido et al., 2020). Recently, several studies have revealed the predictive value of ferroptosisrelated gene (FRG) signatures in various cancers (Liang et al., 2020;Luo and Ma, 2021;Zheng et al., 2021). However, whether these FRGs are correlated with the prognosis of ESCC patients remains unknown.
In this study, we collected the mRNA expression profiles and clinical information of ESCC patients from the Gene Expression Omnibus (GEO) database. Then, we constructed a prognostic risk model based on FRGs to predict the prognosis of ESCC. In addition, we developed a nomogram model combining the risk score and clinical features to assess prognosis. The prognostic value of the risk model and the nomogram was then validated in another independent dataset obtained from the GEO database. Furthermore, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) were used to reveal the underlying biological characteristics and signaling pathways associated with this signature in ESCC. We analyzed the difference in immune cell infiltration in diverse subgroups via CIBERSORT and ImmuCellAI. Finally, we validate the prognostic value and ferroptotic role of FRGs in ESCC through a series of experiments. The overview workflow is presented in Figure 1.

Data Acquisition and Collection of Ferroptosis-Related Genes
information for 60 ESCC cases, was used for validation. A total of 19,631 mRNAs were obtained after annotation. A comprehensive list containing 352 FRGs was retrieved from FerrDb (Zhou and Bao, 2020) (http://www.zhounan.org/ferrdb/index.html) and previously published literature  and is provided in Supplementary Table S1.

Construction and Validation of the Prognostic Risk Model
Univariate Cox regression analysis was performed to screen potential prognostic FRGs, and FRGs with p < 0.1 were considered statistically significant and incorporated into the subsequent least absolute contraction and selection operator (LASSO) regression analysis. LASSO regression analysis was applied to minimize the potential risk of overfitting and select optimal prognostic genes. Ten-times cross-validations were utilized to determine the best penalty parameter lambda (λ). Finally, we performed a multivariate Cox regression for these genes to generate a prognostic signature. The prognostic risk score was determined using a linear combination of the regression coefficient (β) in the multivariate Cox regression model and the expression levels of the genes. In the formula, risk score (β1*expression of gene X1) + (β2*expression of gene X2) + (βi*expression of gene Xi). Patients were divided into lowrisk and high-risk groups based on the median risk score. The difference in OS between the two groups was compared with Kaplan-Meier survival curves and log-rank tests. Furthermore, receiver operating characteristic (ROC) curves were plotted, and area under the curve (AUC) values were calculated to evaluate the predictive power of the gene signature using the "survivalROC" package in R. Moreover, the same analytical methods were performed for the validation cohort to evaluate the prognostic capacity of the gene signature.

Nomogram Construction and Evaluation
Univariate and multivariate Cox analyses were used to identify independent prognostic factors such as age, sex, alcohol use, grade, stage and risk score. Subsequently, a nomogram was constructed based on the results of multivariate Cox analysis to predict the 1-, 3-and 5-years OS of patients with ESCC. The predictive ability of the nomogram was then assessed by calibration curves using the "foreign" package. In addition, ROC curves and the corresponding AUC values were generated using the "survivalROC" package. Decision curve analysis (DCA) was performed to evaluate the clinical benefit that the nomogram can obtain compared to a single independent prognostic predictor. Then, the same operation was conducted for the validation set GSE53622 to evaluate the predictive ability of the nomogram.

Interaction Network and Functional Enrichment Analyses
An interaction network of prognostic FRGs was performed at the STRING website (http://string-db.org). The differentially expressed genes (DEGs) with a false discovery rate (FDR) < 0. 05 and a |log fold change (FC)| > 1 between the low-risk and highrisk groups were identified using the "limma" package in R. Then, the "clusterProfiler" package was utilized to conduct GO and KEGG analyses based on DEGs. GSEA was employed to identify the biological processes, molecular functions and signaling pathways enriched in the low-risk and high-risk groups. The KEGG gene set (C2. cp.kegg.v7.0 symbols. gmt) and GO gene set (C5. go.v7.0 symbols. gmt) were downloaded from the Molecular Signatures Database (MSigDB). The pathways with the following criteria were regarded as significantly enriched: nominal (NOM) p-value < 0.05 and FDR q-value < 0.25.

Correlation Between the Risk Score and Immune Infiltration
To explore the abundance of infiltrating immune cells in the highand low-risk groups, the CIBERSORT algorithm was used to score the infiltration abundance of each immune cell in the samples to evaluate the proportion of 22 types of immune cells in each sample. We then compared the infiltration levels of 22 types of immune cells between the low-and high-risk subgroups. In addition, ImmuCellAI was used to estimate the abundance of immune cell infiltration and predict the response of each sample to immune checkpoint blockade (ICB) therapy (Miao et al., 2020).

Identification of Candidate Small Molecules
The Connectivity Map (CMap) database (http://www. broadinstitute.org) was used to predict the potential compounds that might inhibit or induce biological states encoded by specific markers (Lamb et al., 2006). To explore the potential activity of small molecules from the CMap database in different subgroups, the prognostic ferroptosis-related genes were uploaded to the CMap database for mode-of-action (MoA) analysis.

Patients and Tissue Samples
Paraffin-embedded samples of ESCC, which were diagnosed clinically and pathologically, were collected from 96 patients between 2004 and 2008 in Meizhou People's Hospital, China. None of the patients received radiotherapy or chemotherapy before surgery, and none of them had multiple cancers in other organs. Prior informed consent was obtained from all patients, and this study was approved by the research Ethics Committee of Meizhou People's Hospital. The following clinicopathological parameters were collected from the medical records: age, sex, histological grade, depth of invasion, and clinical stage. The histopathological diagnosis was based on the World Health Organization criteria. Tumor staging was determined according to the 6th edition of the tumor-node-metastasis (TNM) classification of the International Union Against Cancer.
In addition, 20 ESCC tissues and their paired adjacent noncancerous esophageal epithelial tissues were attained from ESCC patients between 2019 and 2021 in The First Affiliated Hospital of Sun Yat-sen University. No patients received radiotherapy or chemotherapy before surgery, and none of them had multiple cancers in other organs. After surgical removal, fresh tissues were Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 697524 immediately snap-frozen in liquid nitrogen and stored at −80°C until total RNA extraction and analysis. Informed consent forms were obtained from all patients included in the present study.

Immunohistochemical Staining and Scoring
The paraffin-embedded samples of ESCC and normal esophageal tissues were cut into 5-μm-thick sections and placed on pathological slides for immunohistochemical staining. Tissue sections were heated at 100°C in citrate buffer solution (pH 6.0) for 10 min to facilitate antigen retrieval. Then, the sections were incubated with rabbit antibody against SRC (1:400, CST, United States) overnight in 4°C followed by incubation with secondary antibody (Dako REAL EnVision, United States). Immunoreacted cells were visualized using diaminobenzidine, and nuclei were counterstained with hematoxylin. Phosphatebuffered saline (PBS) was substituted for the primary antibody as a negative control. Sections were independently evaluated microscopically by two pathologists without knowledge of the clinicopathological features. SRC expression level was determined by integrating the percentage of positive tumor cells and the intensity of positive staining. The intensity of staining was scored as follows: negative (score 0), weak (score 1), moderate (score 2), and strong (score 3). The extent of staining was scored according to the percentage of positive stained tumor cells in the field: <5 (score 0), 5-25% (score 1), 26-50% (score 2), 51-75% (score 3), and 76-100% (score 4). The product of the intensity and extent score was considered as the final histochemistry score (H-score), yielding a range from 0 to 12. When the staining was heterogeneous, each component was scored independently and summed for the results. We defined final score 0-7 as low expression and 8-12 as high expression for SRC immunohistochemical staining.

Cell Viability Assay
Cell viability was detected using Cell Counting Kit-8 (CCK-8, Dojindo, Japan) assays according to the manufacturer's instructions. Briefly, the cells were seeded into a 96-well plate at a density of 5,000 cells/well. After treatment with different drugs at various concentrations for the indicated times, 10 µL of CCK-8 reagent was added to each well, and the cells were cultured for another 2 h. At the end of the incubation, the absorbance at 450 nm was analyzed with a microplate reader (BioTek, United States). All experiments were performed in triplicate.
Assessment of Lipid Peroxidation Using BODIPY-C11 (581/591) 2 × 10 5 cells per well were seeded on 12-well dishes overnight. After treatment with the indicated concentration of IKE and/or KX2-391 for 18 h to induce ferroptosis, cells were incubated with BODIPY-C11 (581/591) (1 μM) for 30 min at 37°C before they were harvested by trypsinisation. Subsequently, cells were resuspended in 300 μL of fresh PBS and analyzed using the 488-nm laser of flow cytometer (FACS Canto II, BD Biosciences) for excitation. Data was collected from the FITC detector for oxidized BODIPY-C11. A minimum of 10,000 cells were analyzed per sample. Data was analyzed using FlowJo V10 Software.

Western Blotting
Western blot analysis was performed to examine the expression of various proteins as described in our previous report (Ye et al., 2017). GAPDH was used as the control. The primary antibodies used included the following: GAPDH, SRC, p-SRC Tyr416 and PTGS2 (1:1,000 dilution; Cell Signaling Technology, United States).

Statistical Aanalysis
The association between SRC expression and clinicopathological variables was assessed using a chi-square test or Wilcoxon rank sum test. The Kaplan-Meier method and log-rank tests were used to compare the overall survival. Multivariate analyses of variables were conducted using a Cox proportional hazards regression model. The difference between two groups was analyzed by a twotailed Student's t test, while values were compared among multiple groups using one-way ANOVA. All analyses were considered statistically significant when p < 0.05 was obtained.  Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 697524 5 were subjected to univariate Cox regression analysis. Thirty-nine FRGs were significantly associated with the prognosis of ESCC (p < 0.1). To narrow the number of FRGs, LASSO regression analysis was then performed on these 39 FRGs. As a result, a total of 15 FRGs markedly related to prognosis of ESCC were obtained (Figures 2A,B). Subsequently, multivariate Cox regression analysis was used to select the best characteristic gene set and construct a regression model. Finally, the prognostic signature was constructed based on 10 ferroptosisrelated genes (SRC, FADS2, GLUD1, POLG, ANO6, SLC2A6, ANGPTL7, PTGS2, ALOXE3 SLC38A1) ( Figure 2C). The formula of the risk model was as follows: risk score 0.741 × SRC + 0.633×FADS2 + 0.544 × GLUD1 + 0.527×POLG + 0.477 × ANO6 + 0.295×SLC2A6 + 0.252 × ANGPTL7−0.202 × PTGS2−0.24 × ALOXE3−0.57 × SLC38A1. The median risk score was used to dichotomize patients in the GSE53624 cohort into low-risk (n 60) and high-risk (n 59) groups. The heatmap showed that SRC, FADS2, GLUD1, POLG, ANO6, SLC2A6 and ANGPTL7 were highly expressed, while PTGS2, ALOXE3, and SLC38A1 were downregulated in high-risk cases ( Figure 2D). Kaplan-Meier curves revealed that the prognosis of patients in the low-risk group was significantly better than that of patients in the high-risk group (p < 0.001, Figure 2E). As shown in Figures 2F,G, the distribution of risk score and survival status indicated that a low risk score was beneficial to survival. The AUC values of the time-dependent ROC curves at 1, 3 and 5 years were 0.815, 0.833 and 0.833, respectively, indicating that the prognostic signature had great specificity and sensitivity for predicting the OS of ESCC patients ( Figure 2H).

Validation of the Prognostic Signature in the Independent Cohort
We validated the predictive ability of the prognostic signature in the independent cohort from the GSE53622 dataset. The patients were divided into low-risk (n 30) and high-risk (n 30) groups using the same prognostic risk model ( Figure 3A). Consistent with the above results, patients in the low-risk group had longer OS than those in the high-risk group ( Figure 3B, p < 0.01). The risk score distribution and survival status were similar to those in the training dataset GSE53624 (Figures 3C,D). The AUC values for 1-, 3-and 5-years survival were 0.492, 0.69 and 0.701, respectively ( Figure 3E). In view of these results, the 10-FRG signature based on the training set exhibited a certain power in predicting the OS of ESCC patients.

Construction and Validation of the Predictive Nomogram
To determine whether the predictive ability of the risk score was independent of other traditional clinical characteristics (including age, sex, tobacco use, alcohol use, tumor location, grade, T stage, N stage and TNM stage), we performed univariate and multivariate Cox regression analyses on these variables using a training set. The results showed that age (HR 1.030), grade (HR 1.667), N stage (HR 2.768) and risk score (HR 0.254) were independent prognostic factors ( Figures 4A,B). Based on these four independent predictive factors, we constructed a predictive nomogram to quantify the prediction results of individual survival probability at 1, 3 and 5 years ( Figure 5A). We then performed time-dependent ROC curve analysis to evaluate the predictive capacity of the nomogram. The AUCs for 1-, 3-and 5-years OS were 0.730, 0.797 and 0.806, respectively, in the training cohort, while the AUCs for 1-, 3-and 5-years OS in the validation set were 0.670, 0.758 and 0.779, respectively ( Figures 5B,C). The C-index for the nomogram was 0.717 (95% CI: 0.663-0.772). The calibration curves of both the training set and validation set showed high consistency between the actual proportion of 1-, 3-and 5-years OS and the nomogram-predicted probability ( Figures 5D,E). Finally, we performed DCA to assess the value of the nomogram in clinical decision making. We found that compared to a single independent predictive factor, the nomogram could obtain the optimal net benefit in both the training set and validation set ( Figures 5F,G). Overall, these results demonstrated that the developed nomogram preforms well in predicting OS.

Functional Analyses of Risk Model
We explored the prognostic FRGs interaction at the STRING online website, and the gene network demonstrated the SRC, PTGS2, CDH1 SMAD3, HSPA5 and PIK3CA were the hub genes ( Figure 6).
A total of 62 DEGs (|logFC| > 1, FDR <0.05) were identified between the high-and low-risk groups in the GSE53624 cohort. To elucidate the biological functions and pathways that were associated with the risk score, the DEGs were used to perform GO enrichment and KEGG pathway analyses. According to GO analysis, the DEGs were enriched in several differentiation-related biological processes, such as epidermis development, skin development, and keratinocyte differentiation ( Figure 7A). KEGG analysis revealed that the DEGs were enriched in ferroptosis-related signaling pathways, including alpha-linolenic acid (LA) metabolism and LA metabolism ( Figure 7B).
GSEA was further conducted between the two groups in the GSE53624 cohort. The results showed that the pathways of cancer, DNA replication, cell cycle and TGF-β signaling from the KEGG database were enriched in the high-risk group ( Figure 7C). In Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 697524 8 contrast, cytochrome P450, arachidonic acid (AA) metabolism, LA metabolism and the JAK-STAT signaling pathway, which play vital roles in ferroptosis, were enriched in the low-risk group ( Figure 7C). Among the GO terms, β-catenin TCF complex assembly, autophagosome organization, phosphatidylinositol-5-phosphate binding and WNT-activated receptor activity were enriched in the high-risk group ( Figure 7D). We found that biological processes associated with differentiation and immunity, such as keratinization, epidermal cell differentiation, humoral immune response and cell killing, were enriched in the low-risk group ( Figure 7D). In addition, chemokine-related molecular functions, such as chemokine receptor binding, chemokine activity and chemokine receptor (CCR) binding, were enriched in the low-risk group ( Figure 7D).

Relationship Between Risk Score and Immune Cell Infiltration
CIBERSORT was implemented to assess the abundance of 22 kinds of immune cell infiltrates in tumor samples ( Figure 8A). As shown in Figure 8B, plasma cells were downregulated in the high-risk group of the training cohort (p < 0.05). In addition, the risk score was negatively correlated with plasma cell abundance and positively correlated with CD8 T cell abundance ( Figure 8C). We also used ImmuCellAI to compare the infiltration levels of 24 kinds of immune cells between low-and high-risk groups of ESCC samples. Figure 8D shows that the proportion of natural T regulatory cells (nTregs) in the high-risk group was significantly higher than that in the low-risk group. Furthermore, ImmuCellAI was also applied to predict the response to ICB therapy and showed that the low-risk group had a better response to ICB than the high-risk group (p < 0.05, Figure 8E). In addition, the risk score was higher in patients with lymph node metastasis and advanced TNM stage (p < 0.05, Figure 8F).

Connectivity Map Analysis Identified Potential Compounds Suitable for Different Molecular Subtypes
Based on the risk score of the predictive signature, patients in the high-risk group tended to have poor survival outcomes, while Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 697524 10 those in the low-risk group tended to have better survival outcomes. Therefore, we attempted to identify potential compounds or inhibitors for ESCC patients in different risk groups. The 10 prognostic FRGs were uploaded to the CMap database for MoA analysis. Three FRGs, including FADS2, PTGS2 and SRC, were screened as potential targets for ESCC. A total of 75 compounds with 39 MoAs were enriched ( Figure 9A). Both SRC and PTGS2 were overexpressed in cancer tissues, and the corresponding inhibitors may be useful for ESCC treatment ( Figure 9B, p < 0.001). Because SRC was highly expressed in the high-risk group ( Figure 9C, p < 0.001), SRC inhibitors, such as bosutinib, dasatinib, PP-1 and PP-2, were predicted to be promising compounds for the treatment of patients in the high-risk group. Cyclooxygenase inhibitors targeting PTGS2 could be suitable for treating low-risk patients.

Clinical Experimental Validation
We performed the PCR and IHC validation in clinical specimens following the steps described above. We evaluated the expression of SRC, the most significantly prognostic FRG according to the coefficient value, in 96 ESCC tissues and 57 adjacent normal esophageal tissues by IHC. Representative tissues with IHC staining are shown in Figure 10A. The expression of SRC presented mainly in the cytoplasm, and thus, brown cytoplasm immunoreactivity for SRC was identified as positive staining. The median H-score of SRC expression in ESCC was 8, thus we defined H-score 0-7 as low expression and 8-12 as high expression for SRC at protein level. As shown in Figures  10A,B, SRC was positively expressed in all ESCC tissues, of which 57 cases (59.4%) show high expression and 39 cases (40.6%) show low expression. In contrast, most of the normal esophageal squamous epithelium showed negative expression of SRC, and none of the normal tissues showed high expression. Thus, it revealed that the SRC expression level in ESCC was significantly higher than in normal tissue (p < 0.001, Figure 10B). In survival analysis, Kaplan-Meier curves showed that patients with high SRC expression survived significantly shorten than patients with low SRC expression (log-rank test, p 2.666e-6, Figure 10C). The mean survival time of patients with low expression was 39.46 ± 5.29 months, but it decreased to 13.20 ± 1.17 months in patients with high SRC expression. Multivariate analysis revealed that high SRC expression (HR 3.36, 95%CI 1.900-5.940, p 0.000) was an independent prognostic factor for poor prognosis in ESCC. The correlations between SRC expression and clinicopathological variables was assessed and shown in Table 1. High expression of SRC was found to significantly correlate with histological grade (p 0.005). Meanwhile, we also evaluated the mRNA expression of SRC in 20 pairs of human ESCC tissues and corresponding non-cancerous tissues. The RT-PCR results showed that SRC was up-regulated in 11 cases (55.0%, Figure 10D). Moreover, SRC was significantly up-regulated in all ESCC cell lines compared with normal esophageal epithelial cell line HET-1A (p < 0.001, Figure 10E).

SRC Protects Esophageal Squamous Cell Carcinoma Cells Against Ferroptosis
To investigate the potential role of SRC in ferroptosis, ESCC cell lines including KYSE150 and ECA9706 were treated with RSL3 and IKE, two well-established ferroptosis agonists. The immunoblotting results of phosphorylated SRC (Tyr416) showed that SRC was activated in response to low concentration of IKE or RSL3 but inhibited under high dose of IKE or RSL3 in ECA9706 cells ( Figure 11A). Whereas the expression level of total SRC was not affected when treated with RSL3 or IKE ( Figure 11A). In addition, low dose of RSL3 and IKE inhibited the expression of phosphorylated SRC (Tyr416) in KYSE150 cells, which was more sensitive to ferroptosis than ECA9706 cells ( Figure 11A). Moreover, ferroptosis agonists significantly increase the expression level of PTGS2 in a dose-dependent manner ( Figure 11A). To evaluate whether SRC protects ESCC cells against ferroptosis, ESCC cells were treated with the SRC inhibitor KX2-391 and assayed for viability after 48 h. SRC inhibition decreased cell viability, and that loss of viability was rescued by ferroptosis inhibitor liproxstatin-1 ( Figure 11B). Moreover, SRC inhibition sensitized ESCC cells to IKE and increased IKE-induced lipid ROS generation (Figures 11C,D). However, PTGS2 inhibition with celecoxib did not inhibit cell viability and ferroptosis sensitivity ( Figures 11E,F).

DISCUSSION
A series of studies have shown that ferroptosis is essential for eradicating cancer cells and that sensitivity to ferroptosis varies in different types of cancers (Xu et al., 2019). Ferroptosis initiation and execution lie at the intersection of glutathione metabolism, lipid peroxidation of PUFAs, iron metabolism and mitochondrial function (Stockwell et al., 2017). A previous study reported that upregulation of GPX4 and downregulation of HMOX1 were poor prognostic factors for ESCC (Shishido et al., 2020). Another study found that DNAJB6 level was negative related to lymph node metastasis in ESCC patient (Jiang et al., 2020). Overexpressing DNAJB6a showed tumor-suppressive effects in vitro and in vivo. However, the prognostic value of FRGs in ESCC has yet to be comprehensively clarified. In this study, we constructed a prognostic risk model based on 10 FRGs, which comprised 3 protective genes (ALOXE3, SLC38A1 and PTGS2) and 7 riskrelated genes (SRC, ANGPTL7, ANO6, SLC2A6, GLUD1, POLG and FADS2). Therefore, patients with ESCC can be classified into low-risk and high-risk groups for discrimination of survival outcomes. Patients in the low-risk group showed better survival than those in the high-risk group. ROC cures also suggested the prediction ability of 10-FRG signature.
Among the 10 genes, ALOXE3, SLC38A1 and PTGS2 were negatively correlated with the risk score and downregulated in the high-risk group. A previous study showed that silencing arachidonate lipoxygenase E3 (ALOXE3) in HT-1080 cells made cells resistant to erastin, which supports that ALOXE3 is required for erastin-induced ferroptosis (Yang et al., 2016). Another study revealed that ALOXE3 is markedly downregulated in human glioblastoma (GBM) and that ALOXE3 deficiency renders GBM cells resistant to ferroptosis, promoting GBM cell survival (Yang et al., 2021). Gao et al. found that knockdown of SLC38A1 markedly blocked ferroptosis, suggesting that SLC38A1 positively regulates ferroptosis (Gao et al., 2016). PTGS2, a gene encoding cyclooxygenase-2 (COX-2), was upregulated upon treatment with either erastin or RSL3 (Yang et al., 2014). The role of PTGS2 in ferroptosis was also validated in our study. Treatment with either IKE or RSL3 increased the expression of PTGS2. However, ferroptotic cell death by IKE was not affected by celecoxib treatment, suggesting that PTGS2 upregulation is simply a downstream marker of ferroptosis. This is consistent with the previous study (Yang et al., 2014). Therefore, these three protective genes are positively related to ferroptosis.
In contrast, SRC, ANGPTL7, ANO6, SLC2A6, GLUD1, POLG and FADS2 are positively related to the risk score and serve as risk-related genes. A previous study revealed that activation of SRC can protect cancer cells from ferroptosis by suppressing the expression of ACSL4, an enzyme that enriches membranes with PUFAs and is required for ferroptosis, revealing the inhibitory role of SRC in ferroptosis (Brown et al., 2017). ANGPTL7 may inhibit ferroptosis because its expression is downregulated during ferroptosis induced by erastin or RSL3 (Yang et al., 2014). Depletion of FADS2 increased erastin-induced ferroptosis by reducing the expression of the ferroptosis-associated regulator GLUD1 at the mRNA level, suggesting the negative role of FADS2 and GLUD1 in ferroptosis regulation (Jiang et al., 2017). A recent study indicated that DNA polymerase γ (POLG) expression is upregulated in pancreatic ductal adenocarcinoma patients and serves as a therapeutic target to induce ferroptosis by zalcitabine . Thus, most risk-related genes have been demonstrated to negatively regulate ferroptosis, promoting cancer growth and development. For most cancers, including ESCC, pathologic TNM staging has been identified as a prognostic indicator and helps to guide the selection of therapeutic strategies. Nevertheless, due to differences in molecular and genetic characteristics, the clinical outcomes and prognoses of cancers vary even among patients with the same TNM stage who receive similar treatments (Gerlinger et al., 2012). Consequently, survival prediction based only on the TNM stage is not fully satisfactory to physicians, suggesting that prognosis prediction should be based on both the TNM stage and molecular characteristics of cancer. In our present study, the AUC value for 1-year survival was 0.492 in validation cohort, indicating the predictive capability of the 10-FRG signature for OS was limited. To personalize prognosis prediction and help to better design the treatment strategy of patients, we combined the risk score with clinicopathological characteristics to construct a nomogram for ESCC. In the univariate Cox regression analysis, age, tobacco use, grade, N stage and the risk score were significantly associated with OS. Furthermore, multivariate Cox analysis showed that the risk score was an independent prognostic factor. We constructed a nomogram based on risk score and clinical factors. The predictive ability of the nomogram was further verified by time-dependent ROC curves and calibration plots in both the training and validation sets. The result showed a higher AUC (1-, 3-and 5-year OS was 0.67, 0.758 and 0.779, respectively) than the 10-FRG risk score alone in validation dataset. The predictive power of the nomogram was demonstrated to be better than that of either the risk score or clinical variable alone. Therefore, physicians can apply the nomogram to improve the accuracy of identifying high-risk patients and realize accurate treatment.
During ferroptosis, enzymatic lipid peroxidation is mediated by the activity of the lipoxygenase (LOX) family (Hassannia et al.,  . LOXs are nonheme iron-containing enzymes that catalyze the deoxygenation of free and esterified PUFAs to generate various lipid hydroperoxides, which can cause the destruction of the lipid bilayer and affect membrane function (Hassannia et al., 2019;Chen et al., 2021). In mammalian cells, LA and AA are the most abundant PUFAs serving as substrates for LOXs (Hassannia et al., 2019). In this study, both AA metabolism and LA metabolism were enriched in the low-risk group, suggesting that ESCC patients in the low-risk group are susceptible to ferroptosis. Thus, ferroptosis inducers may be more valuable to patients with low risk. In contrast, pathways related to carcinogenesis and development, such as colorectal cancer, DNA replication and the cell cycle, were positively correlated with the high-risk group. This result suggested that the risk score was positively associated with malignancy of ESCC, which was also supported by the findings that the risk score was higher in patients with lymph node metastasis and advanced TNM stage. Biological processes involving tumor differentiation and suppression were enriched in the low-risk group. This result was consistent with the result of GO enrichment analysis.
Chemokines are a large family of small, secreted proteins that interact with cell surface G protein-coupled receptors to stimulate the migration of leukocytes during normal immune function (Hughes and Nibbs, 2018). Consequently, chemokines play a central role in the development and homeostasis of the immune system (Hughes and Nibbs, 2018). Chemokine-related molecular functions, including chemokine receptor binding, CCR binding and chemokine activity, were enriched in the low-risk group, suggesting that chemokines are active in patients with low-risk scores.
Increasing evidence has shown that ferroptosis is associated with tumor immunity (Friedmann Angeli et al., 2019). Wang et al. reported that immunotherapy-activated CD8 + T cells induce ferroptosis in tumor cells in vivo (Wang et al., 2019). CD8 + T cellderived IFN-γ downregulated the expression of SLC3A2 and SLC7A11 by inhibiting its transcription, thus promoting lipid peroxidation and ferroptosis in tumor cells (Wang et al., 2019). Furthermore, depletion of cystine or cysteine in combination with ICB synergistically enhanced T cell-mediated antitumor immunity and induced ferroptosis in a mouse model (Wang et al., 2019). Immunohistochemical studies showed that the CD8 + T cell signature was negatively associated with system Xcexpression, suggesting that sensitivity to ferroptosis was parallel to anticancer immunity (Wang et al., 2019). In our study, the CD8 + T cell level in ESCC was positively correlated with the risk score, indicating that the low-risk group may be more sensitive to ICB. This was also supported by the results of ImmuCellAI analysis, which revealed that the response rate to ICB therapy was significantly higher in the low-risk group than in the high-risk group. Thus, ICB and ferroptosis inducers may be suitable for patients with low risk. Furthermore, we used CMap to discover personalized treatment options for ESCC patients and found that SRC inhibitors can be utilized for the high-risk group. The protein encoded by SRC is a tyrosine-protein kinase that plays a role in the regulation of embryonic development, cell growth and multiple fields of tumorigenesis (Roskoski, 2015). SRC has been indicated as a promising therapeutic target in the treatment of solid tumors, including ESCC (Roskoski, 2015). For example, bosutinib, dasatinib, and ponatinib are SRC/multikinase inhibitors that are approved by the FDA for the treatment of chronic myelogenous leukemia (Chen et al., 2015). Dasatinib enhances cisplatin sensitivity in ESCC cells via suppression of the PI3K/AKT and Stat3 pathways (Chen et al., 2015). Bosutinib effectively induces apoptosis in ESCC cells by inhibiting Src/Abl signaling (Ha et al., 2020). In our present study, we demonstrated SRC was overexpressed in ESCC tissues compared with normal esophageal mucosa. In addition, high expression of SRC was significantly associated with advanced histological grade and poor prognosis in ESCC. This was consistent with the results of bioinformatic analysis in our study. Moreover, we also revealed that ferroptosis inducers activated SRC while SRC inhibition sensitized ESCC cells to ferroptosis. Thus, we propose that SRC plays an essential role in ESCC progression and ferroptosis, and it may be used as a potential prognostic marker and therapy target for ESCC.
In summary, we constructed and validated a 10-FRG signature-based risk model and a nomogram that could be used to predict the prognosis of ESCC. Moreover, our study provides a new understanding of ferroptosis in the context of ESCC carcinogenesis and progression and offers important ideas for developing ferroptosis inducers for the treatment of ESCC. Because our results are based on microarray technology, a series of in vitro and in vivo experiments are required to advance the clinical application of our 10-FRG signature to improve the survival rate of ESCC patients.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
AH designed the study and revised the manuscript. JY, YW, HC, WD, and RL collected and assembled the data. JY and YW performed data analysis and interpretation. LS was responsible for IHC. JY drafted the manuscript. All the authors read and approved the final manuscript.