Seven Fatty Acid Metabolism-Related Genes as Potential Biomarkers for Predicting the Prognosis and Immunotherapy Responses in Patients with Esophageal Cancer

Background: Esophageal cancer (ESCA) is a major cause of cancer-related mortality worldwide. Altered fatty acid metabolism is a hallmark of cancer. However, studies on the roles of fatty acid metabolism-related genes (FRGs) in ESCA remain limited. Method: We identified differentially expressed FRGs (DE-FRGs). Then, the DE-FRGs prognostic model was constructed and validated using a comprehensive analysis. Moreover, the correlation between the risk model and clinical characteristics was investigated. A nomogram for predicting survival was established and evaluated. Subsequently, the difference in tumor microenvironment (TME) was compared between two risk groups. The sensitivity of key DE-FRGs to chemotherapeutic interventions and their correlation with immune cells were investigated. Finally, DEGs between two risk groups were measured and the prognostic value of key DE-FRGs in ESCA was confirmed in other databases. Results: A prognostic model was constructed based on seven selected DEG-FRGs. TNM staging and CD8+ T cells were significantly correlated with high-risk groups. Low-risk groups exhibited more infiltrated M0 macrophages, an activation of type II interferon (IFN-γ) responses, and were found to be more suitable for immunotherapy. Seven key DE-FRGs with prognostic value were found to be considerably influenced by different chemotherapy drugs. Conclusion: A prognostic model based on seven DE-FRGs may efficiently predict patient prognosis and immunotherapy response, helping to develop individualized treatment strategies in ESCA.


Introduction
ESCA is the eighth most common cancer worldwide causing the sixth highest cancerrelated mortality rates [1]. Despite the range of available treatment options, the overall five-year survival rate remains less than 20%, with most deaths being associated with distant metastases and the emergence of resistance to chemoradiotherapy [2,3]. Therefore, exploration of novel therapeutic options and development of prognostic models for the management of ESCA are crucial.
The reprogramming of energy metabolism is a hallmark of cancer development, promoting cell growth and proliferation [4,5]. Accumulating evidence suggests that fatty acid metabolism plays a crucial role in metabolic reprogramming, affecting cell membrane formation, energy storage, and the production of signaling molecules [6,7]. Previous study has reported that fatty-acid-metabolism-related genes are associated with malignancy, prognosis, and immune phenotype in gliomas [8]. In cervical cancer patients, enhanced lipolysis and fatty acid synthesis promote lymphatic spread via the activation of nuclear factor kB (NF-kB) signaling [4,9]. Activated fatty acid oxidation improves the survival of acute myeloid leukemia cells [10]. In addition to influencing the effectiveness of chemoand radiation therapy, altered fatty acid metabolism has also been suggested to effect

Construction and Evaluation of a Predictive Risk Score Model
Univariate Cox regression analysis was performed to identify DE-FRGs associated with prognosis, using the "survival" and "survminer" packages in R. Additional LASSO regression analysis was conducted to narrow further the field of key DE-FRGs. Subsequently, multivariate Cox regression analysis was performed to develop a prognostic risk score model for predicting overall survival (OS) of ESCA patients. ESCA patients were divided into high-risk (n = 80) and low-risk (n = 80) groups based on the median risk score. Kaplan-Meier (KM) survival curves were plotted to analyze the difference in overall survival (OS) and progression free survival (PFS) of patients assigned into the high-risk and low-risk groups. Finally, the receiver operating characteristic (ROC) curve was drawn and the area under the curve (AUC) for 1-, 3-, and 5-year OS was calculated through the "survival ROC" package in R to assess the predictive accuracy of the prognostic risk score model. A p-value < 0.05 was used as the filter condition.

The Association between Risk Score and Clinical Characteristics
The Limma R package was utilized to explore the association between risk score and clinical characteristics, including age, gender, and TNM pathology stage [14]. This was followed by univariate and multivariate Cox regression analyses to identify independent prognostic factors using the "survival" package (p < 0.05) [16].

