Upregulated Glycolysis Correlates With Clinical and Transcriptomic Signicances in Head and Neck Squamous Cell Carcinoma

Altered metabolism is an emerging hallmark of cancer. Cancer cells preferentially utilize glycolysis for energy production, termed “aerobic glycolysis.” In this study, we performed a comprehensive analysis of the glycolysis status in the tumor microenvironment (TME) of head and neck squamous cell carcinoma (HNSCC) using data from The Cancer Genome Atlas database. We rst divided 520 patients with HNSCC into two groups based on the mRNA expression of 16 glycolysis-related genes. The glycolysis-high signature positively correlated with human papillomavirus-negative tumor type, advanced T factor, and unfavorable prognosis. The gene set enrichment analysis revealed upregulation of several pathways, including interferon-alpha response, myc targets, hypoxia, epithelial-mesenchymal transition, transforming growth factor-β signaling, and interleukin 6-Janus kinase-signal transducer and activator of transcription 3 signaling, in the glycolysis-high group. Immune cell enrichment analysis revealed decreased inltration of T cells, dendritic cells, and B cells in the glycolysis-high group, suggesting impaired tumor antigen presentation, T cell activation, and antibody production in TME. Moreover, the expression of TGFB1, CD274, and PDCD1LG2, which facilitate immunosuppression in the TME, was upregulated in the glycolysis-high group. Collectively, these ndings suggest the potential of glycolysis monitoring as a biomarker for tumor progression and immunosuppression in the TME of HNSCC. revealed the of various representing in the glycolysis-high group. Immune cell enrichment analyses showed decreased inltration of T cells, DCs, and B cells and increased inltration of monocyte and neutrophil in the glycolysis-high group. Moreover, the expression of immunosuppressive genes, including TGFB1, CD274, and PDCD1LG2, was upregulated in the glycolysis-high group. These results suggest the potential of monitoring glycolysis as a biomarker for tumor progression and immunological status in the TME of


Introduction
Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide, with a reported incidence of 890,000 new cases and 450,000 deaths annually 1,2 . The major risk factors of HNSCC include exposure to tobacco-derived carcinogens and excessive alcohol consumption 3 . Moreover, prior infection with oncogenic strains of human papillomavirus (HPV) has recently been recognized as a risk factor of oropharyngeal cancers, especially in the USA and Western Europe 4 . As the 5-year survival rate of HNSCC is still 66%, new therapeutic modalities have been explored. Recent breakthroughs in cancer therapy have resulted in novel immunotherapy regimens targeting the immune checkpoint programmed death 1 (PD-1) in patients with HNSCC 5,6 . However, the survival bene t is limited to only 20%-30% patients, and new biomarkers that predict the e cacy of anti-PD-1 therapy have been investigated. Accumulating evidence has indicated that in HNSCC, tumor cells orchestrate a highly immunosuppressive tumor microenvironment (TME) via the production of immunosuppressive mediators, recruitment of various stromal cells, expression of immune checkpoint ligands, and downregulation of human leukocyte antigen expression 7,8 , which is referred to as immune evasion 9 . As immune checkpoint inhibitors, including anti-PD-1 antibodies, target the interaction between tumor cells and T cells, a comprehensive understanding of complex immune networks in the TME is needed to establish a new biomarker for immunotherapies.
Reprogramming of energy metabolism is an emerging hallmark of cancers 10,11 . Even during aerobic conditions, cancer cells preferentially utilize glycolysis for producing energy, termed "aerobic glycolysis" 12,13 . Increased glycolysis correlates with aggressive tumor progression, treatment resistance, and unfavorable prognosis in various cancers, including HNSCC [14][15][16][17] . The upregulation of glucose transporter-1 (GLUT1), which transports glucose into the cytoplasm, increases glucose utilization and is associated with aggressive behavior in several cancers 18 . Increased glucose uptake via GLUT1 has been investigated in clinics using positron emission tomography, which monitors the uptake of a radiolabeled glucose analog, 18 F-uorodeoxyglucose. Moreover, both oncogenic signaling pathways and hypoxia induce aerobic glycolysis in cancer cells. Activation of oncogenes such as RAS, MYC, and MAPK and PI3K-AKT signaling as well as TP53 mutation, a tumor suppressor gene, are associated with increased glycolysis [19][20][21][22][23] . The hypoxic TME also alters metabolic reactions in tumor cells by upregulating the expression of glucose transporters and enzymes of the glycolytic pathway, including transcription factors, hypoxia-inducible factor (HIF) 1α and HIF2α [24][25][26] . However, the primary mechanism that regulates glycolysis in cancer cells still remains unclear. Furthermore, the relationship between upregulated glycolysis and clinicopathological features is elusive in various cancer types. Notably, a few studies have revealed the interplay between tumor glycolysis and immune evasion 15,27,28 . A comprehensive analysis that focuses on the relationship between glycolysis and immunological signi cance in the TME is warranted.
In the present study, we performed a comprehensive analysis of glycolysis in HNSCC by analyzing transcriptome and clinical data regarding HNSCC obtained from The Cancer Genome Atlas (TCGA) database. Based on the mRNA expression of 16 glycolysis-related genes, we segregated patients with HNSCC into two groups and performed a comprehensive analysis to compare the clinical and biological signi cance between the glycolysis-high and glycolysis-low groups.

