A novel defined hypoxia-related gene signature to predict the prognosis of oral squamous cell carcinoma

Introduction

As the most common oral cancer, oral squamous cell carcinoma (OSCC) has become a global health challenge (1). Risk factors for OSCC include smoking, drinking, viral infections (Epstein-Barr virus, human papillomavirus and herpes simplex virus), chewing betel nuts, occupational exposure to carcinogens, immune deficiency, radiation, diet, and genetic susceptibility (2).

Exploring the tumor microenvironment (TME) has attracted more attention to studying the mechanism of solid tumors (3). In TME, abnormal tumor growth consumes more oxygen and insufficient blood supply, making hypoxia a common phenomenon and an important sign in solid tumors (4). Many studies have confirmed that hypoxia is closely related to treating tumors (5,6). Hypoxia can induce changes in gene expression and might be an independent predictor of tumors (7). Furthermore, hypoxia has been shown to regulate different factors in TME. Specifically, hypoxia can induce immune cells and immune cytokines, impair the immune response against cancer (8), and further affect immunotherapy (9). In head and neck cancer (HNC), the correlation between hypoxia and the efficacy of immune checkpoint inhibitors has been confirmed, and subsets of HNC have been identified based on the characteristics of hypoxia immunity (10-12).

In OSCC, hypoxia can induce epithelial-mesenchymal transition (EMT) and promote tumor progression (13). Similarly, increasing the level of hypoxia can affect the sensitivity of OSCC to chemotherapeutic drugs and reduce the efficacy (14). In general, hypoxia has not been fully explored in OSCC’s TME. In-depth exploration of the mechanism of hypoxia and OSCC, and the selection of hypoxia-related prognostic targets, are of great significance for understanding tumors and improving prognosis.

In this study, we first downloaded genes related to the hypoxia index from the Molecular Signatures Database (MSigDB). According to the single-sample Gene Set Enrichment Analysis (ssGSEA), different data from The Cancer Genome Atlas (TCGA)-OSCC and GEO database, the enrichment score (ES) is calculated for the sets samples (GSE25099, GSE30784, and GSE41613) and defined as hypoxia potential index (HPI). HPI was used to represent the hypoxic state in tumor samples. Next, we compared the difference in HPI between OSCC tumor samples and normal samples, and the relationship between HPI and patient survival. We also further explored the relationship between the impact of different HPI samples on immune cell infiltration in TME and the relationship between gene mutations. We found that the increase in HPI is related to the higher TME infiltration of eosinophils, macrophages, neutrophils, helper cells of T and Th1 cells, and the mutation of NSD1 may indicate the high HPI of OSCC. Finally, we constructed a prognostic model of 6 sets of hypoxia-related genes (PGK1, JMJD6, S100A4, SLC2A3, DDIT4 and HK1) in OSCC. These results help us better understand the role between immune cells and hypoxia and discover new therapeutic targets related to hypoxia in OSCC. In this study, we used a more rigorous prediction model construction approach to screen prognostic genes by cox risk regression, random survival forest variable hunting, and for genes with prognostic value, we used a permutation approach to screen the optimal combination for predicting patient prognosis. In addition, we also investigated the relationship between hypoxia and mutations, tumor-related pathway. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/atm-21-4990).

Methods

Data download and preprocessing

