Comprehensive Analysis of Hexokinase 2 Immune Infiltrates and m6A Related Genes in Human Esophageal Carcinoma

Background: Hexokinase 2 not only plays a role in physiological function of human normal tissues and organs, but also plays a vital role in the process of glycolysis of tumor cells. However, there are few comprehensive studies on HK2 in esophageal carcinoma (ESCA) needs further study. Methods: Oncomine, Tumor Immune Estimation Resource (TIMER), The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) database were used to analyze the expression differences of HK2 in Pan-cancer and ESCA cohort, and to analyze the correlation between HK2 expression level and clinicopathological features of TCGA ESCA samples. GO/KEGG, GGI, and PPI analysis of HK2 was performed using R software, LinkedOmics, GeneMANIA and STRING online tools. The correlation between HK2 and ESCA immune infiltration was analyzed TIMER and TCGA ESCA cohort. The correlation between HK2 expression level and m6A modification of ESCA was analyzed by utilizing TCGA ESCA cohort. Results: HK2 is highly expressed in a variety of tumors, and its high expression level in ESCA is closely related to the weight, cancer stages, tumor histology and tumor grade of ESCA. The analysis results of GO/KEGG showed that HK2 was closely related to cell adhesion molecule binding, cell-cell junction, ameboidal-type cell migration, insulin signaling pathway, hif-1 signaling pathway, and insulin resistance. GGI showed that HK2 associated genes were mainly involved in the glycolytic pathway. PPI showed that HK2 was closely related to HK1, GPI, and HK3, all of which played an important role in tumor proliferation. The analysis results of TIMER and TCGA ESCA cohort indicated that the HK2 expression level was related to the infiltration of various immune cells. TCGA ESCA cohort analyze indicated that the HK2 expression level was correlated with m6A modification genes. Conclusion: HK2 is associated with tumor immune infiltration and m6A modification of ESCA, and can be used as a potential biological target for diagnosis and therapy of ESCA.


INTRODUCTION
Recent studies show that Esophageal carcinoma (ESCA) ranks seventh in terms of incidence and sixth in mortality overall (Sung et al., 2021). Despite substantial improvements in the diagnosis and treatment of esophageal diseases, the prognosis for patients with ESCA remains poor (Bray et al., 2018). The occurrence and development of ESCA is an extremely complex biological process, and the expression of many genes has changed. Therefore, further research on the molecular mechanism of ESCA can provide new theoretical value for the diagnosis and treatment of tumors.
Hexokinase (HK) is an enzyme capable of phosphorylating hexose and is a rate-limiting enzyme in the glycolytic pathway. Four subtypes of human HK have been discovered, which are encoded by the genes HK1, HK2, HK3, and HK4, respectively (Zhang et al., 2017). HK1 is widely expressed in mammalian tissues, HK2 is usually expressed in insulin-sensitive tissues such as fat, bone and heart muscle, and HK3 is expressed at a low level, while HK4 expression is limited to pancreas and liver (Smith, 2000;Zhang et al., 2017). HK2 not only plays a role in physiological function of human normal tissues and organs, but also plays an important role in the glycolysis of tumor cells . Currently, researchers have found that expression of HK2 is increased in a variety of tumors and promotes the development of tumors (Shen et al., 2020;Wang et al., 2020). Our previous study found that HK2 was highly expressed in ESCA, but no more studies were conducted on the biological function of HK2 (Liu et al., 2020b).
Tumor immunotherapy and N6-methyladenosine (m6A) modification are hot spots in tumor treatment, which are extensively used in the investigation and therapy of ESCA (Baba et al., 2020;Wu S. et al., 2020;Wu X. et al., 2020;Yang et al., 2020). However, there are less investigates on the overall understanding of HK2 in ESCA, particularly the correlation between HK2 and ESCA immune cell infiltration and m6A modification.
On this project, we processed The Cancer Genome Atlas (TCGA) ESCA cohort and performed bioinformatics analysis using the R software, online website and other databases. The differences of HK2 expression in different cancer tissues were studied, and the mRNA and protein expression of HK2 in ESCA were verified by cell assay and immunohistochemistry (IHC). The co-expression gene networks of HK2 in ESCA were analyzed, and the possible biological mechanisms and signal pathways involved in associated genes were analyzed. Eventually, the correlation between the expression difference of HK2 and the tumor immune cell infiltration and m6A modification of ESCA will be explored, which will help to study the potential pathogenesis of ESCA.