Results
Glycolysis-high gene signature positively correlated with locally advanced tumors and unfavorable prognosis We rst performed non-supervised hierarchical clustering of 520 patients with HNSCC based on data on the expression of 16 glycolysis-related genes obtained from TCGA database ( Figure 1a). Then, patients were divided into two groups; both groups were compared in terms of clinical parameters, including HPV status, primary lesion, T factor, N factor, M factor, and tumor-node-metastasis (TNM) stage (Table 1). In the glycolysis-low group, 42 (35.2%) patients were HPV-positive, whereas in the glycolysis-high group, 55 (13.7%) patients were HPV-positive. Regarding primary lesions, the proportion of oral cavity lesions was higher in the glycolysis-high group (71.1%) than in the glycolysis-low group (55.4%); however, the proportion of oropharynx lesions was lower in the glycolysis-high group (5.5%) than in the glycolysis-low group (17.6%). In the glycolysis-low group, the proportion of patients with Tis, T1, and T2 tumors (51.3%) was higher than that of patients with T3 and T4 tumors (48.7%). However, in the glycolysis-high group, the proportion of patients with Tis, T1, and T2 tumors was lower (37.2%) than that of patients with T3 and T4 tumors (62.8%). No differences in the N factor, M factor, and TNM stage were observed between the groups. Additionally, we performed survival analyses to compare survival between the groups ( Figure  1b). The glycolysis-high group exhibited shorter overall survival (OS) than the glycolysis-low group, especially in patients who were HPV-positive. Additionally, the glycolysis-high group exhibited shorter disease-free survival than the glycolysis-low group.

Glycolysis-high gene signature correlated with various hallmark pathways
To investigate the biological differences between the glycolysis-low and glycolysis-high groups, we compared the groups and identi ed 2366 differentially expressed genes (DEGs), including 459 upregulated and 1907 downregulated genes in the glycolysis-high group (Figure 2a, Suppl. Table 1). To explore the biological functions of DEGs and related pathways, we conducted Gene set enrichment analysis (GSEA; Figure 2b, 2c). In the glycolysis-high group, 10 hallmark pathways were upregulated and two pathways were downregulated (FDR< 0.10). Various pathways representing cancer hallmarks, such as interferon-alpha (IFN-α) response, myc targets, hypoxia, epithelial-mesenchymal transition (EMT), transforming growth factor (TGF)-β signaling, and interleukin 6 (IL6)-Janus kinase (JAK)-signal transducer and activator of transcription 3 (STAT3) signaling, were upregulated in the glycolysis-high group.

