Glycolysis Changes the Microenvironment and Therapeutic Response Under the Driver of Gene Mutation in Esophageal Adenocarcinoma

Background: Esophageal cancer is one of the most leading and lethal malignancies. Glycolysis and the tumor microenvironment (TME) are responsible for cancer progressions. We aimed to study the relationships between glycolysis, TME, and therapeutic response in esophageal adenocarcinoma (EAC). Materials and Methods: We used the ESTIMATE algorithm to divide EAC patients into ESTIMATE high and ESTIMATE low groups based on the gene expression data downloaded from TCGA. Weighted gene co-expression network analysis (WGCNA) and Gene Set Enrichment Analysis (GSEA) were performed to identify different glycolytic genes in the TME between the two groups. The prognostic gene signature for overall survival (OS) was established through Cox regression analysis. Impacts of glycolytic genes on immune cells were assessed and validated. Next, we conducted the glycolytic gene mutation analysis and drug therapeutic response analysis between the two groups. Finally, the GEO database was employed to validate the impact of glycolysis on TME in patients with EAC. Results: A total of 78 EAC patients with gene expression profiles and clinical information were included for analysis. Functional enrichment results showed that the genes between ESTIMATE high and ESTIMATE low groups (N = 39, respectively) were strongly related with glycolytic and ATP/ADP metabolic pathways. Patients in the low-risk group had probabilities to survive longer than those in the high-risk group (p < 0.001). Glycolytic genes had significant impacts on the components of immune cells in TME, especially on the T-cells and dendritic cells. In the high-risk group, the most common mutant genes were TP53 and TTN, and the most frequent mutation type was missense mutation. Glycolysis significantly influenced drug sensitivity, and high tumor mutation burden (TMB) was associated with better immunotherapeutic response. GEO results confirmed that glycolysis had significant impacts on immune cell contents in TME. Conclusion: We performed a comprehensive study of glycolysis and TME and demonstrated that glycolysis could influence the microenvironment and drug therapeutic response in EAC. Evaluation of the glycolysis pattern could help identify the individualized therapeutic regime.


INTRODUCTION
Esophageal cancer is the eighth most common malignancy and the sixth cause of cancer death globally, which accounts for more than 570,000 new cases and 500,000 deaths annually (Bray et al., 2018;Wang et al., 2018). Esophageal adenocarcinoma (EAC) is the predominant pathological type in western countries, with an increasing proportion from 35 to 61% over the past 30 years (Alsop and Sharma, 2016). The global incidence rate of EAC is approximately 0.7/100,000 person-years, and the 5-year survival rate is merely less than 20%, although multidisciplinary treatments have been applied, including esophagectomy, radiation, and chemotherapy (Arnold et al., 2015;Markar et al., 2017;Smyth et al., 2017;Zhao et al., 2019). Considering the chemotherapeutic resistance, several targeted agents have been applied in patients with EAC, such as imatinib (Mayr et al., 2012). Unfortunately, the efficacy is still not satisfactory. Recently, immunotherapy-targeting PD-1 has revolutionized the therapy in cancer patients. However, not all patients with EAC respond to immunotherapy (Däster et al., 2020). Therefore, there is an urgent necessity to better understand the molecular characteristics and genetic features that could help predict accurate survival and identify suitable patients who will benefit from immunotherapy.
The tumorigenesis and development of EAC is a highly complex biology, involving the tumor cell-intrinsic and cellextrinsic factors (Quante et al., 2018;Talukdar et al., 2018). Genetic alterations are the primary mechanisms that drive the initiation and progression of EAC, not only conferring tumor cells infinite proliferative abilities but also reprogramming metabolic pathways to adapt to the hostile environment, such as aerobic glycolysis (Hochwald and Zhang, 2017;Talukdar et al., 2018). The seminal discovery of tumor glycolysis has been considered a hallmark of cancer, proposed by Otto Warburg in 1923 (Warburg andMinami, 1923). The glycolytic phenotype renders cancer cells selective advantages by unlimited growth and attenuated apoptosis (Xu et al., 2017). In addition, it is gradually evident that elevated glycolysis is closely related to the immune escape by changing the microenvironment and inhibiting the functions of immune cells (Jiang et al., 2019a). Mounting evidence from cell-based assays has linked glycolysis to TME, and preclinical investigations have demonstrated the effectiveness targeting glycolysis in some cancers (Lim et al., 2017;Kornberg et al., 2018;Jiang et al., 2019a;Jiang et al., 2019b;Kang et al., 2020).
Genetic mutation results in the rewire of the glucose metabolism decreased cancer cell apoptosis and immortal growth. Consequently, these events bring about the component reconstruction in the tumor microenvironment (TME), thus changing the purity of the tumor. Reciprocally, the intimate interactions between glycolytic cells and the extracellular matrix further exacerbate the remodeling of TME, including the stromal and immune cells. It is well accepted that the tumor is highly dependent on TME, which is preponderant on prognosis and impacts the therapeutic efficacy profoundly, such as the immune checkpoint therapy (Wu and Dai, 2017;Alsina et al., 2018;Taube et al., 2018;Hinshaw and Shevde, 2019;Li et al., 2020). However, there are few studies exploring the associations between glycolysis and TME in EAC, and far less is known about how genetic mutations orchestrate the glycolysis under these aberrant TME conditions. Herein, we investigated the effects of glycolysis on immune cells and revealed the genetic mutation diversity. Our study unraveled that glycolysis could influence TME under the driver of genetic mutation and could serve as prognostic biomarkers. Moreover, we constructed the risk score system to predict drug sensitivity and immunotherapeutic response. The results hold great promise in targeting glycolysis and utilizing TME to improve the treatment in patients with EAC.