Ethics Statement
The protocol of this study had been approved by the Ethics Committee of Taihe Hospital Affiliated of Hubei University of Medicine (Shiyan, China) (document NO.2021KS021) and conducted according to the principles stated in the Declaration of Helsinki, and the requirement to obtain informed consent was waived.

Expression of Hexokinase 2 in Pan-Cancer and Esophageal Carcinoma
Oncomine 1 (Rhodes et al., 2004(Rhodes et al., , 2007 and TIMER 2 (Li et al., 2016 databases were used to analyze the expression level of HK2 in different tumors. Oncomine database used Student's t-test to compare the transcription level of HK2 in clinical cancer samples and normal control group, and selected data with multiple change > 2 and P-value < 0.01. We also downloaded the RNA sequencing data of ESCA from TCGA 3 (Tomczak et al., 2015) and GEO (GSE38129) 4 (n = 60; GSE23400, n = 106) cohort to analyze the difference of HK2 expression between ESCA and normal tissues. UALCAN database 5 (Chandrashekar et al., 2017) was used to analyze the relationship between the expression level of HK2 and the clinicopathological characteristics of ESCA patients. Eventually, we validated the mRNA and protein expression of HK2 in ESCA and control sample by qRT-PCR and IHC assay according to the method described previously (Liu et al., 2020a(Liu et al., , 2021a, as detailed in Supplementary Materials.

Enrichment Analysis of Hexokinase 2 Gene Co-expression Network
The co-expression network of HK2 in TCGA ESCA cohort was analyzed using LinkedOmics 6 database. Pearson correlation coefficient was used for statistical analysis, and volcanic maps and heat maps were used for display. The rank criterion was an FDR < 0.05. The ClusterProfiler software package of R was used to analyze the Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway that co-expressed genes may participate in, and the ggplot2 software package was used to visually analyze the data.

Analysis of Gene-Gene Interaction and Protein-Protein Interaction of Hexokinase 2 Gene
We use GeneMANIA database 7 (Warde-Farley et al., 2010) to query and generate a list of genes that have similar functions with the HK2 gene and constructed an interactive network to illustrate the relationship between genes. The PPI of HK2 protein was calculated and predicted by using the STRING database 8 (Szklarczyk et al., 2019). The TIMER ESCA cohort shows that HK2 mRNA expression levels in different tumor types. (C) Differential expression of HK2 in ESCA and normal tissues was analyzed by TCGA ESCA cohort (D) Difference expression of HK2 in ESCA and normal tissues was analyzed by GSE38129 cohort. (E) Difference expression of HK2 between ESCA and matched normal tissues in GSE23400 cohort. (F) Difference expression of HK2 in different ESCA cell lines and human normal epithelial cell lines. Immunohistochemical staining showed the expression of HK2 in ESCA samples (G) and paracarcinoma samples (H). (I) The mean HK2 IHC score in ESCA tissue (2.32 ± 0.621) was significantly higher than that of matched peritumoral tissue (0.84 ± 0.618). * P < 0.05; * * P < 0.01; * * * P < 0.001; * * * * P < 0.0001.

Correlation Analysis of Hexokinase 2 and Immune Infiltrating Cells
The correlation between HK2 and immune infiltrating cells in ESCA samples was evaluated using the TIMER database. Immune infiltrating cells include B cell, neutrophil, CD4 + T cell, macrophage, CD8 + T cell, and dendritic cell. The somatic copy number alteration (SCNA) module of the TIMER tool was used to associate the genetic copy number variation (CNV) of HK2 with the relative abundance of tumor infiltrating cells. The CIBERSORT (Newman et al., 2015) software package of R was utilized to analyze the differences of 22 immune cells between high and low HK2 expression groups in ESCA samples. In addition, we also analyzed the correlation between HK2 and immune cell markers in ESCA samples utilizing TIMER databases and TCGA ESCA cohort. Immunofiltration marker genes refer to the previous study (Liu et al., 2021b).

Correlations of Hexokinase 2 Expression With m6A Associated Genes in Esophageal Carcinoma
The R software package was utilized to evaluate the correlation between the expression of HK2 and the expression of m6A associated gene in the TCGA ESCA cohort, including ZC3H13, YTHDF3, HNRNPC, METTL14, HNRNPA2B1,  . The R software was utilized to analyze the differences of 20 m6A associated genes between high and low HK2 expression groups in ESCA samples. The data were analyzed visually by ggplot2 software package.

Transcriptional Levels of Hexokinase 2 in Esophageal Carcinoma Patients
The ESCA cohort of TCGA and GEO was utilized to analyze the differential expression of HK2 in ESCA tissue samples and normal tissue samples. Both TCGA and GEO cohort analysis indicated that the HK2 expression level in ESCA tissue samples was remarkable higher than that in normal tissue samples (Figures 1C-E). To further verified the accuracy of data investigation, we conducted qRT-PCR and IHC experiments, respectively. The qRT-PCR results indicated that the HK2 mRNA expression in the two ESCA cell lines (ECA109 and KYSE-150) was remarkable higher than that in the normal human esophageal epithelial cells (Het-1A) ( Figure 1F). The results IHC staining revealed that HK2 was mainly expressed in the cytoplasm. The HK2 IHC score in tumor tissue samples was remarkable higher than that in paracarcinoma tissue samples (2.32 ± 0.621 vs. 0.84 ± 0.618) (Figures 1G-I, P < 0.0001). These results believe that HK2 plays a potential role in the occurrence and development of ESCA.
We analyzed clinical data from the ESCA cohort using the UALCAN database to better understand the clinical correlation between HK2 expression and ESCA. The results indicated that the expression of HK2 in cancer stage 2 was higher than that in cancer stage 3 (P < 0.05). The expression of HK2 in normal weight patients was higher than that in extreme weight patients (P < 0.05). The expression of HK2 in esophageal squamous cell carcinoma was higher than that in esophageal adenocarcinoma (P < 0.05). The expression of HK2 in tumor grade 3 was lower than that in tumor grade 2 (P < 0.01) and 1 (P < 0.05), respectively. However, the expression level of HK2 did not differ in different gender and age groups (Figure 2).

Enrichment Analysis of Hexokinase 2 Gene Co-expression Network
The co-expressed genes related to HK2 expression in TCGA ESCA cohort were analyzed using LinkedOmics database. The analysis results are shown in Figure 3A, 1682 genes were positively associated with the HK2 expression, and 2,855 genes were negatively associated with the HK2 expression (FDR < 0.05). The heat map indicated the top 50 important genes that are positively and negatively associated with HK2 expression, respectively ( Figure 3B). Supplementary Table 1 shows all co-expressed genes.
GO function and KEGG pathway research of HK2 coexpressed genes were carried by R language. Under the condition of p.adj < 0.1, HK2 co-expressed genes were participated in 25 biological processes (BP), 6 cell component (CC), 19 molecular function (MF) and 6 KEGG. The bubble graph shows the top 15 messages of GO and KEGG, including 5 messages of BP, CC, and MF. The results of GO functional annotations showed that HK2 co-expressed genes were mainly participated in the cell adhesion molecule binding, cell-cell junction, and ameboidal-type cell migration ( Figure 3C). KEGG pathway analysis showed that these genes were mainly related to the insulin signaling pathway, hif-1 signaling pathway, and insulin resistance ( Figure 3D). Supplementary Table 2 shows all GO function and KEGG pathway research.

Analysis of Gene-Gene Interaction and Protein-Protein Interaction of Hexokinase 2 Gene
A GGI network composed of 21 genes was constructed, and the HK2 gene was surrounded by 20 nodes. Their functionality was analyzed using the GeneMANIA database ( Figure 4A). These nodes represent genes closely related to HK2 in terms of physical interactions, shared protein domains, predicted, co-localization, pathway, coexpression, and genetic interactions. Analysis showed that there were five genes most related to HK2, which were TIGAR (TP53 induced glycolysis regulatory phosphatase), PTGDS (prostaglandin D2 synthase), NUDCD3 (NudC domain containing 3), FNTA (farnesyltransferase) and HKDC1 (hexokinase domain containing 1). Among them, TIGAR, PTGDS, NUDCD3 and FNTA are related to the physical interaction with HK2, respectively. Both HKDC1 and HK2 have predicted function and shared protein domains. The results also showed that these genes were most strongly associated with glycolysis (FDR = 4.53E-08). In addition, these genes are associated with carbohydrate catabolic process, carbohydrate kinase activity, single-organism carbohydrate catabolic process, cellular glucose homeostasis, glucose catabolic process and glucose homeostasis.

Correlation Analysis Between Hexokinase 2 Expression and Immune Marker Sets
To explore the relationship between HK2 and various immune infiltrating cells of ESCA, we used TIMER databases and TCGA ESCA cohort to study the associated between HK2 and immune marker genes of various immune cells ( Table 2). Both analyses indicated that the expression of HK2 was associated with the immune marker genes of CD8 + T Cell, Th1, T cell exhaustion and dendritic cell. The scatter plots showed the correlation between HK2 expression level and these three immune marker genes, respectively (Figure 7).

Hexokinase 2 Expression Is Correlated With m6A RNA Methylation Regulators in Esophageal Carcinoma
The TCGA ESCA cohort was analyzed to study the association between the expression of HK2 and the expression of 20 m6A related genes in ESCA. Analysis indicated that the expression of HK2 was remarkable positively associated with 4 m6A associated genes in ESCA, including IGF2BP2 (r = 0.200, P = 0.010), HNRNPA2B1 (r = 0.170, P = 0.035), YTHDC2 (r = 0.220, P = 0.004) and YTHDF2 (r = 0.200, P = 0.010) (Figure 8A). The scatter plot indicates the correlation between HK2 and m6A associated genes (Figure 8B). In addition, 162 cases of ESCA were divided into high expression group (n = 81) and low expression group (n = 81) according to the expression level of HK2. We further evaluated the expression abundance of m6A associated genes between high and low HK2 expression groups to study whether m6A modification was different between the two groups ( Figure 8C). Analysis indicated that compared with the low expression group, the expression of YTHDF3 and YTHDF2 were increased in the high expression group of HK2 (P < 0.05). These analyses suggest that HK2 is strongly associated to m6A modification in ESCA.

DISCUSSION
As one of the isozyme subtypes of HK, HK2 is the first ratelimiting enzyme in the catalytic glycolysis pathway and mainly distributed in the cytoplasmic region (Smith, 2000;Zhang et al., 2017). HK2 overexpression can promote the process of glycolysis of tumor cells and provide necessary energy for the proliferation and migration of tumor cells, thus promoting the occurrence and development of tumor cells (Shen et al., 2020;Wang et al., 2020;Yu et al., 2021). Chen et al. (2019) discovered that high expression of HK2 was remarkable associated with the degree of malignancy and poor prognosis of gallbladder tumors. Downregulation of HK2 expression could significantly inhibit the proliferation, migration and invasion of gallbladder cancer cells, and at the same time reduce glucose consumption and cellular lactic acid production. DeWaal et al. (2018) discovered that HK2 knockout inhibited the glycolytic effect of HCC cell, promoted the level of oxidative phosphorylation, and at the same time increased sensitivity of cancer cells to metformin. These results suggest that HK2 can be used as a potential target for cancer gene therapy.
Studies have shown that HK2 is highly expressed in a variety of tumors and plays an important role in the development and progression of tumors. Nevertheless, there are few researches on the comprehensive investigation of HK2 in ESCA. On the project, we predicted the expression differences of HK2 in cancers through bioinformatics investigation, and verified the HK2 expression in ESCA through in vitro experiments. Bioinformatics is widely used in tumor research Hu et al., 2021;Liao et al., 2021). Investigation of the Oncomine database discovered that HK2 was highly expressed in 8 types of cancer, and investigation of the TIMER discovered that HK2 was highly expressed in 11 types of cancer. According to GEO and TCGA ESCA cohort analysis, the expression level of HK2 in ESCA tissue samples was remarkable higher than that in normal tissue samples. We also detected the expression of HK2 in ESCA sample and normal sample by qRT-PCR and IHC. The investigation results were consistent with the above analysis. At present, we also found that the expression of HK2 was associated to weight, cancer stages, tumor histology and tumor grade. In summary, HK2 may serve as a possible diagnostic and therapeutic biomarker of ESCA. Nevertheless, previous researches on the role of HK2 in ESCA are limited to the energy metabolism of glycolytic pathway. Liu et al. (2019) found that the stem cell characteristics and abnormal metabolic reprogramming of esophageal cancer stem cells were dependent on the Hsp27-AKT-HK2 pathway, and ESCA patients with overexpression of Hsp27 and HK2 in tumor tissues had the worst prognosis. Mao et al. (2018) found that inhibiting the expression of HK2 could significantly inhibit the glycolytic effect of ESCA and affect the proliferation of cancer cell. Nevertheless, there are few reports on other biological functions of HK2 in ESCA. In this study, the co-expressed genes related to HK2 expression in TCGA ESCA cohort were analyzed using LinkedOmics database. The GO and KEGG function enrichment investigation of 100 co-expressed genes associated to the expression of HK2 indicated that the co-expression of HK2 was significantly correlated to the cell adhesion molecule binding, cell-cell junction, and ameboidal-type cell migration. KEGG pathway investigation indicated that HK2 co-expression was significantly correlated to the insulin signaling pathway, hif-1 signaling pathway, and insulin resistance. All these biological functions and pathways are associated to the occurrence and development of ESCA. GGI investigation discovered 20 genes related to HK2, and these genes were mainly closely related to glycolysis. The three genes most related to HK2 are TIGAR, PTGDS, and NUDCD3, respectively. Studies have shown that inhibition of TIGAR expression can reduce the proliferation of ESCA tumors, and patients who had high expression of TIGAR had poorer prognosis (Chu et al., 2020). PPI analysis showed that HK2 had the strongest correlation with the proteins of HK1, GPI and HK3, and inhibition of the expression of HK1 could significantly inhibit the proliferation of ESCA cells (Li et al., 2014). HK1, HK3, and GPI are all key genes in glycolysis, which can promote the proliferation and development of tumor cells by promoting the glycolysis effect of tumor cells. We discovered that the above three proteins had the highest association with HK2. These researches may suggest an association between HK2 and glycolysis and proliferative development of ESCA. These researches indicated that HK2 is not only involved in the glycolysis of ESCA, but also may play a variety of biological functions in the development and development of ESCA.
Tumor cell infiltration immunity has been shown to be associated with tumor progression and prognosis of ESCA (Baba et al., 2020;Wu X. et al., 2020). We used TIMER database and TCGA cohort to analyze the association between HK2 and ESCA immune infiltrating cells. TIMER database analysis showed that HK2 expression was correlated with B cell, CD4 + T cell and macrophage. In addition, the relationship between six coexpressed genes of HK2 and immune infiltrating cells was also analyzed. The co-expressed genes were MXD1, MAPK6, TGFA, LOC728743, RASSF4, and PBXIP1. We found that the expression of these six genes was also closely related to the infiltration of ESCA immune cells, which were all related to the infiltration of B cell, CD4 + T cell and macrophage. These results suggest that HK2 and its co-expressed genes may be participated in the immune response of ESCA tumor microenvironment, especially on B cell, CD4 + T cell and macrophage. In addition, it was also discovered that HK2 CNV has a remarkable correlation with the infiltration level of B cell, CD4 + T cell, and dendritic cell. According to the differential expression of HK2, the proportion of 22 tumor immune cells in ESCA was assessed by CIBERSORT research. We discovered five types of immune cells, including regulatory T cells, resting NK cells, M0 macrophages, activated mast cells and resting mast cell, whose proportions varied significantly according to the expression level of HK2. Furthermore, through the analysis of TIMER database and TCGA ESCA cohort, we discovered that the expression of HK2 was remarkable associated with the gene markers of CD8 + T Cell, Th1, T cell exhaustion and dendritic cell. We believe that high expression of HK2 in ESCA FIGURE 8 | Association of HK2 expression with m6A associated genes in esophageal carcinoma (ESCA). (A) The association between the expression level of HK2 and the expression of m6A associated genes in TCGA ESCA cohort was analyzed. (B) The association between HK2 and m6A associated genes was displayed by scatter plot, include IGF2BP2, HNRNPA2B1, YTHDC2, and YTHDF2. (C) The differential expression of m6A associated genes in the high and low HK2 expression groups was analyzed. *P < 0.05; **P < 0.01; ns, no significance.
patients may trigger an immune response. These results suggest that HK2 plays a remarkable role in the immune regulation of ESCA. Nevertheless, more studies are needed to further certify our speculation.
As part of epigenetics research, m6A is the most common and plentiful modification in post-transcriptional modification of RNA, which can affect the progress of tumor by regulating the biological functions associated with tumor. Shen et al. (2020) found that m6A associated gene METTL3 stabilizes the expression of HK2 and SLC2A1 in colorectal cancer through the m6A-IGF2BP2/3 dependent mechanism, and enhances the glycolysis ability of tumor cells, thus promoting the progress of colorectal cancer. Wang et al. (2020) discovered that METTL3 enhanced the stability of HK2 through YTHDF1-mediated m6A modification, thus promoting the Warburg effect of cervical cancer, and finally promoting the malignant proliferation of tumors. On the project, we attempted to analyze whether the expression level of HK2 is associated to the m6A modification in ESCA. We discovered that the expression of HK2 was remarkable associated with IGF2BP2, HNRNPA2B1, YTHDC2, and YTHDF2. We also discovered that the expression levels of YTHDF3 and YTHDF2 were remarkable increased in the group with high expression of HK2. We believe that HK2 gene may have a regulatory relationship with m6A modification, and this regulatory mechanism is involved in glycolysis and malignant proliferation of ESCA.

CONCLUSION
In conclusion, we verified that HK2 is overexpressed in ESCA, and its expression level is associated to clinical case characteristics. The expression level of HK2 is strongly associated to the degree of immune cell infiltration. HK2 is associated with m6A modification, which may enhance the stability of HK2 through m6A modification, and thus promote the glycolytic effect and malignant proliferation of ESCA. HK2 can be used as a potential biological target for diagnosis and therapy of ESCA.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Taihe Hospital Affiliated of Hubei University of Medicine. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
X-SL conceived the project and wrote the manuscript. X-SL, Y-JC, F-YL, and R-MW participated in data analysis.