The Cancer Genome Atlas (TCGA) mRNA sequence data [htseq-FPKM in log2(x + 1) transformed] together with clinical information of OSCC were obtained from the UCSC Xena browser (GDC hub: https://gdc.xenahubs.net). For validation, microarray profiles of GSE41613 containing clinical annotations were extracted by GEOquery R package. The mentioned clinical traits are demonstrated in Table 1 . The batch effects normalized mRNA data of pancancer with clinical information were downloaded from UCSC Xena browser. The hypoxia gene set containing 200 classical hypoxia-associated genes was obtained from gene set enrichment analysis (GSEA) (http://www.gsea-msigdb.org/). Expression data of OSCC cell lines [TPM in log2(x + 1) transformed) were downloaded from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) (22). Drug sensitivity data (area under the curve—AUC) of OSCC cells from the Cancer Therapeutics Response Portal (CTRP v.2.0) and PRISM Repurposing dataset (19Q4) were acquired from the dependency map (DepMap) portal (https://depmap.org/portal/). The ICB treatment cohort GSE91061 (23) was downloaded from the GEO database [FPKM in log2(x + 1)] transformed and used for subsequent validation. The CheckMate 009 (CM-009), CheckMate 010 (CM-010), and CheckMate 025 (CM-025) (24) were combined together to investigate the significance of our risk score [FPKM in log2(x + 1)]. We also downloaded RNA-seq (count values) data of IMvigor210 cohort (25) with clinical information by the “IMvigor210CoreBiologies” R package and transformed it into FPKM values. The log2(FPKM + 1) was calculated on expression data for further comparison.

Unsupervised clustering was applied to recognize different hypoxia patterns and classify OSCC patients for further analysis. A consensus hierarchical clustering algorithm based on the expression of 34 prognostic hypoxia genes was conducted by the “ConsensuClusterPlus” R package with Euclidean distance and Ward.D2’s linkage (number of bootstraps=50, item subsampling proportion = 0.8, feature subsampling proportion = 0.8).

Univariate Cox regression analysis was conducted to identify prognostic hypoxia genes and clinical events. Multivariate Cox regression analysis was performed to recognize independent prognostic factors. The Kaplan–Meier survival curve was applied to analyze the prognostic significance between distinct groups.

The hypoxia-associated gene sets were downloaded from GSEA. The single-sample gene set enrichment analysis (ssGSEA) algorithm in “GSVA” R package was conducted to calculate the hypoxia score of each OSCC patient.

The MAF file of OSCC containing the detailed mutation information of the training set was downloaded from TCGA (https://portal.gdc.cancer.gov/) and further processed. The “maftool” R package was performed to analyze gene mutant features between two OSCC subclasses.

The “Limma” R package was applied to identify differential genes between two clusters with a standard of |log FC| >1.2 and adjusted P-value <0.05. Further gene ontology (GO) function enrichment of selected genes was performed by ClueGO in Cytoscape.

The immune score and the tumor purity were calculated by the ESTIMATE algorithm (26). The CIBERSORT algorithm was applied to evaluate the LM22 gene signatures in OSCC subtypes (27). What is more, the Epic algorithm was also used to calculate the contents of immune cell infiltration in the microenvironment (28).

k-Nearest neighbor (k-NN) imputation was performed to impute the missing AUC values of the CTRP and PRISM datasets. Before imputation, drugs with more than 20% of missing data were excluded. Furthermore, the “pRRophetic” R package was performed to measure the AUC values of samples by ridge regression.

Considering the difference of each platform, before developing or validating the risk score, we conducted z-scale of the mRNA data in each platform (TCGA, GSE41613, GSE91061, CM cohorts, and IMvigor210). Then, the “glmnet” R package was performed to filter the prognosis-related hypoxia genes by LASSO Cox regression analysis with a 10-fold cross-validation. After identifying the significant genes, their regression coefficients (β) were estimated by multivariate Cox regression via LASSO, and we calculated the risk score of each OSCC patient by the formula as follows:

Univariate Cox and multivariate Cox regression analyses of some clinical traits were first performed and finally determined a sum of four independent prognostic factors for further establishment. Afterward, a nomogram with the four factors was developed for predicting 1- and 3-year OS of OSCC patients. The calibration plot was performed to estimate the accuracy and consistency of the prognostic models. Survival net benefits of each variable were estimated with decision curve analysis (DCA) by “stdca.R.”

Principal components analysis (PCA) was applied to verify the hypoxia patterns of different subtypes. Potential ICB response was predicted by the tumor immune dysfunction and exclusion (TIDE) algorithm (29). The “upsetR” R package was used to visualize the intersections between promising agents in different subtypes.

R 4.0.2 (https://www.r-project.org/) was mainly used for statistical analysis. Student’s t-test or one-way analysis of variance was used to analyze differences between groups in variables with a normal distribution. Categorical variables between two groups were compared using chi-square test. A two sided P-value <0.05 was considered statistically significant.

As depicted in Figure 1A , a brief flowchart was demonstrated to introduce our study. Considering the critical role of hypoxia condition in the tumor microenvironment, we summarized a sum of 188 classical hypoxia-stimulated genes available from GSEA and estimated their prognostic value for further classification ( Table S1 ). Univariate Cox proportional hazards model was conducted and finally filtered 34 genes with significant risks on survival of patients in the training set ( Figures S1A, B ). Hence, based on the expression similarity of the 34 hypoxia-related gene signature, the consensus clustering method was used to cluster the samples. We selected k = 2 as the optimal number of clusters, which could divide all samples into two groups with less correlation between groups in the training and testing cohorts ( Figures 1B, C ). Then, PCA was conducted to compare the transcriptional profile between these two clusters in the two cohorts, suggesting a significant distinction between these two subgroups ( Figures 1D, E ). In order to evaluate the clinical relevance of this clustering, the survival analysis between the two subclasses was conducted. In these two sets, cluster 2 was consistently associated with worse prognosis, highlighting the potential clinical utility of this hypoxia-associated subtyping ( Figures 1F, G ).

To better understand the hypoxia status of the two clusters, we conducted the ssGSEA algorithm to calculate the scores of some hypoxia-associated processes. As expected, patients in cluster 2 were enriched in higher hypoxia condition in the training and testing cohorts ( Figure 2A ). What is more, a total of nine hypoxia-associated key genes were also verified to be highly expressed in cluster 2, which was consistent with the aforesaid ssGSEA result ( Figure 2B ). Hence, we could define cluster 2 as a “high hypoxia subclass” compared with cluster 1.

Recent studies have reported the hypoxia phenotype associated with gene mutations (30). We further investigated the difference of gene mutations among these two clusters. As illustrated in the waterfall plot, differently mutated genes were detected between the two clusters and GNPTAB was finally identified as the most differentially highly mutated gene in cluster 2 ( Figure 3A ) (P < 0.01). Furthermore, based on the oncodriveCLUST algorithm, we predicted HRAS as the driver gene candidate in cluster 1 and MAST4 in cluster 2 ( Figure 3B ). What is more, tumor mutational burden (TMB) was significantly increased in cluster 2 ( Figure 3C ).

To obtain deeper insights into the molecular characteristics of the two OSCC clusters, we conducted the differentially expressed genes (DEGs) analysis and their GO analysis in the training dataset. With a threshold of |log2 FC| >1.2 and adjusted P-value <0.05, a sum of 55 DEGs were identified for the two clusters. The expressions of DEG between these two clusters were demonstrated by a heatmap ( Figure S2A ). GO analysis based on Cytoscape showed that the cluster-specific genes were significantly enriched in immune cell infiltration, suggesting a distinct immune difference between these two clusters ( Figure S2B ).

To reveal the difference of these two clusters on the tumor microenvironment, we first calculated the immune score and tumor purity both in the training and testing sets based on the ESTIMATE algorithm. We found that the immune score was decreased and purity score was elevated in cluster 2 compared with cluster 1 ( Figures 4A and S3A ). With the significant difference in immune score and purity score identified between clusters, we further compared the relative ratio of 22 kinds of immune cells by the CIBERSORT algorithm. There existed six immune cell populations significantly differently enriched between the two clusters in the training set and nine immune cells in the testing set ( Figures 4B and S3B ). Combined, macrophages M0, activated mast cells, were enriched in cluster 2, while CD8 T cells, resting mast cells, were deficient in both two sets. We further conducted the Epic algorithm to validate our results and found that only CD8 T cells were consistently lacking in cluster 2 in the two cohorts ( Figures 4C and S3C ). CD8 T cell, also known as cytotoxic T cell (CTL), exerted a critical role in antitumor immunity. We further examined two indicators of T-cell killing ability between the two clusters. Similarly, cluster 2 also exhibited lower CYT score and IFNG expression than cluster 1 in the training set and testing set, which was consistent with previous studies that showed an association between high CYT levels and higher patient OS ( Figures 4D and S3D ). Taken together, it was the lower composition of CD8 T cells and their disability of killing tumor cells that led a worse prognosis in cluster 2.

After investigating the distinct molecular and biological characters between these two clusters, we sought to explore specific treatment options for each cluster. Considering the vital role of CD8 T cells in immunotherapy and their significant differences between the two clusters, we further assessed their immunotherapy response based on the TIDE method. In both training set and testing set, the TIDE score was significantly lower in cluster 1 compared with cluster 2, indicating patients in cluster 1 might be more sensitive to ICB therapy ( Figure 5A ). For cluster 2 patients, we hoped to seek for traditional chemotherapeutics to achieve targeted therapy. After the filtering procedure described in the Material and Methods, we finally obtained 16 OSCC cells with 913 drugs in the PRISM and 22 OSCC cells with 465 drugs in the CTRP dataset. The pRRophetic package with a built-in ridge regression model was then applied to predict the drug response of clinical samples in the training set based on their expression profiles, and the estimated AUC value of each compound in each sample was thus obtained. We finally identified four agents simultaneously with lower AUC values in cluster 2 both in the PRISM- and CTRP-predicted datasets ( Figures 5B, C and S4 ). To further filter a more therapeutically significant drug in OSCC, we took their clinical phase and experimental evidence from the literature into account. Finally, we identified only bortezomib as the optimal drug that has the potential for cluster 2 treatment ( Figure 5D ).

To establish a signature for clinical implications, it is of great significance to filter the most representative genes of each cluster. Considering HIF1A serving as the key transcription factor in hypoxia, we intersected the DEGs between the two clusters with 4,748 potential targets of HIF1A in OSCC and found a sum of 6 candidate genes in the intersection ( Figure 6A ), identified as “Clustering-specific hypoxia-related genes.” To obtain the most powerful prognostic markers, the LASSO Cox regression analysis was conducted ( Figure 6B ). A total of five gene signatures were generated and the coefficients were estimated by multivariate Cox regression via LASSO ( Table S2 ). There existed a transcriptional difference between the two clusters ( Figure 6C ). After calculating the risk scores of the signature based on the regression coefficients, we intriguingly found that cluster 2 possessed a higher score in the two cohorts ( Figures 6D, E ). Further survival analysis revealed that patients in the high-score group exhibited significantly worse prognosis than OSCC patients or cluster 1 patients with low-score ( Figures 6F, G ). Although there was no significant survival difference between high and low scores in cluster 2 in the training set (P = 0.1) and testing set (P = 0.13), it was still obvious that a high hypoxia score was associated with the tendency toward worse prognosis ( Figures 6F, G ). The results were consistent with the above data that cluster 2 conferred the poorer prognosis. In order to determine the prognostic significance of the signature in other organ sites, we conducted the survival analysis of our hypoxia score across 33 TCGA cancer types. Similarly, the hypoxia risk score also served as an unfavorable prognostic biomarker for pancancer ( Figure 6H ). What is more, the predicted AUC values of bortezomib from CTRP and PRISM were also decreased in the high hypoxia score group, validating its promising clinical value for high-risk OSCC patients ( Figures 6I, J ).

To verify whether the hypoxia-related signature was an independent prognostic factor, univariate and multivariate Cox regression analyses were conducted ( Figures 7A, B ). The results in univariate Cox regression revealed that risk score, age, and angiolymphatic and perineural invasion had a significant association with the OS of OSCC patients. In multivariate Cox regression, risk score, age, and angiolymphatic and perineural invasion were identified as independent prognostic factors of OSCC. Then, we applied these four independent factors to establish a nomogram for predicting OSCC 1- and 3-year OS ( Figure 7C ). With the score increasing, the OS of patients decreased. Moreover, the calibration plots at 1 and 3 years approached 45 degrees, indicating a great performance of the established nomogram ( Figure 7D ). Meanwhile, DCA was performed to compare the clinical usability and benefits of the nomogram with that of the age and angiolymphatic and perineural invasion. As shown in Figure 7E , compared with age and angiolymphatic and perineural invasion, the 1-year DCA curves of the new nomogram showed larger net benefits across a range of death risk.

Immunotherapy has been proven relevant to improve survival in the treatment of multiple tumor types. Thus, identification of patients who will benefit most from ICB treatment is of great necessity. Our analysis revealed that the TIDE was significantly increased in the high hypoxia score group, indicating its crucial role in regulating immune response ( Figure 8A ). Based on three immunotherapy cohorts, we identified that patients with a high hypoxia score group always exhibited clinical disadvantages and markedly shortened survival (P = 0.026 in GSE91061, P = 0.039 in CM009+010+025 cohorts, and P = 0.029 in IMvigor210) ( Figures 8B, C, E ). In CM009+010+025 cohorts, the chi-squared test conducted between low and high hypoxia score groups demonstrated significantly better therapeutic outcomes in low score patients ( Figure 8D ). Similarly, patients with high hypoxia scores exhibited less treatment effectiveness in the IMvigor210 cohort ( Figure 8F ). We also compared the hypoxia score levels in the three immune subtypes of IMvigor210. The immune-inflamed subtype showed significantly the lowest risk score, which further confirmed our analysis above ( Figure 8G ). In addition, TMB was significantly decreased in the high-score group ( Figure 8H ). In all, our results strongly suggested that hypoxia score was associated with the response to immunotherapy and could further effectively predict the prognosis of patients.