Data Acquisition
Gene expression data and clinical information were downloaded from the Cancer Genome Atlas (TCGA) database (https://portal. gdc.cancer.gov/). The mRNA expression profiles were log2 normalized for further analysis. Clinical information included gender, age, stage, survival status, and follow-up time.

Tumor Microenvironment and Glycolysis
TME is composed of resident stromal cells and infiltrating immune cells (IICs), reflecting tumor purity. With the increase of stromal cells and IICs, the tumor purity becomes lower. The stromal score, immune score, and ESTIMATE score were calculated by applying the ESTIMATE algorithm (Chakraborty and Hossain, 2018). The ESTIMATE score is the comprehensive parameter of the stromal and immune scores. Then, patients with EAC were classified into ESTIMATE high and ESTIMATE low groups according to the median of the ESTIMATE score. The differently expressed genes (DEGs) were screened by the weighted gene co-expression network analysis (WGCNA) with the false discovery rate (FDR) ≤0.05 and log2 fold change (log2FC)| >2.
To explore whether glycolysis affects tumor purity, we performed gene set enrichment analysis (GSEA) between the ESTIMATE high and ESTIMATE low groups. Five glycolysisrelated gene sets, namely, Hallmark, BioCarta, KEGG, GO, and Reactome, were downloaded from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb) and analyzed using the GSEA software (version 4.1.0). The permutation number was set as 1,000 for every phenotype. The gene sets were considered statistically significant when the nominal (NOM) p-value ≤0.05, FDR ≤0.05, and normalized enrichment score |(NES)| >1. Finally, the intersection genes (IGs) from the WGCNA and GSEA gene sets were identified for further analysis.

Gene Interaction Analysis and Enrichment Analysis
Gene interaction analysis was performed through the "corrplot" package in R software (version 4.1.0). Hub genes were screened with "cytoHubba" in Cytoscape software.
The IGs based on TME and glycolysis were analyzed for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) by the "clusterProfiler" R package. GO analysis has three functional parts, including the biology process (BP), cellular component (CC), and molecular function (MF).
Next, we performed functional similarity analysis, which was measured through the "GOSemSim" R package (Wang et al., 2007). Functional similarity could be used for the purpose of assessing the intimacy and relationship between each gene and its partners by evaluating the function and location.