TCGA is a cancer genomics database (https://cancergenome.nih.gov/), which includes genome information from many cancer samples and corresponding normal samples (15,16). Transcriptome sequencing data (FPKM and count expression profile) from patients with TCGA-HNSC, corresponding clinical information, and mutation data, which were downloaded through the University of California Santa Cruz Xena (UCSC Xena; https://xena.ucsc.edu/) platform (17). Samples with hypopharynx, larynx, oropharynx, tonsil tissue types were excluded, and patient data with incomplete follow-up information were excluded. Finally, 326 OSCC samples were included in our study. The GEOquery R package was used to download OSCC sequencing data (GSE25099: 57 tumor samples and 22 normal samples; GSE30784: 17 dysplasia samples, 45 normal samples and 167 tumor samples; GSE41613: 76 tumor samples) from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). Among them, GSE41613 has complete patient follow-up information. The study design and characteristics of the samples used at each analysis stage are shown in Figure 1. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Figure 1 The general design of the research and the flow of the study.

The HPI of hypoxia index

Next, we continue to explore the relationship between HPI and common tumor-related pathways. The NFkB signaling pathway (18), the transforming growth factor (TGF)-beta signaling pathway (19-21), epithelial-to-mesenchymal transition (22,23), glycolysis (24,25), and angiogenesis (26,27) as common tumor-related pathways play a vital role in various cancers. The sets of hallmark genes for hypoxia (200 genes) and other tumor-related pathways [TNFA signaling via NFKB (200 genes), TGF beta signaling (54 genes), Epithelial Mesenchymal Transition (200 genes), Glycolysis (200 genes), and Angiogenesis (36 genes)] were downloaded from the Molecular Signatures Database (MSigDB). The enrichment score (ES) of these gene sets was calculated by ssGSEA in the ‘GSVA’ package R (28,29), and the hypoxia gene list enrichment score was defined as the hypoxia potential index (HPI) to computationally dissect hypoxia levels/trends of tissue samples. Also, the HPI of the GEO datasets was evacuated by ssGSEA.

Survival analysis and intergroup difference of HPI

Calculate the HPI in each sample of TCGA-OSCC and GSE41613, and sort the samples. ‘Survival’ and ‘survminer’ are used to determine the best cut-off value to group samples and draw survival curves for high HPI groups and low HPI groups. Welch’s t-test was used to compare the intergroup difference between tumor and normal samples, and the ‘ggstatsplot’ package was used for plotting.

Tumor sample clustering based on HPI

To further explore the relationship between HPI and OSCC subtypes, we performed a gene expression differential analysis to identify differentially expressed genes between the HPI_high group and the HPI_low group. Count expression profile and edgeR R package were used for differential analysis of gene expression. A total of 141 differential expression genes (|logFC|>1.2, P<0.01) were selected for a consistent cluster analysis. Different clustering variables (k) determine the cluster with the highest correlation within the group and the lowest correlation between the clusters. The ConsensusClusterPlus package and TCGA-OSCC data were used for the identification of subgroups (30). The ‘survival’ package and the ‘survminer’ package were used for survival analysis between subgroups.

HPI-related Gene Set Enrichment Analysis (GSEA) and Gene Ontology (GO) function enrichment

GSEA helps us determine the correlation between HPI and multiple pathways in OSCC (29). In the TCGA-OSCC tumor samples, we sorted the samples according to HPI and identified differentially expressed genes between the high-HPI and low-HPI samples. Functional enrichment analysis of GSEA and GO based on HPI-related differential genes was completed with the ClusterProfiler package (31), and the GOplot package is used to visualize the visualization of biological data.

The relationship between HPI and tumor-related pathways and immune cell infiltration

The study by Bindea et al. provided marker genes for 24 immune cells (32); based on the marker genes, the infiltration level of 24 immune cells in samples from this study was calculated by ssGSEA (29). Wilcoxon rank sum tests were used to compare the difference in immune cell infiltration between samples from the high HPI group and the low HPI group. Similarly, the enrichment scores of tumor-related pathways (TNFA signaling via NFKB, TGF BETA signaling, Epithelial Mesenchymal Transition, Glycolysis, and Angiogenesis) and immune cell infiltration levels were also compared.

Correlation between hypoxia and tumor-related pathways, hypoxia, and immune cells

Whether hypoxia in OSCC is related to common tumor-related pathways is worthy of further exploration. The Pearson correlation coefficient was used in the TCGA-OSCC and GSE41613 datasets to identify the correlation between hypoxia and tumor-related pathways. Furthermore, we also analyzed the correlation between different immune cells and hypoxia in OSCC.

The relationship between HPI and gene mutation

From UCSC Xena (xenabrowser.net), we downloaded the somatic mutation data of TCGA-HNSC. Samples of somatic mutation data in the samples we selected before were used for the next analysis. The independence test was used to calculate the correlation between gene mutation and the hypoxia index.

Construction of hypoxia gene prognosis model

There were 200 genes related to hypoxia. First, we used univariate Cox hazard analysis to identify hypoxic genes with prognostic values, genes with P value <0.05 were selected as prognostic genes. We further screened the prognostic genes using the random survival forest variable hunting (RSFVH) algorithm, after passing the gene combinations, the combination with significant P value and less number of genes was taken using log-rank test and defined as prognostic signature. By “rms” R package, the nomogram was constructed based on prognostic genes which selected by last step, the nomogram was used to calculate each gene’s risk score and predict 1- and 3-year survival of OSCC patients in the TCGA dataset. The calibration plots were used to evaluate the performance of the prognostic nomogram in the TCGA-OSCC and GSE41613 datasets.

Construction of a protein-protein interaction (PPI) network of hypoxia-related genes

The String database (https://string-db.org/) is a searchable database of known PPIs and predicted PPIs. After univariate Cox hazard analysis, genes with P<0.05 were selected for the construction of PPI.

Statistical analysis

All analyses were performed by R version 4.0.5. Welch’s t test, Wilcoxon rank-sum test were used to compare differences between groups and Kaplan-Meier method was used to compare the prognostic value of different subgroups. A P value <0.05 was considered statistically significant.

Results

Survival prediction of HPI in OSCC

According to the best HPI cutoff value, tumor samples in TCGA-OSCC and GSE41613 were divided into the HPI high group and the HPI low group. Compare the prognosis between the two groups and draw the patient’s survival curve (Figure 2A,2B). The results showed that in samples of TCGA-OSCC and GSE41613, low HPI was associated with a better survival time (GSE41613: P=0.00022; TCGA-OSCC: P=0.013). In the GSE25099 data set, we compared the difference in HPI between tumor tissue and normal tissue (Figure 2C). The results suggest that HPI in tumor tissue is higher than in normal tissue (P=0.030). In the GSE30784 data set, we compared HPI differences between normal and tumor (P=4.16e-10), normal and dysplasia (P=0.007), dysplasia and tumor (P=0.019), respectively (Figure 2D). HPI in normal dysplasia and tumor tissues gradually increased, which was consistent with the survival curve.

Figure 2 The survival relationship and the analysis of the difference of HPI in OSCC. (A,B) Survival curves in patients with GSE41613 and TCGA-OSCC tumors; (C,D) in GSE25099 and GSE30784, the difference in HPI between tumor patients and normal patients was shown as violin charts. HPI, hypoxia index.

Identification of hypoxia-related tumor clusters

Next, we downloaded the OSCC sample count expression matrix from UCSC and grouped the samples according to HPI. Then, using the edgeR package, the difference between the OSCC high HPI group and the low HPI group was analyzed, and 141 hypoxia-related differential genes (|logFC| >1.2, P<0.01) were obtained. To further explore the association between differential genes and OSCC subtypes, we performed a consistent cluster analysis in 326 tumor samples in the TCGA cohort. By increasing the clustering variable (k) from 2 to 10, we found that when k=2, the correlation within the group is the highest, while the correlation between the groups is low (Figure 3A). The overall survival time (OS) was compared between the two groups, but no obvious differences were found (P=0.088, Figure 3B). Subsequently, a heat map was drawn to show the relationship between HPI and patient clinical information (tumor grade stage, T, N and M) between the two clusters (Figure 3C). The results suggest that there is no significant cluster between HPI and the patient’s clinical characteristics.

Figure 3 Clustering of tumors based on differential genes of hypoxia. (A) Consistent grouping of 330 tumor patients (k=2); (B) survival curve of two clusters (P=0.088); (C) heat maps of the clinical characteristics (tumor grade, T, N, and M) of patients in the two sets of high HPI and low HPI samples. HPI, hypoxia index.

Enrichment analysis of differential genes

We first performed GSEA on the hypoxia-related gene set in TCGA. We found that 141 DEGs are enriched in Base excision repair, DNA replication, Fanconi anemia pathway, Mismatch repair, Ribosome and Ubiquinone, and other terpenoid quinone biosynthesis pathways (Figure 4A). The results of the GSEA results suggest that OSCC may be involved in the damage repair process in tumor tissues, which is similar to the role of hypoxia in other tumors (33-35). Through the GO annotation, analyze the potential functions of differential genes. Muscle system process, extracellular matrix structural constituent, and collagen binding were the most common GO terms for biological processes, cellular components, and molecular functions, respectively (Figure 4B).

Figure 4 Enrichment analysis of differential genes. (A) GSEA of hypoxia-related genes; (B) GO function enrichment analysis of hypoxia-related differential genes. GSEA, Gene Set Enrichment Analysis; GO, Gene Ontology.

The relationship between HPI and immune cells

In addition, we compared the relationship between immune cell infiltration and HPI in the TCGA-OSCC and GSE41613 datasets. The results showed that in TCGA-OSCC and GSE41613, the common differentially infiltrated immune cells include eosinophils, macrophages, neutrophils, helper T cells, and Th1 cells (Figure 5A,5B). Infiltration of these immune cells was significantly increased in the high HPI group. However, the five tumor-related pathway enrichment scores in the high HPI group and the low HPI group were also compared using the Wilcoxon rank sum test. They showed significant differences in the different HPI groups; the HPI high group showed significantly higher enrichment scores (Figure 5C,5D).

Figure 5 Relationship between HPI and immune cell tumor infiltration. (A,B) In GSE41613 and TCGA-OSCC, the difference in immune cell infiltration between the high HPI group and the low HPI group; (C,D) In GSE41613 and TCGA-OSCC, the relationship between tumor-related pathways and immune cell infiltration. *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001. HPI, hypoxia index; TCGA, The Cancer Genome Atlas; OSCC, oral squamous cell carcinoma.

Correlation between hypoxia and tumor-related pathways, hypoxia, and immune cells

Given that both immune cell tumor infiltration and tumor-related pathways are significantly associated with hypoxia, we continue to study the correlation between hypoxia and these pathways. The results showed that hypoxia and tumor-related pathway ES (TNFA SIGNALING VIA NFKB, TGF BETA SIGNALING, EPITHELIAL MESENCHYMAL TRANSITION, GLYCOLYSIS, and ANGIOGENESIS) in TCGA-OSCC and GSE41613 have significant correlations (Figure 6A). The correlation between HPI and immune cell tumor invasion showed a correlation between HPI and eosinophils, macrophages, neutrophils, and Th1 cells (Figure 6B).

Figure 6 The correlation between hypoxia, tumor-related pathways, and immune cells. (A) Correlation between hypoxia and tumor-related pathways; (B) the correlation between hypoxia and immune cell TME infiltration. TME, tumor microenvironment.

Relationship between HPI and gene mutation

Gene mutations are frequent in tumors and can affect the treatment of tumors (36,37). The relationship between hypoxia and gene mutations in OSCC tumors is still unclear. Here, we perform independent tests on HPI and mutant genes to screen for highly related mutant genes to HPI. The results showed that NSD1 was correlated with the HPI level in OSCC (P=5.55e−04, Figure 7).

Figure 7 The relationship between hypoxia and gene mutation. In patients with OSCC, mutations in the NSD1 gene are significantly associated with hypoxia. OSCC, oral squamous cell carcinoma.

Construction of a prognostic model of hypoxia-related genes

We performed a univariate Cox hazard analysis on the 200 hypoxia genes and found that 36 genes were associated with the prognosis of OSCC patients (P<0.05, Figure 8A). Then we selected 10 genes (PGK1, HS3ST1, ALDOC, JMJD6, ADM, S100A4, SLC2A3, DDIT4, TIPARP, and HK1) using the random forest-supervised classification algorithm (Figure 8B). The prognostic risk model that these genes constitute - The log10 plog-rank value is sorted, and the results show that the prognostic model is composed of 6 genes (PGK1, JMJD6, S100A4, SLC2A3, DDIT4 and HK1) has the highest ranking (Figure 8C). A total of 326 patients from TCGA-OSCC were used to establish a predictive nomogram for the six genes to predict overall survival at 1 and 3 years based on the Cox regression model (Figure 8D). The calibration curves estimating survival probability at 1 and 3 years showed good agreement between the nomogram-predicted and observed values (Figure 8E,8F).

Figure 8 Construction of a prognostic model of genes related to hypoxia. (A) The volcano plot shows the P value of the hypoxia index gene in the univariate Cox regression analysis. (B) Random survival forest analysis screened 10 genes. (C) After Kaplan-Meier analysis of 2¹⁰−1=1,023 combinations, we rank the survival P value of the first 20 gene sets, and the gene set composed of 6 genes has the best survival assessment. (D,E) nomogram to show the relationship between the six genes and the clinical characteristics of OSCC patients. (F) Draw a correction curve to evaluate the accuracy of the risk score prediction. OSCC, oral squamous cell carcinoma.

Construction of the PPI network of hypoxia genes

Finally, according to the result of the univariate Cox hazard analysis on hypoxia genes, 36 genes with P<0.05 were used to construct a PPI network (Figure 9). We found that the six hub genes (PGK1, JMJD6, S100A4, SLC2A3, DDIT4, and HK1) were used to construct a prognostic model. PGK1, SLC2A3 and HK1 have high interactions with other genes in the hypoxia gene set. This result further indicates that the hub gene we have screened may play a crucial role in OSCC.

Figure 9 Construction of the PPI network of hypoxia genes. A total of 36 genes with P<0.05 were used to construct a PPI network. PPI, protein-protein interaction.

Discussion

As the most common oral cancer, OSCC has a higher recurrence rate and risk of metastasis, lowering the overall survival rate of patients and becoming an important risk that endangers the life safety of patients (38,39). Hypoxia has been found in previous studies to play an important role in the prognosis and treatment of breast cancer (40,41), nonsmall cell lung cancer (42), and ovarian cancer (43). In this study, we have comprehensively explored the various mechanisms of action of hypoxia indicators in OSCC.

In the hypoxic microenvironment, the ability to invade and metastasize tumors is strengthened (6). In oral cancer, changes in the expression of genes related to hypoxia have been described (44). Higher HPI is associated with a poor prognosis in patients with OSCC (45) and mediated apoptosis through two pathways (both intrinsic (mitochondrial) and extrinsic (cell death receptor-mediated) pathways) (46). We first compare the HPI in the OSCC samples and verify the higher HPI in tumor tissues through different datasets. Survival analysis showed that higher HPI was associated with a worse patient prognosis in the TCGA and GEO data sets. In this study, in addition to conventional cox risk regression, we also used random survival forest variable hunting, and for genes with prognostic value, we used permutation to filter the optimal combination for predicting patient prognosis. In addition, we also explored the relationship between hypoxia and immune cells, tumor-related pathways, and explored the mechanism of the role of hypoxia in tumor prognosis from various aspects. Furthermore, we elaborated on the impact of HPI on OSCC from three aspects: immune cells, gene mutations, and tumor prognosis.

Hypoxia in tumor tissues has been proven to reduce the effect of immunotherapy (47) and may enhance myeloid-derived suppressor cells (MDSC)-mediated T cell activation (48,49). When analyzing the TME infiltration of HPI and immune cells, we found that eosinophils, macrophages, neutrophils, helper T cells, and Th1 cells had higher infiltrations in the high HPI group. Previous studies have shown that the main immune cells associated with hypoxia in OSCC are Macrophages (50,51) and Neutrophils (50,52). Eosinophils, T helper cells, and Th1 cells have also been reported in OSCC (53,54), but they are also related to hypoxia, our first discovery. Together, the increase in other hypoxia-related phenotypes (TNFA signaling via NFKB, TGF BETA signaling, Epithelial Mesenchymal Transition, Glycolysis and Angiogenesis) is also related to increased immune cell infiltration. However, we found that the NSD1 gene mutation in OSCC was significantly related to the hypoxia index. The NSD1 mutation has shown a close relationship with tumor malignancy in colorectal cancer (55) and head and neck cancer (56,57) and is a possible marker of a poor prognosis of tumors. However, there have been no studies in OSCC that describe the relationship between hypoxia and NSD1 mutations. Finally, we constructed a survival prediction model for hypoxia-related genes in OSCC, and determined the relationship between six hypoxia-related genes (PGK1, JMJD6, S100A4, SLC2A3, DDIT4 and HK1) and prognosis, and evaluated them with a calibration curve. Tumor cells tend to produce ATP through aerobic glycolysis, rather than oxidative phosphorylation (Warburg Effect). Phosphoglycerate Kinase 1 (PGK1) is the key rate limiting enzyme in the aerobic glycolysis process (58) and has been verified as hypoxia-related OSCC prognostic genes (59). Jumonji Domain Containing 6 (JMJD6) inhibition may be the way to treat the oral cavity (60). Calcium binding protein A4 (61) (S100A4), Solute Carrier Family 2 Member 3 (62) (SLC2A3), hexokinase 1 (63) (HK1) in OSCC are related to a poor prognosis of patients. DNA Damage Inducible Transcript 4 (DDIT4) is the first to be related to the poor prognosis of OSCC tumors.

Conclusions

In general, we rely on hypoxia-related HPI to explore the immunity, gene mutations, and patient prognosis of OSCC patients. In addition to determining the TME infiltration of immune cells (Eosinophils, Macrophages, Neutrophils, T cells and Th1 cells) associated with HPI, we also constructed a prediction model based on the six best hypoxia-related genes (PGK1, JMJD6, S100A4, SLC2A3, DDIT4, and HK1).

Acknowledgments

We thank the TCGA program and GEO database for the RNA-sequence and clinical data of patients with OSCC.

Funding: None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://dx.doi.org/10.21037/atm-21-4990

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://dx.doi.org/10.21037/atm-21-4990). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Warnakulasuriya S. Global epidemiology of oral and oropharyngeal cancer. Oral Oncol 2009;45:309-16. [Crossref] [PubMed]
2. Argiris A, Karamouzis MV, Raben D, et al. Head and neck cancer. Lancet 2008;371:1695-709. [Crossref] [PubMed]
3. Gottesman MM. Mechanisms of cancer drug resistance. Annu Rev Med 2002;53:615-27. [Crossref] [PubMed]
4. Shao C, Yang F, Miao S, et al. Role of hypoxia-induced exosomes in tumor biology. Mol Cancer 2018;17:120. [Crossref] [PubMed]
5. Brown JM, Giaccia AJ. The unique physiology of solid tumors: opportunities (and problems) for cancer therapy. Cancer Res 1998;58:1408-16. [PubMed]
6. Harris AL. Hypoxia--a key regulatory factor in tumour growth. Nat Rev Cancer 2002;2:38-47. [Crossref] [PubMed]
7. Roma-Rodrigues C, Mendes R, Baptista PV, et al. Targeting Tumor Microenvironment for Cancer Therapy. Int J Mol Sci 2019;20:840. [Crossref] [PubMed]
8. Multhoff G, Vaupel P. Hypoxia Compromises Anti-Cancer Immune Responses. Adv Exp Med Biol 2020;1232:131-43. [Crossref] [PubMed]
9. Facciabene A, Peng X, Hagemann IS, et al. Tumour hypoxia promotes tolerance and angiogenesis via CCL28 and T(reg) cells. Nature 2011;475:226-30. [Crossref] [PubMed]
10. Seiwert TY, Burtness B, Mehra R, et al. Safety and clinical activity of pembrolizumab for treatment of recurrent or metastatic squamous cell carcinoma of the head and neck (KEYNOTE-012): an open-label, multicentre, phase 1b trial. Lancet Oncol 2016;17:956-65. [Crossref] [PubMed]
11. Ferris RL, Blumenschein G Jr, Fayette J, et al. Nivolumab for Recurrent Squamous-Cell Carcinoma of the Head and Neck. N Engl J Med 2016;375:1856-67. [Crossref] [PubMed]
12. Brooks JM, Menezes AN, Ibrahim M, et al. Development and Validation of a Combined Hypoxia and Immune Prognostic Classifier for Head and Neck Cancer. Clin Cancer Res 2019;25:5315-28. [Crossref] [PubMed]
13. Joseph JP, Harishankar MK, Pillai AA, et al. Hypoxia induced EMT: A review on the mechanism of tumor progression and metastasis in OSCC. Oral Oncol 2018;80:23-32. [Crossref] [PubMed]
14. Pérez-Sayáns M, Suárez-Peñaranda JM, Pilar GD, et al. Hypoxia-inducible factors in OSCC. Cancer Lett 2011;313:1-8. [Crossref] [PubMed]
15. Cancer Genome Atlas Research Network. The Cancer Genome Atlas Pan-Cancer analysis project. Nat Genet 2013;45:1113-20. [Crossref] [PubMed]
16. Akbani R, Ng PK, Werner HM, et al. A pan-cancer proteomic perspective on The Cancer Genome Atlas. Nat Commun 2014;5:3887. [Crossref] [PubMed]
17. Goldman M, Craft B, Swatloski T, et al. The UCSC Cancer Genomics Browser: update 2015. Nucleic Acids Res 2015;43:D812-7. [Crossref] [PubMed]
18. Mitchell S, Vargas J, Hoffmann A. Signaling via the NFκB system. Wiley Interdiscip Rev Syst Biol Med 2016;8:227-41. [Crossref] [PubMed]
19. Colak S, Ten Dijke P. Targeting TGF-β Signaling in Cancer. Trends Cancer 2017;3:56-71. [Crossref] [PubMed]
20. Xie F, Ling L, van Dam H, et al. TGF-β signaling in cancer metastasis. Acta Biochim Biophys Sin (Shanghai) 2018;50:121-32. [Crossref] [PubMed]
21. Zhao M, Mishra L, Deng CX. The role of TGF-β/SMAD4 signaling in cancer. Int J Biol Sci 2018;14:111-23. [Crossref] [PubMed]
22. Dongre A, Weinberg RA. New insights into the mechanisms of epithelial-mesenchymal transition and implications for cancer. Nat Rev Mol Cell Biol 2019;20:69-84. [Crossref] [PubMed]
23. Zhang Y, Weinberg RA. Epithelial-to-mesenchymal transition in cancer: complexity and opportunities. Front Med 2018;12:361-73. [Crossref] [PubMed]
24. Akram M. Mini-review on glycolysis and cancer. J Cancer Educ 2013;28:454-7. [Crossref] [PubMed]
25. Abbaszadeh Z, Çeşmeli S, Biray Avcı Ç. Crucial players in glycolysis: Cancer progress. Gene 2020;726:144158 [Crossref] [PubMed]
26. Soleimani A, Rahmani F, Ferns GA, et al. Role of the NF-κB signaling pathway in the pathogenesis of colorectal cancer. Gene 2020;726:144132 [Crossref] [PubMed]
27. Park JA, Kwon YG. Hippo-YAP/TAZ signaling in angiogenesis. BMB Rep 2018;51:157-62. [Crossref] [PubMed]
28. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
29. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
30. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
31. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
32. Bindea G, Mlecnik B, Tosolini M, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013;39:782-95. [Crossref] [PubMed]
33. Shu S, Wang Y, Zheng M, et al. Hypoxia and Hypoxia-Inducible Factors in Kidney Injury and Repair. Cells 2019;8:207. [Crossref] [PubMed]
34. Darby IA, Hewitson TD. Hypoxia in tissue repair and fibrosis. Cell Tissue Res 2016;365:553-62. [Crossref] [PubMed]
35. Moen I, Stuhr LE. Hyperbaric oxygen therapy and cancer--a review. Target Oncol 2012;7:233-42. [Crossref] [PubMed]
36. Kaur RP, Vasudeva K, Kumar R, et al. Role of p53 Gene in Breast Cancer: Focus on Mutation Spectrum and Therapeutic Strategies. Curr Pharm Des 2018;24:3566-75. [Crossref] [PubMed]
37. Kim S, Jeong S. Mutation Hotspots in the β-Catenin Gene: Lessons from the Human Cancer Genome Databases. Mol Cells 2019;42:8-16. [PubMed]
38. Zhang L, Meng X, Zhu XW, et al. Long non-coding RNAs in Oral squamous cell carcinoma: biologic function, mechanisms and clinical implications. Mol Cancer 2019;18:102. [Crossref] [PubMed]
39. Rosebush MS, Rao SK, Samant S, et al. Oral cancer: enduring characteristics and emerging trends. J Tenn Dent Assoc 2011;91:24-7; quiz 28-9. [PubMed]
40. Ueda S, Saeki T, Osaki A, et al. Bevacizumab Induces Acute Hypoxia and Cancer Progression in Patients with Refractory Breast Cancer: Multimodal Functional Imaging and Multiplex Cytokine Analysis. Clin Cancer Res 2017;23:5769-78. [Crossref] [PubMed]
41. Chen A, Sceneay J, Gödde N, et al. Intermittent hypoxia induces a metastatic phenotype in breast cancer. Oncogene 2018;37:4214-25. [Crossref] [PubMed]
42. Lin SH, Koong AC. Breathing New Life Into Hypoxia-Targeted Therapies for Non-Small Cell Lung Cancer. J Natl Cancer Inst 2018; [Crossref] [PubMed]
43. Daponte A, Ioannou M, Mylonis I, et al. Prognostic significance of Hypoxia-Inducible Factor 1 alpha(HIF-1 alpha) expression in serous ovarian cancer: an immunohistochemical study. BMC Cancer 2008;8:335. [Crossref] [PubMed]
44. Eckert AW, Kappler M, Schubert J, et al. Correlation of expression of hypoxia-related proteins with prognosis in oral squamous cell carcinoma patients. Oral Maxillofac Surg 2012;16:189-96. [Crossref] [PubMed]
45. Heo K, Kim YH, Sung HJ, et al. Hypoxia-induced up-regulation of apelin is associated with a poor prognosis in oral squamous cell carcinoma patients. Oral Oncol 2012;48:500-6. [Crossref] [PubMed]
46. Nagarajah NS, Vigneswaran N, Zacharias W. Hypoxia-mediated apoptosis in oral carcinoma cells occurs via two independent pathways. Mol Cancer 2004;3:38. [Crossref] [PubMed]
47. Chouaib S, Noman MZ, Kosmatopoulos K, et al. Hypoxic stress: obstacles and opportunities for innovative immunotherapy of cancer. Oncogene 2017;36:439-45. [Crossref] [PubMed]
48. Noman MZ, Desantis G, Janji B, et al. PD-L1 is a novel direct target of HIF-1α, and its blockade under hypoxia enhanced MDSC-mediated T cell activation. J Exp Med 2014;211:781-90. [Crossref] [PubMed]
49. Barsoum IB, Smallwood CA, Siemens DR, et al. A mechanism of hypoxia-mediated escape from adaptive immunity in cancer cells. Cancer Res 2014;74:665-74. [Crossref] [PubMed]
50. Mughees M, Sengupta A, Khowal S, et al. Mechanism of tumour microenvironment in the progression and development of oral cancer. Mol Biol Rep 2021;48:1773-86. [Crossref] [PubMed]
51. Okubo M, Kioi M, Nakashima H, et al. M2-polarized macrophages contribute to neovasculogenesis, leading to relapse of oral cancer following radiation. Sci Rep 2016;6:27548. [Crossref] [PubMed]
52. Siddiqui S, Singh A, Faizi N, et al. Cell cannibalism in oral cancer: A sign of aggressiveness, de-evolution, and retroversion of multicellularity. J Cancer Res Ther 2019;15:631-7. [Crossref] [PubMed]
53. Stasikowska-Kanicka O, Wągrowska-Danilewicz M, Danilewicz M. T cells are involved in the induction of macrophage phenotypes in oral leukoplakia and squamous cell carcinoma-a preliminary report. J Oral Pathol Med 2018;47:136-43. [Crossref] [PubMed]
54. Deepthi G, Kulkarni PG, Nandan K SR. Eosinophils: An imperative histopathological prognostic indicator for oral squamous cell carcinoma. J Oral Maxillofac Pathol 2019;23:307. [Crossref] [PubMed]
55. Zhunussova G, Afonin G, Abdikerim S, et al. Mutation Spectrum of Cancer-Associated Genes in Patients With Early Onset of Colorectal Cancer. Front Oncol 2019;9:673. [Crossref] [PubMed]
56. Cancer Genome Atlas Network. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature 2015;517:576-82. [Crossref] [PubMed]
57. Pan C, Izreig S, Yarbrough WG, et al. NSD1 mutations by HPV status in head and neck cancer: differences in survival and response to DNA-damaging agents. Cancers Head Neck 2019;4:3. [Crossref] [PubMed]
58. McCarrey JR, Thomas K. Human testis-specific PGK gene lacks introns and possesses characteristics of a processed gene. Nature 1987;326:501-5. [Crossref] [PubMed]
59. Zhao C, Zhou Y, Ma H, et al. A four-hypoxia-genes-based prognostic signature for oral squamous cell carcinoma. BMC Oral Health 2021;21:232. [Crossref] [PubMed]
60. Lee CR, Lee SH, Rigas NK, et al. Elevated expression of JMJD6 is associated with oral carcinogenesis and maintains cancer stemness properties. Carcinogenesis 2016;37:119-28. [Crossref] [PubMed]
61. Hu FW, Lee SS, Yang LC, et al. Knockdown of S100A4 impairs arecoline-induced invasiveness of oral squamous cell carcinomas. Oral Oncol 2015;51:690-7. [Crossref] [PubMed]
62. Feitosa SG, Viana KF, Luna ECM, et al. Immunohistochemical Evaluation of GLUT-3 and GLUT-4 in Oral Epithelial Dysplasia and Oral Squamous Cell Carcinoma. Asian Pac J Cancer Prev 2018;19:1779-83. [PubMed]
63. Yamada T, Uchida M, Kwang-Lee K, et al. Correlation of metabolism/hypoxia markers and fluorodeoxyglucose uptake in oral squamous cell carcinomas. Oral Surg Oral Med Oral Pathol Oral Radiol 2012;113:464-71. [Crossref] [PubMed]

(English Language Editor: J. Chapnick)