Correlation between glycolysis activity and immune cell in ltration
We next calculated the enrichment scores of 16 basic immune cell types to assess the correlation between glycolysis activity and immune cell in ltration. The immune score and the enrichment scores of CD4+ T cells, CD8+ T cells (not signi cant), B cells, and dendritic cells (DCs) were lower in the glycolysis-high group than in the glycolysis-low group. However, the enrichment scores of Th2 cells, NKT cells, monocytes, and neutrophils were higher than those in the glycolysis-low group than in the glycolysis-high group (Figure 3a). We also examined the correlation between each glycolysis gene and immune cell enrichment scores, which indicated the same trends ( Figure 3b).

Glycolysis activity correlated with the immunosuppressive microenvironment
We compared the expression of immune-related genes to emphasize immunological differences between glycolysis groups (Figure 4a). The expression of immunosuppressive genes, including TGFB1, CD274, and PDCD1LG2, was higher in the glycolysis-high group than in the glycolysis-low group (Figure 4b). In contrast, the expression of immune checkpoint molecules, including CTLA4, PDCD1, and TIGIT, was lower in the glycolysis-high group than in the glycolysis-low group (Figure 4b). No difference in the expression of the immunostimulatory genes IFNG and GZMB was observed between glycolysis groups (Suppl. Figure  1). Similar trends were observed between each glycolysis gene and immune-related genes (Figure 4c).

Discussion
Altered metabolism is an emerging hallmark of cancer 10,11 . Various cancers exhibit upregulated glycolysis and downregulated mitochondrial oxidative phosphorylation. Immune evasion has recently attracted attention due to recent advances in cancer immunotherapy. In the present study, we performed a comprehensive analysis of the transcriptome and clinical data regarding HNSCC obtained from TCGA database. The glycolysis-high signature positively correlated with HPV-negative tumor type, advanced T factor, and unfavorable prognosis. GSEA revealed the upregulation of various pathways representing cancer hallmarks in the glycolysis-high group. Immune cell enrichment analyses showed decreased in ltration of T cells, DCs, and B cells and increased in ltration of monocyte and neutrophil in the glycolysis-high group. Moreover, the expression of immunosuppressive genes, including TGFB1, CD274, and PDCD1LG2, was upregulated in the glycolysis-high group. These results suggest the potential of monitoring glycolysis as a biomarker for tumor progression and immunological status in the TME of HNSCC.
Non-supervised hierarchical clustering analysis revealed a positive correlation between glycolysis and advanced T factors. In normal cells, mitochondrial oxidative phosphorylation is widely used to e ciently generate energy. However, cancer cells often rely on aerobic glycolysis, resulting in remarkably upregulated glycolysis, although not an e cient strategy to produce adenosine 5'-triphosphate (ATP). Accelerated cell proliferation has been implicated in the dependence of cancer cells on glycolysis. For rapid proliferation, cells require nutrient metabolism to generate nucleotides, amino acids, and lipids rather than e cient ATP production. In cancer cells, upregulated glycolysis results from altered metabolic pathways that provide biomass for cell proliferation 12,29 . Consequently, the advanced T factor in the glycolysis-high group could be attributed to accelerated cell proliferation, which relies on aerobic glycolysis. Furthermore, we observed a positive correlation between the HPV-negative tumor type and the glycolysis-high signature. Comprehensive analysis of the TCGA database indicated frequent alterations in TP53 and CDKN2A in HPV-negative HNSCC, but not in HPV-positive HNSCC. TP53 is a negative regulator of glycolysis 30 . Therefore, upregulated glycolysis may be a biological feature of HPV-negative carcinogenesis. In addition, upregulated glycolysis is an unfavorable prognostic factor. Upregulation of glycolysis in cancer cells results in the accumulation of lactate in the TME 12 . A high lactic acid concentration in the TME inhibits the cytotoxic activity of cytotoxic T cells, resulting in immune evasion of tumors 31 . A lactate-enriched TME also mediates increased invasion, apoptosis resistance, and increased survival and proliferation of tumor cells [32][33][34] . Therefore, advanced T factor and a lactateenriched TME may be associated with shorter survival in glycolysis-upregulated HNSCCs.
We further compared the two glycolysis groups using GSEA. Various pathways that represent cancer hallmarks were upregulated in the glycolysis-high group. The IFN-α response pathway exhibited the highest normalized enrichment score. IFN-α is a cytokine family that belongs to type I IFNs 35 . Based on its stimulatory effects on innate and adaptive immunity, IFN-α has been widely used in clinical oncology to treat various cancer types, especially hematological malignancies 36 . However, recent studies have demonstrated the pro-tumoral role of IFN-α signaling. In HNSCC, IFN-α signaling promotes the immunosuppressive TME by activating the PD-1-programmed death ligand-1 (PD-L1) axis 37 . Therefore, the upregulation of PD-1 ligand genes, CD274 and PDCD1LG2, in the glycolysis-high group could be related to the upregulated IFN-α response. Although the abundance of IFN-γ in the TME has been recognized as a key mechanism that induces PD-L1 expression in tumor cells 38 , our results suggest that IFN-α could also be an important mechanism for PD-L1 induction in glycolysis-upregulated HNSCC. Moreover, growing evidence suggests that IFN-α signaling contributes to increased migration, drug resistance, and immune evasion in breast cancer 39 . Further investigations focused on the relationship between IFN-α signaling and glycolysis in the TME of HNSCC are needed. Furthermore, GSEA revealed that that the EMT pathway and TGF-β signaling pathway were upregulated in the glycolysis-high signature. EMT is a process in which cells lose their epithelial features and acquire mesenchymal traits, including motility, invasiveness, and various characteristics of stem cells 40 . TGF-β is a well-known cytokine family and a major driver of EMT 41 . Accumulating evidence suggests that TGF-β signaling also enhances glycolysis by activating various glycolytic enzymes, including GLUT1, hexokinase 2, 6phosphofructo-2-kinase/fructose-2,6-biphosphatase 3, 6-phosphofructo-1-kinase, pyruvate kinase M2, and lactate dehydrogenase type A 42 . Consistent with these ndings, TGF-β signaling may play a vital role in the induction of both glycolysis and EMT in HNSCC. Additionally, the signi cant upregulation of myc target pathways and hypoxia signaling was consistent with their roles in the induction of glycolysis 20,24-26 .