Establishment of Prognostic Signatures
First, univariate Cox regression analysis was used to identify IGs which were related to patients' overall survival (OS). Then, statistically significant IGs (p < 0.05) were enrolled into the multivariate Cox regression. Finally, patients were divided into high-and low-risk groups according to the median of the risk score, in which the risk score was calculated as follows: risk score j n 1 Coefj p Xj, with Coef j indicating the coefficient and Xj representing the relative expression levels of each IG standardized by the z-score.

TME and Gene Mutation
Next, we selected the prognostic glycolysis-related genes and hub genes to investigate their relations with IICs between the two groups through single-sample gene set enrichment analysis (ssGSEA) using the "GSVA" R package. The effects of OGG on immune cells were assessed using the linear regression.
To analyze why glycolysis affects TME, we calculated the glycolytic gene mutation frequency, variant classification, variant type, and single nucleotide variants (SNVs) between the ESTIMATE high and ESTIMATE low groups. Additionally, to fully understand the role of gene mutation in TME, we performed tumor mutation burden (TMB) analyses and explored their relationships with IICs through the simple nucleotide variation data from TCGA and cBioPortal online databases (http://www.cbioportal.org/).

Drug Sensitivity Analysis and Immunotherapy Response
The drug sensitivity of each patient with EAC was predicted by the Genomics of Drug Sensitivity in Cancer database (GDSC; https://www.cancerrxgene.org/). The half-maximal inhibitory concentration (IC50) was calculated through the "pRRophetic" R package, and the IC50 differences between the high-and lowrisk groups were compared (Geeleher et al., 2014).
The response to immunotherapy was estimated using the Tumor Immune Dysfunction and Exclusion website (TIDE; http://tide.dfci.harvard.edu/login/). The TIDE and PDL-1 scores were compared between the high-and low-risk groups.

External Cohort Validation
The impacts of glycolysis on TME in EAC were validated through the Gene Expression Omnibus (GEO) database (https://www. ncbi.nlm.nih.gov/gds/). The study was considered eligible for external cohort validation according to the following criteria: 1) studies with Homo sapiens samples and 2) studies with a sample number more than 50, and 3) studies with detailed experiment information and complete expression profiles. The primary goal of validation was to confirm whether the ESTIMATE algorithm method is suitable for patients with EAC, and the secondary goals were to determine whether glycolysis could influence the components of the microenvironment and affect drug sensitivity. The overall design of this study is shown in Figure 1.

Statistical Analysis
All statistical analyses were performed by the R software (version 4.1.0). The DEG analysis between the ESTIMATE high and ESTIMATE low groups was carried out by applying the unpaired t-test. Cox regression analysis was used to determine the prognostic factors. Kaplan-Meier (K-M) curves and log-rank tests were utilized to assess the prognostic outcome. The Mann-Whitney U test was used to compare the immune score, immune cell infiltrations, and immune signatures. Spearman's correlation analysis was used to evaluate the interactions. p < 0.05 was considered significant.

Identification of Tumor Purity and Glycolysis
A total of 87 EAC samples and gene expression data were available from the TCGA database, including 9 normal and 78 Frontiers in Genetics | www.frontiersin.org December 2021 | Volume 12 | Article 743133 FIGURE 1 | Schematic of the study design. A total of 78 EAC patients were recruited for further study. TME is mainly composed of tumor cells, stromal cells, and immune cells. With higher stromal and immune cells, the tumor purity is low.   Table S2).
Then, GSEA was conducted to assess the glycolytic differences between the ESTIMATE high and ESTIMATE low groups. The results showed that the GO glycolytic process (NES −1.54, NOM p 0.050, FDR 0.050) and Reactome glycolysis (NES −1.83, NOM p 0.008, FDR 0.008) were significantly enriched in ESTIMATE high group patients ( Figures 2D,E). There were no significant enrichments in the Hallmark (NES −1.44, NOM p 0.067, FDR 0.067), BioCarta (NES −1.23, NOM p 0.227, FDR 0.227), and KEGG (NES −1.20, NOM p 0.228, FDR 0.228) pathways ( Figures 2F-H). There were 117 DEGs between GO and Reactome glycolysis gene sets. After screening, a total of 34 IGs were selected from WGCNA and GSEA for further analysis (Table 1) ( Figure 2I). The details about IGs in each EAC sample are shown in Supplementary Table S3.

Gene Interaction Networks and Functional Enrichment Analysis
To explore the correlation between the IGs, we calculated their coefficients. Gene interaction analysis showed that GAPDH and TPI1 had the strongest positive correlation (coef 0.81), whereas PRKACB and RAE1 had the strongest negative correlation (coef −0.48) ( Figure 3A). To explore the IG functions, we performed GO and KEGG enrichment analyses using R packages. GO results showed that IGs were significantly enriched in the glycolytic and ATP/ADP metabolic pathways. In addition, nuclear-, glucose-, and carbohydrate-related activities were closely associated with CC and MF terms ( Figure 3B). KEGG results demonstrated that carbon, gluconeogenesis, and HIF−1 signaling pathways were enriched ( Figure 3C).
Based on the GO analysis and semantic similarities, we ranked the genes by average functional similarities between IGs and their partners, with the cutoff value at 0.75. The box plots and raincloud plots are demonstrated in Figures 3D,E. From the pictures, we can clearly see that NUP43, NUP37, and DDIT4 had strong similarities and weak correlation with BPGM.

Prognostic IG Signatures
Univariate Cox regression analysis revealed that NUP88, RAE1, SEH1L, NUP37, and NUP43 were significantly associated with patients' OS (all p < 0.05) ( Figure 4A). After multivariate Cox regression analysis, three IGs (NUP88, SEH1L, and NUP37) were used to develop the risk score based on the following formula: risk score 0.637 * expression of NUP88 + 0.494 * expression of SEH1L + 0.657 * expression of NUP37. Also, the three genes were all risk genes with hazard ratio (HR) > 1. A total of 78 patients with EAC were classified into low-and high-risk groups according to the median risk score (n 39). The K-M survival plot showed that patients in the low-risk group had significant probabilities to survive longer than those in the high-risk group (median time 1.75 vs 0.745 years, p < 0.001) ( Figure 4B).
In order to evaluate the prognostic values of clinical information in OS, we integrated the patients' clinical features with IGs. Univariate Cox regression analysis showed that the tumor stage (HR 3.308, p < 0.001) and risk score (HR 1.954, p 0.013) were significantly associated with OS ( Figure 4C). Multivariate Cox regression analysis results demonstrated that tumor stage (HR 7.971, p < 0.001), metastasis (HR 0.167, p 0.033), and risk score (HR 2.822, p 0.002) were independent risk factors for OS ( Figure 4D). In addition, the distributions of each patient and their survival statuses are shown in Figures 4E-G. We can clearly see that patients in the low-risk group had a better prognosis than those in the high-risk group.

Effect of Glycolytic Genes on IICs
We selected three prognostic genes and five hub genes (NUP88, SEH1L, NUP37, GCK, NUP62, NUP155, NUP205, and NUP214) to assess whether glycolysis affects the IICs in TME. The results demonstrated that NUP62, NUP155, NUP205, and SEH1L had  Frontiers in Genetics | www.frontiersin.org December 2021 | Volume 12 | Article 743133 6 significant impacts on the IIC expression level, especially on the T-cells and mast cells (all p < 0.05). The details are shown in Figure 5A.
To fully explore the relationships between these genes and IICs, we performed Spearman correlation analysis by using the "Limma" package. The results showed that NUP62 was strongly associated with T-cells, dendritic cells, and antigen-presenting cells (all p < 0.05) ( Figure 5B). NUP155 had a close relationship with T-cells, mast cells, and MHC class I activity (all p < 0.05) ( Figure 5C). NUP205 showed close relationships with T-cells, mast cells, and type I IFN response (all p < 0.05) ( Figure 5D). SEH1L exhibited significant associations with dendritic cells, antigen-presenting cells, and chemokine receptors (CCR) (all p < 0.05) ( Figure 5E). Collectively, these findings suggested that the IGs had profound effects on immune cells and immunological functions.

Genetic Mutation and the Tumor Microenvironment
Gene mutations were analyzed for further investigation into the mechanisms that TME is affected by gene alteration. The genetic mutations are significantly different between the ESTIMATE high and ESTIMATE low groups. In the ESTIMATE high group, the five most common mutant genes were TP53, TTN, HMCN1, DNAH5, and SYNE1, and the most common mutational type was missense mutation (Figures 6A,C). In the ESTIMATE low group, the most common mutant genes were TP53, TTN, MUC16, SYNE1, and PCLO. The most common mutational type was also missense mutation ( Figures 6B,F). The frequencies of missense mutation were significantly lower than that in the ESTIMATE high group, indicating that the glycolytic level and tumor purity were different from those of the ESTIMATE high group (Figures 6A,B). Single-nucleotide polymorphism (SNP) had the highest frequency in the variant type ( Figures 6D,G). G > A was the most frequent type in the SNV class ( Figures 6E,H). The results imply that these mutant genes drive a higher glycolytic level, consequently changing the tumor purity in the microenvironment.
To broaden the understanding of the glycolytic gene mutations, the cBioPortal database was applied to validate these findings. We selected the most common glycolytic genes (TP53, TTN, HMCN1, DNAH5, SYNE1, MUC16, and PCLO) for further verification. Consistent with the above findings, the results from the cBioPortal database showed that TP53, TTN, and MUC16 possessed the highest mutation frequencies (87, 42, and 26%, respectively), and missense mutation was the commonest type ( Figure 7).

Tumor Mutation Burden Analysis and Immunotherapy Response
TMB refers to the total number of mutations per mega base in tumor tissues. By analyzing the SNP data downloaded from TCGA, we calculated the TMB frequency in each patient with EAC. The range of TMB is from 0.053 to 41.053 in EAC. Moreover, we further analyzed the effect of TMB on survival in patients with EAC. A total of 78 patients were classified into high-and low-TMB groups according to the median TMB. As shown in the K-M curves, patients in the high-TMB group had significantly higher mortality than those in the low-TMB group (p 0.05) ( Figure 8A). To exhibit the relationships between TMB, glycolysis, and TME (ESTIMATE score), we applied the Sankey diagram to visualize their correlations with the "ggalluvial" package in R. The result is shown in Figure 8B.
Low TMB usually implies poor response for immunotherapy (Chan et al., 2019). To explore whether TMB will influence the immunotherapy response, we compared the PDL-1 and TIDE scores between high-and low-TMB groups. As a result, patients in the high-TMB group significantly responded to anti-PDL-1 therapy (p 0.040) ( Figure 8C). In Figure 8D, patients with high TMB had a higher TIDE score than those with low TMB (p 0.000). In addition, the correlations between the immune checkpoints are also explored in Figure 8E.

Drug Sensitivity Analysis
We compared the IC50 differences of chemotherapeutic and targeted drugs between high-and low-risk score groups, including bexarotene ( Figure 9A), camptothecin ( Figure 9B), gemcitabine ( Figure 9C), imatinib ( Figure 9D), methotrexate ( Figure 9E), and vorinostat ( Figure 9F). The results demonstrated that there were higher IC50 levels of bexarotene and imatinib in the high-risk score group, which indicated that patients with a low-risk score were more sensitive to the two drugs. Oppositely, the IC50 levels of camptothecin, gemcitabine, methotrexate, and vorinostat were higher in the low-risk score group, implying that patients in the high-risk score group were more sensitive to the four drugs.

Validation of Glycolysis and Its Impact on TME
For external cohort validation, GSE12898 was employed, which consisted of 75 EAC samples. Consistent with the classification result of the ESTIMATE algorithm in TCGA, this method successfully divided 75 patients with EAC into the ESTIMATE high and ESTIMATE low groups (n 38 and 37, respectively). There were 3,428 genes between the two groups, including 52 significantly different glycolysis-related genes (Supplementary Table S4). The correlation analysis showed that 28 glycolytic genes had significant impacts on immune cells in the microenvironment, and the majority were B-cell and T-cell subtypes. The details are shown in Table 2. Furthermore, the drug sensitivities were also analyzed between the ESTIMATE high and ESTIMATE low groups. The results showed that the IC50 level of bexarotene ( Figure 10A) was higher in the ESTIMATE high group. However, the IC50 levels of camptothecin ( Figure 10B), gemcitabine ( Figure 10C), and vorinostat ( Figure 10F) were lower in the ESTIMATE high group, implying that the patients in the ESTIMATE high group were more sensitive to the four drugs. There were no significant differences regarding IC50 in imatinib ( Figure 10D) and methotrexate ( Figure 10E). Taken together, glycolysis directly changed TME and indirectly influenced drug sensitivity.

DISCUSSION
The abilities of cancer cells to switch metabolisms and evade the immunity system in TME are well-documented characteristics in tumors. Elevated glycolysis is commonly observed in cancer progression and is associated with significant disruptions of a previous finely tuned microenvironment (Chang et al., 2015;Jiang et al., 2019a). As the tumors develop, they constantly interact with neighboring cells, such as stromal cells and immune cells, under the driver of genetic mutations, thus altering their phenotypes and functions (Butturini et al., 2019). In the context of these intricate crosstalks between malignant cells and non-malignant cells, the impacts of increased glycolysis and a dysregulated TME on immune response and effective therapy are of vital importance (Roma-Rodrigues et al., 2019). Although the research studies focusing on the tumor glycolysis and TME have exploded exponentially in recent years, the underlying mechanisms of how they act both independently and synergistically still remain elusive. In this study, we explored systematically the links between glycolysis and TME in patients with EAC as well as the relations with genetic mutations. In the present study, we found that the glycolytic level was higher in the ESTIMATE high group, reasonably reflecting the fact that glycolysis may change the tumor purity. By constructing the predictive survival model based on IG signatures, we discovered that three genes (NUP88, SEH1L, and NUP37) may serve as independent prognostic biomarkers for OS in EAC. In addition, we revealed that gene mutation types and frequencies were distinct between the different ESTIMATE score groups, lending us a hypothesis that genetic alteration may drive TME changes. In addition, our data suggested that glycolysis could influence drug sensitivity and immunotherapeutic response. Several studies have confirmed the close relationships between glycolysis and EAC (Lynam-Lennon et al., 2014;Harada et al., 2020;Kang et al., 2020). Consistent with these findings, our functional enrichment analysis results showed that IGs were strongly enriched in the glycolytic processes, ATP generation, the HIF−1 signal pathway, RNA activities, and so on. The presence of aerobic glycolysis under normal conditions efficiently promotes tumor cell growth by the following mechanisms: 1) the rate of glycolysis is accelerated at 100 times compared to oxidative phosphorylation to compensate for the mathematical disadvantage in terms of net ATP production (glucose oxidative phosphorylation although mitochondria generate 18 times ATP compared to glycolysis) (Pfeiffer et al., 2001); 2) glycolysis provides sufficient and essential intermediates, such as NADPH and ribose-5-phosphate, which are indispensable for biosynthesis to meet avid proliferative requirements (Boroughs and DeBerardinis, 2015); 3) lactate, the obligatory product derived from glycolysis, could activate the HIF−1 signal pathway to induce vascular endothelial growth factor (VEGF) expression to stimulate angiogenesis in the microenvironment (Sonveaux et al., 2012). In addition, lactate is crucial to reorganize the tumor physical architectures in TME, and the accumulation of extracellular lactate is detrimental for normal healthy cells, such as immune cells (Romero-Garcia et al., 2016;Cassim and Pouyssegur, 2019). The study about how RNA interacts precisely with glycolysis is still in its infancy. However, Hua Q et al. gave us a hint that long non-coding RNA may promote glycolysis by sponging miRNA (Hua et al., 2019). Hence, it is reasonable to speculate that certain mutational genes regulated and modified by RNA may be involved in the glycolysis in EAC (Hochwald and Zhang, 2017). The prognostic signature of glycolysis in EAC was established based on three genes (NUP88, SEH1L, and NUP37). NUP88, located at chromosome 17p13, encodes Nup88 protein (Zhao et al., 2012). Nup88 is a nucleoporin comprising nuclear pore complexes (NPCs) and plays critical roles in maintaining the spindle stability and preventing aneuploidy formation during mitosis (Hashizume et al., 2010). Unanimously, the GO results in our study also illustrated that glycolytic genes had intimacy with the nuclear pore, emphasizing the importance of nuclear proteins in glycolysis. The Nup88 overexpression is highly associated with tumor development and decreased survival, suggesting that NUP88 acts as an oncogene (Martínez et al., 1999;Naylor et al., 2016). In line with these studies, our results also proved that NUP88 was a risk factor for OS (HR > 1). Another prognostic gene is NUP37, which shares similarities with NUP88 and is also a member of NPC. It exerts primary functions of sustaining NPC integrity and modulates the cell cycle . Previous studies showed that the elevated expression of NUP37 is associated with worsened survival rates in liver cancer (Uhlen et al., 2017). In addition, a latest study by Huang L et al. demonstrates that NUP37 silencing induces inhibition of lung cancer cell proliferation (Huang et al., 2020). These findings were in agreement with our results, pointing out that NUP37 played an oncogenic role in OS (HR > 1). Nonetheless, the function of NUP37 in EAC has never been explored and needs further experiments to confirm in vitro and in vivo. SEH1L, also known as Seh1, is a part of NPC as well. However, the field is still in its infancy, and only a handful of animal models have been developed to investigate the role of SEH1L. Studies have shown that Seh1 could promote oligodendrocyte differentiation Raices and D'Angelo, 2019). However, little is known about how it works in EAC. Undeniably, more research studies are warranted to understand the specific effects of SEH1L on EAC. The notion that glycolysis has profound impacts on immune cells in TME is well recognized. The investigation linking metabolic demands and immune cells was first documented in neutrophils and macrophages (Alonso and Nungester, 1956;Newsholme et al., 1986). Immune cells hold a resemblance with tumor cells to engage glycolysis, which require rapid energy sources to produce immune-mediators for migration and phagocytosis. Our study demonstrated that NUP88 had significant correlation with B-cells and mast cells (all p < 0.05) and NUP37 had positive correlation with T-cells and negative correlation with dendritic cells (DCs) and macrophages (all p < 0.05). These findings enforced the concept that glycolysis contributes to dramatic alterations of immune cells in TME, hence influencing the immune response and immune-based treatments. Glycolytic tumor cells compete with T-cells for glucose and impair T-cell activation (Peng et al., 2016). In addition, IFN-γ secreted by Th1 cells is sensitive and is unstable to lactate. Noteworthily, this could polarize the T-cells toward differentiating into the Th2 subpopulation that favors tumor progression by inhibiting the antitumor effect of M1 macrophages (Kareva and Hahnfeldt, 2013). Paradoxically, direct evidence from in vivo experiments supports that glycolysis could promote Th1 cell differentiation through an epigenetic mechanism (Peng et al., 2016). This calls for further studies aiming to shed light on the characteristics of glycolysis on T-cells. The divergent effects of glycolysis on DC are also observed. Glycolysis produces excessive lactate and lowers the pH in TME. Lactate together with decreased pH suppresses DC differentiation and consequently abolishes the antigenpresenting functions to T-cells (Harmon et al., 2020). However, inhibition of glycolysis also blocks DC maturation through HIF−1α. This highlights the multiple roles of DC in glycolysis and needs to be carefully interpreted in a contextdependent manner.
Gene mutations facilitate tumor cell metabolic plasticity to create favorable microenvironments beneficial to uncontrolled proliferation. Deepening our knowledge about the differences of gene alterations between different TMEs may allow for the development of therapeutic strategies. Spurred by this promising target, we explored upon this issue and showed that ESTIMATE high and ESTIMATE low groups manifested different gene mutation profiles, in which TP53 and TTN were the most prevalent mutant genes. Moreover, the most common type is missense mutation. Mutant TP53 endows tumor cells with adaptabilities to cope with the harsh microenvironment by providing adequate nutrients, thereby escaping from antitumor immune attack. Experimental mice models in the study by Basu S et al. testified that mutant TP53 could rewire the tumor glycolytic metabolism and enhance metastasis in TME (Basu et al., 2018). Mutant TP53 has additional impacts on TME beyond changing  (Kieser et al., 1994). Second, mutant TP53 could regulate chemokines that are involved in the homeostatic microenvironment (Yeudall et al., 2012). Last but not the least, mutant TP53 could reprogram the infiltrating immune cells and reshape the microenvironment (Cooks et al., 2018). TTN is located on chromosome 2q31, consisting of 364 exons, and it is the longest described coding gene (Chauveau et al., 2014). TTN ranks the third most abundant protein in both cardiac and skeletal muscle tissues, followed by actin and myosin (Chauveau et al., 2014). However, its mutation is not a rare entity in various cancers. Consistent with prior investigation by Cheng X et al., we also found that TTN missense mutation was a rather frequent type (Cheng et al., 2019). In addition, TTN and TP53 comutation is often accompanied during tumorigenesis and may serve as a prognostic biomarker either alone or in combination [61][62][63] . Current guidelines recommend adjuvant chemotherapy for patients at an advanced stage. However, how to select suitable patients who will benefit more from the chemotherapeutic regime is the prior concern. Our data demonstrated that patients with a high-risk score were more sensitive to camptothecin, gemcitabine, methotrexate, and vorinostat, suggesting that targeting glycolysis may alleviate the chemotherapeutic resistance. Despite immunotherapy bringing about a breakthrough for cancer patients, only a minority of patients could reap survival benefits actually. The TMB analysis in our study will accurately and effectively identify which patients will respond to immunotherapy in patients with EAC. Collectively, the proposed risk score system in our study has potency to help clinicians devise an individualized treatment strategy.
The strength of the present study is such that we performed a systematic analysis about glycolysis and TME in EAC for the first time based on the National Public Database, which provides robust data and statistical support. This study draws a close link between tumor glycolysis and the microenvironment and tentatively explains the mechanisms from the viewpoint of genetic alteration. Meanwhile, there are several limitations. First of all, TME is a complex mixture of parenchymal cells, the extracellular matrix, and numerous cytokines except for tumor and immune cells. These are not available from the public database and may greatly affect the analysis. Second, the results are not validated in in vitro and in vivo experiments. Last, the methods proposed in this study may not be applicable to all tumors as a result of heterogeneity. Notwithstanding its limitations, our study does provide an overview of glycolysis and TME in EAC, and this lays the foundation for further basic research in the area of metabolism and the microenvironment.
In summary, we found that glycolysis could change the microenvironment under the driver of genetic mutation and influence the immunotherapy in EAC. New efforts target that EAC should incorporate the idea that the glycolytic metabolism could reshape TME. Further studies are necessary to confirm our conclusion.

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 authors.

AUTHOR CONTRIBUTIONS
LZ contributed to conception, design, and data acquisition of the study. FY contributed to interpretation and data analysis. XL contributed to data collection and interpreted the results. QL conceived and designed the study protocol. CZ reviewed and approved the final version of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Natural Science Foundation of China (NSFC, No. 81773266) and the Key Discipline Group Construction Project of Pudong Health Bureau of Shanghai, China (No. PWZxq 2017-13) to QL.