From the OSCC database (http://www.cancergenome. nih.gov/), TCGA RNA-seq dataset and the respective clinical data of 257 oral squamous cell carcinoma samples (e.g., palate samples, gingival samples, tongue root samples, unspecified oral samples, mouth floor samples, unspecified tongue samples, and others) were downloaded (Zollinger et al., 2018; Chen et al., 2020). The baseline and clinical data, including sex, age, tumor grade, race, survival status, survival time and pathologic stage, were collected. The login numbers of all specimens are listed in Supplementary file S1. We also searched the human expression profile of oral cancer from NCBI GEO (accession number: GSE41613). The public domain names and direct Web links are summarized in Supplementary File S2.

CIBERSORT is a deconvolution method applied for the characterization of cell composition of a complex tissue (Peng et al., 2021). LM22 is a leukocyte gene signature matrix that contains the expression levels of 547 genes for distinguishing 22 human hematopoietic cell phenotypes (Mills et al., 2021). In this study, CIBERSORT combined with LM22 signature matrix were used to quantify the infiltration levels of twenty-two immune cells in OSCC. Subsequently, consensus clustering was performed using the ConsensuClusterPlus R package (Yang et al., 2017). The cumulative distribution function curve of a consensus matrix was used to determine the optimal number of clusters. Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression (ESTIMATE) data were then employed to assess the stromal score and immune score of each OSCC sample. The Kaplan-Meier survival curve was performed for different TIME clusters using survminer R package.

OSCC patients were stratified into three TIME clusters on the basis of tumor immune cell infiltration. To filter the differentially expressed genes (DEGs) among different TIME subtypes, the genes related to TIME clusters were identified using the limma R package (Brorson et al., 2020). The threshold values of adjusted p-value < 0.05 and |log2 (fold change) | >1 were applied.

First, consensus clustering was used to categorize the OSCC patients into different TIME gene clusters based on the expression of all DEGs by ConsensuClusterPlus R package. Then, the prognostic value of ferroptosis-related genes was computed by K-M analyses. p-value < 0.05 was deemed as statistically significant. Next, Boruta algorithm was employed to evaluate the dimensionality reduction of the distinct TIME gene signatures, and principal component 1 was extracted as the signature score via principal component analysis. Finally, a gene expression grade index was applied to measure the TIME score of each specimen: TIME score = ∑ PC 1 x − ∑ PC 1 y

The genetic alteration data of OCSC patients were derived from the TCGA data portal (https://www.cancer.gov/tcga/) and applied as the training dataset. We identified the OSCC driver genes by maftools in R software and calculated the sum of somatic alterations in OSCC driver genes of each sample as the TMB score (Mayakonda et al., 2018). We also estimated the correlation between the TIME scores and the TMB scores (Zhang L et al., 2020). Next, the Kaplan-Meier survival curve for overall survival (OS) was constructed based on the high and low TMB score patterns using survminer R package. Another Kaplan-Meier survival curve was generated among the four different groups (high TMB score and high TIME score; high TMB score and low TIME score; low TMB score and high TIME score; low TMB score and low TIME score).

TIMER is a novel algorithm used for characterizing the cell composition of a complex tissue and to estimate the level of tumor-infiltrating immune cells according to the corresponding gene expression profiles (Li et al., 2020). In this study, TIMER method was used to predict the tumor purity and immune score of each OSCC specimen from two TIME score clusters and to verify the correlation between the tumor-infiltrating immune cell level and the immune checkpoint gene expression level (CTLA4, PD-L1 and PD-1) estimated by CIBERSORT method (Wang et al., 2020).

GSEA database (http://software.broadinstitute.org/gsea/index.jsp) was employed to assess the difference in potential biological processes between high-risk and low-risk groups (p < 0.05).

All statistical tests were bilateral, and were carried out using R version 3.6.1 and Bioconductor (https://www.bioconductor.org/). The levels of significance (p-value) were set as <0.05, <0.01 < 0.001 and <0.0001.

Firstly, the enrichment levels of TIME-infiltrating cells in OSCC tumor tissues were quantify using the CIBERSORT and ESTIMATE algorithms (Schelker et al., 2017; Yoshihara et al., 2013). Unsupervised clustering of 304 tumor samples from the meta-cohort (GSE41613 and tumor genome map [TCGA]) matched with tumor immune microenvironment (TIME) was carried out using the ConsesusClusterPlus package in R software, and the OSCC patients were divided into different subtypes through consensus clustering (Figures 1A–C; Supplementary File S3). (Gan et al., 2018) To compare the immune cell compositions of the TIME, we performed a heatmap for the three clusters (Figure 1D). Among the three distinct subtypes of TIME cell infiltration, the TIME cluster A exhibited a high level of M0 macrophage infiltration; the TIME cluster B demonstrated a significant increase in the amount of resting memory CD4⁺ T cells; and the TIME cluster C showed significant increases in the infiltration levels of CD8⁺ T cells, Follicular helper T cells and activated memory CD4⁺ T cells. The clinical characteristics of the three TIME clusters are presented in Table 1. Notably, OSCC patients in cluster C had a longer survival time (Log-rank test/Log rank, p = 0.048) than those in cluster A and B (Figure 1E), suggesting that TIME is associated with the prognosis of diffuse OSCC patients. Subsequently, the correlation coefficient heatmap was constructed, and the findings revealed that the infiltration of CD8⁺ T cells was positively correlated with the infiltration of CD4⁺ T cells (Figure 1G). Furthermore, we also analyzed the expression levels of the crucial immune checkpoints PD-1 and CTLA-4 in each TIME subtype (Goltz et al., 2018; Gabrych et al., 2019). The patients in TIME cluster C had significantly higher expression levels of CTLA-4 and PD-1 compared to cluster A and B.

To explore differentially expressed genes among distinct immunophenotypes, limma R package was employed to compare the differences in mRNA expression levels among different subtypes. As a result, a total of 360 DEGs were identified in all OSCC samples. Then, an unsupervised clustering analysis was conducted on these 360 DEGs (Supplementary Files S4, S5), which classified the samples into two genome clusters, namely, gene clusters A and B (Figures 2A–C; Supplementary File S6). The 360 gene signatures were positively related to the gene cluster (TIME gene signature A), and the remaining were designated as TIME gene signature B (Supplementary File S7) (Kaymak et al., 2021). The expression of these DEGs is visualized in a heatmap (Figure 2D) by using the clusterProfiler R package (Zeng et al., 2019).

To further evaluate the prognosis of OSCC patients with different TIME gene clusters, the survival analysis in R programming language with a threshold of log-rank<0.05. It was observed that the patients in gene cluster B exhibited a better prognosis, whereas those in gene cluster A had a poor outcome (log rank test, p = 0.007; Figure 2E).

GO function analysis was carried out for TIME signature genes A and B respectively (Figures 2F,G; Supplementary Files S8, S9). Interestingly, the results showed that gene cluster B was markedly associated with higher immune scores (Figure 2H). The immune scores in the two TIME gene clusters were in accordance with the distribution patterns of TIME-infiltrating cells. The gene cluster A displayed a remarkable increase in the filtration of M0 macrophages, resting mast cells and neutrophils, while the gene cluster B exhibited a higher infiltration of activated CD4⁺ memory T cells, CD8⁺ T cells and M1 macrophages.

Additionally, the two TIME signature gene clusters demonstrated remarkable difference in the expression levels of PD-1 and CTLA-4 (Figure 2I). The PD-1/CTLA-4 expression level of TIME gene cluster B was relatively high, while that of TIME gene cluster A was slightly lower.

Principal component analysis was conducted to measure the TIME score of each sample, that is, the sum of the TIME score A and from TIME signature gene A and B, respectively, in order to obtain the quantitative index of TIME landscape (Supplementary File S10). The patients were stratified into high and low TIME score groups according to the optimal cut-off values obtained from survminer R package.

Figure 3A shows the distribution patterns of OSCC patients in the two gene clusters. The prognostic implication of TIME scores was evaluated, and the immune tolerance conditions of each group in the TCGA cohort were assessed before analyzing the prognostic values of TIME score in the TCGA cohort and other independent datasets. The immune checkpoint-related gene signatures of LAG3, HAVCR2, CD27, PDCD1LG2, PDCD1, CD274, TIGIT, CTLA4, ICOS and IDO1 were selected (Song et al., 2021). It was found that all these immune checkpoint-related genes were markedly overexpressed in high TIME score group, as revealed by the Wilcoxon test (Figure 3B). TIME score subtype analysis in the Kaplan-Meier plotter indicated that OSCC patients in high TIME score group exhibited markedly better survival rate than those in low TIME score group (log rank test, p = 0.017; Figure 3C). In addition, the results of GSEA demonstrated that the signaling pathways of B cell receptor, nature killer cell-mediated cytotoxicity, and toll-like receptor were significantly enriched in high TIME score group (Figure 3D).

TMB is an important molecular marker for assessing the outcomes of tumor immunotherapy (Mayakonda et al., 2018). Previous evidence shows that the level of TMB can not only predict the efficacy of immunotherapy, but also accurately predict the efficacy of many targeted and chemotherapeutic drugs (Zhu et al., 2019b). Besides, TMB also has important clinical significance. In this study, we attempted to determine the relationship between TMB and TIME scores. First, the TMBs of OSCC patients with high and low TIME scores were compared (Supplementary File S11). As shown in Figure 4A, the TMB of OSCC patients was higher in low TIME score group than in high TIME score group (Wilcoxon test, p = 0.015). Next, we categorized the patients into high and low TMB groups. As shown in Figure 4B, OSCC patients in low TMB group exhibited better survival rate than those in high TMB group (log-rank test, p < 0.001). Based on these contradictory results, we further explored the combined effect of TMB and TIME scores on OSCC prognosis. Stratified survival analysis showed that TMB scores did not affect TIME scores-based prognostic prediction. However, a significant difference in survival rates was observed between high and low TMB groups stratified by TIME scores (Figure 4C).

In addition, the distribution patterns of somatic alterations in OSCC patients were compared between high and low TIME groups. The R package maftools was used to identify potential driver genes in OSCC patients. The top 20 variant mutated driver genes were selected for further analysis. The annotated mutation files were downloaded from TCGA database and then analyzed. As shown in Figures 4D,E and Table2, the frequencies of NSD1, CASP8 and TP53 mutations were noticeably different between high and low TIME score groups. These results may provide new insights into the mechanisms of gene mutations and TIME compositions as well as novel targets for immunotherapy.

To further evaluate the prognostic implication of TIME scores on OSCC patients, all samples were categorized into different groups by clinical characteristics (e.g., gender, stage and age). The survival curves of TIME scores were constructed using R software with survival package. As shown in Figures 5A,B, TIME scores were positively correlated with OS among patients younger than 60. However, no marked difference was found between high and low score groups with age over 60. The TIME scores were also positively correlated with OS among male patients (Figures 5C,D) and stage Ⅲ/Ⅳ patients (Figures 5E,F).