Immune cell enrichment analysis revealed a negative correlation between glycolysis and in ltration of CD4+ T cells, CD8+ T cells, B cells, and DCs. The interaction between DCs and T cells is an important axis
of the anti-tumor immune response through tumor-antigen presentation and T cell activation 43 . The negative effects of the lactate-enriched TME on T cell proliferation, cytotoxic activity of T cells, and maturation of DCs support our observation 44,45 . Moreover, the comparison of immune-related gene expression between the glycolysis groups revealed higher expressions of CD274, PDCD1LG2, and TGFB1 in the glycolysis-high group, which facilitated T cell exhaustion and dysfunction in the TME 46 . Therefore, the decreased in ltration of T cells and DCs in the glycolysis-high group may be related to increased immune evasion in the glycolysis-elevated TME. The lower gene expression of immune checkpoint molecules on T cells, such as CTLA4, PDCD1, and TIGIT, may re ect lower T cell in ltration in the glycolysis-high group.
Recent advances in cancer immunotherapy targeting solid tumors have mainly focused on T cells; however, the anti-tumor roles of B cells through antigen presentation and antibody production have been recently recognized 47 . The signi cant decrease in B cell in ltration in the glycolysis-high group may also re ect immune evasion in the TME of HNSCC. Moreover, both monocytes and neutrophils represent diverse functions and subsets, facilitating pro-tumoral and anti-tumoral immune responses 48,49 . Although increased in ltration of both monocytes and neutrophils in the glycolysis-high group could be an aspect of immune evasion in the glycolysis-elevated TME, further characterization of these immune cell types would be necessary.
There are limitations to the present study. We performed a comprehensive analysis using bulk RNA sequencing data. In the case of bulk RNA sequencing, gene expression values represent not only mRNA expression in tumor cells but also that in all other cell types existing in the TME. This method should work well for genes that are speci c to certain cell types; however, for genes that are commonly expressed in various cell types, single-cell RNA sequencing would be better suited for a more precise analysis of the TME. Accumulation of public archives of single-cell RNA sequencing data is needed for an in-depth characterization of the metabolic state in the TME.
In conclusion, we elucidated the correlation between upregulated glycolysis and several clinical and transcriptomic signi cances in HNSCC by analyzing data from TCGA database. Glycolysis monitoring has the potential as a biomarker for tumor progression and immunological status in the TME of HNSCC. Hierarchical clustering of patients with glycolysis-related genes All patients with HNSCC underwent non-supervised hierarchical clustering based on z-scores of the expression of 16 glycolysis-related genes. The following 16 glycolysis-related genes were used based on the glycolysis core pathways of the Reactome pathway database: PFKP, PGK1, PGAM1, PGAM4, ENO1, PKM2, LDHA, TPI1, GAPDH, GPI, ALDOA, HK1, SLC2A1, HK2, ALDOC, and PKLR 50 . Patients were divided into two groups-the glycolysis-low and glycolysis-high groups-using the cutree R function based on clustering results. The heat map was illustrated using the pheatmap R package.