Construction of a Nomogram for Patients with ESCA
A nomogram consisting of age, gender, pathologic staging, and prognostic risk score model was constructed using the "survival", "regplot" and "rms" packages, to predict the likelihood of 1-, 3-, and 5-year survival of the studied ESCA patients. Then, independent prognostic analysis was performed to evaluate whether nomogram could independently predict patient prognosis. ROC curves and calibration curves were plotted to confirm the prediction accuracy of the nomogram [17].

Association between the Risk Model and Immune Parameters
The extent of immune cell infiltration was compared between the high-and low-risk groups using the CIBERSORT algorithms [18] and potential differences in immunologic functioning were evaluated using the "Limma", "GSVA", "GSEABase", "ggpubr", and "reshape2" R packages [19]. Finally, the tumor immune dysfunction and exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/login/ 18 March 2022) was applied to predict the effects of the observed changes on potential immunotherapy responses in the high-and low-risk groups [20].

GSCA Analysis
Gene Set Cancer Analysis (GSCA, http://bioinfo.life.hust.edu.cn/GSCA/#/ 20 June 2022) is an online analysis tool for genomic, pharmacogenomic, and immunogenomic gene expression analyses in cancer [21]. We used GSCA to analyze the correlations between the expression of key biomarkers, immune cells, and drug sensitivity.

A Protein-Protein Interaction Network of DEGs in Groups of Different Risk Score Groups
DEGs in 160 ESCA patients between two different risk groups were analyzed using the Limma R package (|logFC| > 1, FDR < 0.05) to identify DEGs [4]. Gene Ontology (GO) enrichment analyses of the detected DEGs was then performed using the "clusterProfiler" R package [22]. A protein-protein interaction (PPI) network of DEGs was then constructed and visualized via the STRING database (STRING, https://string-db.org/ 17 March 2022) and the Cytoscape software (version: 3.7.2), identifying hub genes using cytoHubba [4]. According to the median expression value of the hub genes, all samples were divided into low-and high-expression groups. Kaplan-Meier analysis was carried out to compare the survival characteristics between the two groups [23]. Finally, the correlation between immune cell infiltration and prognostically relevant DE-FRGs was analyzed [24].

Validation of the Expression and Prognostic Value of Seven FRGs
UALCAN (http://ualcan.path.uab.edu/index.html 28 June 2022) is a comprehensive interactive web resource for analyzing cancer OMICS data (TCGA, MET500, CPTAC, and CBTTC), allowing users to identify biomarkers or to perform in silico validation of potential genes of interest. It uses graphical representations of gene expression profiles of proteincoding, miRNA-coding, and lincRNA-coding genes and combines these with survival information [25]. UALCAN database was used to confirm the expression difference of key DE-FRGs and investigate whether their expression was correlated with survival differences between the two groups.

Identifying Fatty-Acid-Metabolism-Related DEGs in ESCA Samples
The database search-derived lists of fatty acid metabolism-related genes were intersected, and duplicates were eliminated. A total of 309 genes were identified as having a known or proposed role in fatty acid metabolism ( Figure 1A). Of these, 108 were differentially expressed between ESCA samples and healthy samples in the TCGA database (Figure 1B,C; Table 1).

Identifying Fatty-Acid-Metabolism-Related DEGs in ESCA Samples
The database search-derived lists of fatty acid metabolism-related genes were intersected, and duplicates were eliminated. A total of 309 genes were identified as having a known or proposed role in fatty acid metabolism ( Figure 1A). Of these, 108 were differentially expressed between ESCA samples and healthy samples in the TCGA database (Figure 1B,C; Table 1).

Establishing and Validating a Prognostic FRG Signature
Univariate Cox regression analysis indicated that 19 differentially expressed DE-FRGs (DE-FRGs) were associated with different clinical outcomes (p < 0.05) ( Figure 2A). LASSO analyses were used to further analyze 19 DE-FRGs to preventing overfitting. This additional analysis identified a signature, consisting of 11 genes, that showed altered expression levels depending on clinical prognosis ( Figure 2B,C). Subsequently, seven DE-FRGs were identified as potential prognostic-related biomarkers of ESCA patients using multivariate Cox regression analysis ( Figure 2D). Kaplan-Meier analysis indicated that compared with the low-risk group, the high-risk group was significantly associated with poor OS and PFS (p < 0.05) ( Figure 2E,F). Furthermore, time-dependent ROC curves were constructed and the corresponding area under the curve (AUCs) figures were calculated to assess the predictive power of the prognostic risk model. The AUC values for 1-, 3-, and 5-year OS were 0.787, 0.829, and 0.937, respectively ( Figure 2G).

Association between Risk Score and Clinical Characteristics
We analyzed the association between patient clinical characteristics and the gene expression-based risk score. This analysis showed no association between the age or gender of the patients and their corresponding risk scores ( Figure 3A,B). However, a higher risk score correlated with significantly more advanced pathologic stage and T stage of the tumor (p = 0.0017, Figure 3C, Figure S1). Univariate and multivariate Cox regression analysis revealed that our risk score and the clinical T stage could serve as independent prognostic factors (p < 0.001) ( Figure 3D,E). ROC curve analysis further supported the high sensitivity and specificity of the predictive score ( Figure 3F,G). acid metabolic genes; (C) The best parameter (lambda) in the LASSO-Cox model; the red dots indi cates the partial probability of deviation values, the gray lines indicates standard error (SE). (D) genes were eventually identified as prognosis related biomarkers based on multivariate Cox regres sion analysis; ESCA patients were divided into high-risk (n = 80) and low-risk (n = 80) groups based on the median risk score. *p < 0.05, **p < 0.01, ***p < 0.001. Black squares represent Hazard Rati (HR) value. (E,F) Kaplan-Meier survival analysis was performed to assess the difference in OS and PFS between the high-risk and low-risk group; (G) ROC curves of risk model for predicting overal survival at 1, 3, and 5 years.

Association between Risk Score and Clinical Characteristics
We analyzed the association between patient clinical characteristics and the gene ex pression-based risk score. This analysis showed no association between the age or gende of the patients and their corresponding risk scores ( Figure 3A,B). However, a higher risk score correlated with significantly more advanced pathologic stage and T stage of the tu mor (p = 0.0017, Figure 3C, Figure S1). Univariate and multivariate Cox regression analysi revealed that our risk score and the clinical T stage could serve as independent prognosti factors (p < 0.001) ( Figure 3D,E). ROC curve analysis further supported the high sensitivity and specificity of the predictive score ( Figure 3F,G).

Establishing and Evaluating a Nomogram for Predicting Survival
A nomogram, integrating age, gender, TNM stage of the tumor, and the risk score model was established to predict OS in ESCA patients ( Figure 4A). The calibration curves at 1 year, 3 years, and 5 years indicated that the nomogram could accurately predict the OS of ESCA patients ( Figure 4B). Based on univariate Cox regression analysis the nomogram model and the TNM tumor stage could predict the prognosis of the patients independently of each other or other clinical parameters ( Figure 4C). Multivariate Cox regression analysis showed that the nomogram model was an independent prognostic factor with an HR of 1.193 (95% CI = 1.060-1.343) ( Figure 4D). ROC analysis also revealed that the nomogram could predict the OS for ESCA patients with remarkable accuracy (AUC: 1 year = 0.736, 3 years = 0.849; Figure 4E,F).

Immunological Features of the Tumor and GSCA Analysis
The expression of seven key FRGs showed strong correlation with the infiltration of the tumors with immune-activating cells ( Figure 5A). The ratios of 22 distinct immune cell types in the tumors of high-and low-risk ESCA patients is shown in Figure 5B. The number of infiltrating CD8+ T cells was higher in the high-risk group, while M0 macrophages were more predominant in the low-risk patients ( Figure 5B). In terms of immune function, our results also showed that type II IFN production was significantly higher in the lowrisk group ( Figure 5C). Furthermore, significantly higher TIDE scores were seen in the low-risk than group, indicating that potentially the high-risk group could benefit more from receiving immunotherapy ( Figure 5D). We also explored the correlation between the

Immunological Features of the Tumor and GSCA Analysis
The expression of seven key FRGs showed strong correlation with the infiltration of the tumors with immune-activating cells ( Figure 5A). The ratios of 22 distinct immune cell types in the tumors of high-and low-risk ESCA patients is shown in Figure 5B. The number of infiltrating CD8+ T cells was higher in the high-risk group, while M0 macrophages were more predominant in the low-risk patients ( Figure 5B). In terms of immune function, our results also showed that type II IFN production was significantly higher in the lowrisk group ( Figure 5C). Furthermore, significantly higher TIDE scores were seen in the low-risk than group, indicating that potentially the high-risk group could benefit more from receiving immunotherapy ( Figure 5D). We also explored the correlation between the expression of the seven key FRGs and drug sensitivity, defined as IC50 values, based on Spearman's correlation analysis. Our results revealed that most drugs effected the seven key FRGs, suggesting that these molecules could be exploited as potential therapeutic drug targets in the management of ESCA ( Figure 5E,F).

DEGs in the Low-and High-Risk Score Groups
Using the "Limma" package, we identified 48 genes that were differentially expressed between the high-and low-risk groups ( Table 2). GO analyses indicated that these DEGs mainly belonged to the response to glucocorticoid, response to corticosteroid, intermediate filament cytoskeleton organization, intermediate filament-based process, skin development, negative regulation of peptidase activity, unsaturated fatty acid biosynthetic process, and regulation of peptidase activity KEGG pathways ( Figure 6A). A PPI network of 48 DEGs was constructed using the STRING database ( Figure 6B). The Top 10 genes of the network, including PPL, MMP9, TGM1, ALOX12, ANXA1, VIL1, IL1RN, GPA33, MLXIPL, and PCK1 of the network were selected using the cytoHubba plugin in Cytoscape ( Figure 6C). Analyzing these against the survival parameters showed that high expression of PPL was significantly associated with favorable OS of ESCA patients ( Figure 6D). The proportions of infiltrating immune cells in samples expressing PPL at high or low levels was analyzed using the "Limma" package. Our result showed that the presence of M0 macrophages was significantly more pronounced in samples with a high PPL high expression. In addition, the proportion of T cell CD4+ resting memory T cells was lower in the immune infiltrates in the high-PPL-expression group ( Figure 6E).

Validating the Expression of the Seven FRGs and Their Prognostic Value
The UALCAN online database was used to analyze how the expression of key FRGs affected survival times. The results indicated that FABP2, HSPH1, and IDH3G were more abundantly expressed in tumors ( Figure 7A-C). Furthermore, a higher abundance of PDHA1 correlated with poor prognosis in ESCA patients ( Figure 7D). Similarly, the increased expression of HSPH1, IDH3G, NUDT7 and PDHA1 was associated with shorter OS in the subgroup of patients suffering from esophageal adenocarcinoma (EAD) ( Figure  7E-H).

Validating the Expression of the Seven FRGs and Their Prognostic Value
The UALCAN online database was used to analyze how the expression of key FRGs affected survival times. The results indicated that FABP2, HSPH1, and IDH3G were more abundantly expressed in tumors ( Figure 7A-C). Furthermore, a higher abundance of PDHA1 correlated with poor prognosis in ESCA patients ( Figure 7D). Similarly, the increased expression of HSPH1, IDH3G, NUDT7 and PDHA1 was associated with shorter OS in the subgroup of patients suffering from esophageal adenocarcinoma (EAD) (Figure 7E-H).

Discussion
There is increasing evidence that metabolic dysregulation plays a critical role in ca cer cell growth, proliferation, angiogenesis, and invasiveness [14,26]. Previous observ tions suggest that abnormal glycolytic metabolism is associated with the physiologi behavior of human malignant tumors [27]. In colorectal cancer, abnormal anaerobic m abolic pathways play an important role in the formation of cancer stem-like cells (CSC promoting the rapid formation, development, and therapy resistance of these tumo

Discussion
There is increasing evidence that metabolic dysregulation plays a critical role in cancer cell growth, proliferation, angiogenesis, and invasiveness [14,26]. Previous observations suggest that abnormal glycolytic metabolism is associated with the physiological behavior of human malignant tumors [27]. In colorectal cancer, abnormal anaerobic metabolic path-ways play an important role in the formation of cancer stem-like cells (CSCs), promoting the rapid formation, development, and therapy resistance of these tumors [28,29]. Fatty acid metabolism is involved in cellular energy production, membrane synthesis, and signal transduction pathways relevant to tumorigenesis and development [14,30]. Deregulated anabolism/catabolism of fatty acids may support cancer cell growth [6]. A previous study reported that fatty acid synthesis may promotes esophageal adenocarcinoma [31]. A recent study has shown that loss of FBP1 promotes migration, proliferation and invasion through regulating fatty acid metabolism in ESCA [32].Although several studies focused on the role of fatty acid metabolism in a variety of tumors, this topic remains unclear in ESCA. The identification of key molecular markers related to fatty acid metabolism, together with the exploration of their role in the development of ESCA, could provide novel insights into the biological behavior of these tumors, potentially highlighting new, more effective, therapeutic strategies.
In the present study, we explore the role of DE-FRGs in ESCA. Interestingly, 19 FRGs that were differentially expressed in patients with ESCA showed a correlation with the clinical/biological behavior of the tumors. Via LASSO and multivariate Cox regression analysis, we were able to identify seven key DE-FRGs, including PDHA1, CD36, IDH3G, HSPH1, FABP2, NUDT7 and SERINC1 that may become useful prognostic biomarkers in the clinical management of ESCA patients. Subsequently, a prognostic risk score model was established that was able to divide patients into distinct high-and low-risk groups showing significant differences in OS and DFS. This prognostic risk score was an independent prognostic factor according to univariate and multivariate Cox regression analyses. Furthermore, the predictive potential of this model was confirmed by combination with clinical characteristics (age, gender, and TNM stage of the tumor) in a risk-assessment nomogram. The risk model presented here might help to identify ESCA patients with a poor prognoses, and could be utilized in the management of the disease. PDHA1(Pyruvate Dehydrogenase E1 Subunit Alpha 1) is a protein-coding gene involved in the pyruvate and thiamine metabolism pathways. In HNSCC cells, LDHA/PDHA1 changes may associated with a broad metabolic reprogramming while intracellular molecules including polyunsaturated fatty acids and nitrogen-metabolism-related metabolites underlie the malignant changes [33]. Previous study has shown that promotion of the tricarboxylic acid cycle (upregulation of PDHA1, PDHB, ACO2, and DLST expression) and inhibition of ME1 expression may inhibit fatty acid synthesis [34]. The inhibition of PDHA1 expression in the LnCap human prostate cancer cells led to the "Warburg effect", resistance to chemotherapy, improved migration, and increased expression of stem cell markers [35]. In ESCC, the low expression of PDHA1 correlates with poor clinical prognosis and can result in metabolic reprogramming, again leading to the Warburg effect increasing malignant potential [36,37]. Decreased PDHA1 protein expression was also found to predict poor prognosis in gastric cancer [38]. HSPH1 encodes a member of the heat shock protein 70 family. Previously, elevated levels of HSPH1 expression were reported in tissues of HNSC patients. Moreover, the overexpression of HSPH1 was associated with poor overall survival (OS) [39]. HSPH1 is one of the most prominently upregulated proteins in several malignancies, with a well-documented involvement in Wnt-and chronic nuclear factor-kappa B signaling [40][41][42]. It has been suggested that the analysis of HSPH1 expression levels may help in predicting the effectiveness of chemotherapeutic approaches acting on the EGFR-TKI pathway in advanced lung adenocarcinoma [40]. The involvement of HSPH1 in anticancer immunity has also been described [41]. NUDT7, Nudix Hydrolase 7, is an enzyme involved in peroxisomal lipid metabolism. In the liver, upregulated expression of NUDT7 can inhibit peroxisomal fatty acid oxidation [43]. In clear cell renal cell carcinomas the alternative splicing of NUDT7 was found to be correlated with overall survival time [44]. Fatty acid binding proteins (FABPs) are key proteins in lipid transport, which can maintain a steady pool of fatty acids in the epithelium by traffic lipids from the intestinal lumen to enterocytes and bind superfluous fatty acids. As a lipid chaperone, FABP2 can also carry lipophilic drugs to improve targeted transport [45]. During the work presented here we identified and validated the prognostic value of seven key DE-FRGs. Of these FABP2, HSPH1, and IDH3G were upregulated and the high expression of four FRGs, including PDHA1, HSPH1, IDH3G, and NUDT7 was associated with poorer OS. These findings suggest that these key DE-FRGs might have a prognostic value in the assessment of ESCA and could represent potential therapeutic targets.
Previous studies have suggested that prominent CD8+ T cell infiltrates were associated with clinical prognosis and immune responses in ESCA [46,47]. In renal cell cancer an increased CD8+ T-cells to Treg ratio is associated with poor prognosis [48]. We observed an increase in the number of infiltrating CD8+ T cells in the samples of high-risk group ESCA patients, indicating an unfavorable prognosis. In contrast, patients with a low-risk score showed an activated type II IFN response. Based on the TIDE algorithm this may indicate suitability for the immunotherapy. Previous studies demonstrated that type II IFN (IFN-γ) could indirectly regulate PD-L1 levels in small cell lung cancer [49]. It was also suggested that IFN-γ production induced elevated PD-L1-mRNA expression resulting in favorable OS [50]. Indeed IFN-γ production is a key driver of PD-L1 expression in both cancer and host cells, and may improve the likelihood of anti-PD-1 therapies being effective [51]. These previous findings support the notion that the favorable prognosis of patients in the low-risk group may be due to more effective immune responses.
It was reported that the expression of PPL was significantly decreased in esophageal cancer tissues and that the proteins were barely detectable in advanced cancer samples [52]. The experimental knockdown of PPL decreased cellular motility, reduced attachment, and generally inhibited malignant progression [53]. However, this finding contrasts with observations that the high expression of PPL is associated with favorable survival in patients with adenoid cystic carcinomas and sarcomas [54]. Our results showed that the expression of PPL gene was upregulated in low-risk group, and high expression of this gene was associated with better prognoses. This result may further support the hypothesis that the low-risk group was significantly associated with a favorable prognosis.

Conclusions
In summary, this first investigation of the role of fatty acid metabolism-related genes in ESCA identified seven genes showing strong association with the prognosis and therapeutic responses to chemotherapy and immunotherapy in this disease. This allowed us to develop a risk score model that could effectively divide patients into high-and low-risk groups. Using this score clinically might aid the development of individualized treatment strategies in the future. However, the prognosis value and immunotherapy effect of seven fatty-acidmetabolism-related genes in patients with esophageal cancer need to be further confirmed in larger clinical studies.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/vaccines10101721/s1, Figure S1: The association between risk score and clinical characteristics; Supplementary