Differentially expressed gene analysis
In the glycolysis-low and glycolysis-high groups, DEGs were identi ed using the ExperimentHub R package and DESeq2 R package. DEGs were ltered using the threshold |log 2 FC|≥1 and adjusted P-value of<0.05. A volcano plot was used to visualize DEGs using the calibration R package.
Gene set enrichment analysis GSEA (GSEA v4, Broad Institute) was performed to identify pathways upregulated in the glycolysis-high group. The normalized enrichment score, P-value, and false discovery rate (FDR) q-value were calculated for each gene set using the Hallmark pathway database. A correlation plot was constructed using the corrplot and pheatmap R packages.

Cell type enrichment analysis
The xCell tool, a gene signature-based method, was employed to evaluate the enrichment of various cell types in 520 HNSCC tissues. Among the 64 cell types, the enrichment scores of 13 basic immune cell types were calculated using the xCell R package.

Statistical analysis
Data were analyzed using R (version 4.0.3; The R Foundation for Statistical Computing, Vienna, Austria) in combination with R studio version 1.3.1093 (R studio, Boston, MA, USA) and GraphPad Prism version 8 (GraphPad Software, San Diego, CA, USA). Student's t-test was used to compare continuous variables between the groups. The chi-square test for independence and Fisher's exact test were used to compare categorical variables. Pearson's correlation coe cient was utilized to evaluate the correlation between two continuous variables. Two-sided P-values of <0.05 were considered statistically signi cant. Survival curves were calculated using the Kaplan-Meier method and compared using the log-rank test.   Differentially expressed genes and upregulated pathways in the glycolysis-high group a, Volcano plot of differentially expressed genes in the glycolysis-high group. Red dots represent upregulated genes (Padj < 0.05, log2FC > 1), whereas blue dots represent downregulated genes (Padj < 0.05, log2FC < -1). b, Upregulated and downregulated hallmark pathways in the glycolysis-high group obtained by GSEA (FDR < 0.10). c, GSEA plots of pathways shown in b. Padj, adjusted P-value; GSEA, gene set enrichment analysis.

Figure 3
Correlation between glycolysis activity and immune cell enrichment scores a, Bar graphs of immune cell enrichment scores in the glycolysis-high and glycolysis-low groups. b, Correlation matrix displaying Rvalues for assessing the correlation between the expression of 16 glycolysis-related genes and immune cell enrichment scores across 520 HNSCC tissues. The glycolysis-related genes were clustered by hierarchical clustering based on R values. HNSCC, head neck squamous cell carcinoma.

Figure 4
Correlation between glycolysis activity immune-related genes a, Heat map of immune-related gene expression in the glycolysis-high and glycolysis-low groups. b, Violin plots of normalized gene expressions that were signi cantly upregulated or downregulated in the glycolysis-high group shown in a. c, Correlation matrix displaying R-values for assessing the correlation between the expression of 16 glycolysis-related genes and that of immune-related genes across 520 HNSCC tissues. Glycolysis-related genes and immune-related genes were clustered by hierarchical clustering based on the R values. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. HNSCC, head neck squamous cell carcinoma.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. TakahashietalSuppl.Table1.xlsx GlycolysisSupple.Figure1.pdf