Normalized transcriptome data and clinical information of OSCC were retrieved from the Cancer Genome Atlas (TCGA) database via the UCSC Xena (https://tcga.xenahubs.net) on 11 March 2020. Moreover, the GSE41613 dataset was downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) database (Lohavanichbutr et al., 2013). The specific clinical information was listed in Supplementary Table S1. Before surgery, all patients did not receive radiotherapy or chemotherapy. Primary tumor site of all patients was the same. The “HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION” gene set containing 200 genes (Supplementary Table S2) was obtained from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).

In the TCGA dataset, 328 OSCC patients with complete clinical information were regarded as the discovery set. In the GSE41613 dataset, raw microarray data of 97 patients that possessed complete survival information were pre-corrected, transformed by log2 and normalized, followed by gene annotation. This dataset was utilized as the validation set.

Prognosis-related EMT genes were firstly screened in OSCC. Univariate cox analyses were employed to screen EMT genes with p < 0.05 in the TCGA dataset. By applying glmnet package, an optimal prognostic model was built with least absolute shrinkage and selection operator (LASSO) Cox regression analysis based on the prognosis-related EMT genes (Friedman et al., 2010). The optimal value of penalty parameter λ was determined via ten-fold cross-validation. The risk score of OSCC samples was computed according to regression coefficient as well as expression value of each gene in this model. The formula of the risk score was as follows: risk score = ∑ i = 1 n L A S S O   c o e f f i c i e n t   o f   g e n e   i * e x p r e s s i o n   o f   g e n e   i . OSCC patients in discovery and validation datasets were separated into high- and low-risk subgroups according to the median value of risk score, respectively. High-risk subgroup was defined as patients who had risk scores greater than the median while low-risk subgroup was defined as patients who had risks scores less than the median. To compare the difference in overall survival (OS) time between two subgroups, Kaplan–Meier curves were depicted by survival package and difference was determined with log-rank tests. Time-dependent receiver operating characteristic (ROC) curves of 1-year, 3-year and 5-year OS time were conducted for calculating the area under the curve (AUC) values to assess the predictive efficiency of the gene signature and other clinical features (age, gender, grade and stage) by applying survivalROC package (Heagerty and Zheng, 2005). Furthermore, univariate cox analyses were presented for evaluating the relationships of OSCC prognosis with the gene signature and clinical features. To validate whether the gene signature could be independently predictive of patients’ prognosis, multivariate cox analyses were performed based on prognosis-related factors with p < 0.05. Hazard ratio (HR) as well as 95% confidence interval (CI) was computed. Factors with HR > 1 were risk factors and those with HR < 1 were protective factors.

For determining whether the gene signature was independent of other clinical features, age, gender, grade, and stage were separated into high- and low-risk subgroups in the discovery dataset. OSCC subjects were stratified into age >65 and <65 subgroups, female and male subgroups, grade I-II/III-IV subgroups, stage I-II/III-IV subgroups. Kaplan-Meier OS analysis was carried out in each subgroup and difference was evaluated by log-rank test.

Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) package (https://sourceforge.net/projects/estimateproject) (Yoshihara et al., 2013) was employed for inferring stromal score, immune score, and tumor purity in each OSCC sample from the TCGA dataset based on gene expression profiles.

By applying CIBERSORT package (https://cibersort.stanford.edu/) (Newman et al., 2015), the infiltration levels of different immune cells were inferred in OSCC tissue specimens from the TCGA dataset according to gene expression profiling. Only specimens with p < 0.05 could be retained for further analyses. LM22 leukocyte gene signature matrix was used as the reference, which contained 547 genes that distinguished 22 human hematopoietic cells as follows: 7 kinds of T cell types, naïve and memory B cells, plasma cells, NK cells as well as myeloid subsets.

Gene Set Enrichment Analysis (GSEA) software was utilized for exploring the activated signaling pathways through comparing high- and low-risk subgroups from the TCGA dataset (Subramanian et al., 2005). The KEGG gene set (c2. cp.kegg.v7.1. symbols) was used as reference. 1,000 gene-set permutations were carried out. The terms with normalized enrichment score |NES| > 1.5 and FDR <0.05 were chosen as distinct pathways activated in high- or low-risk phenotypes, which were used for multiple GSEA gene sets.

Mutation Annotation Format (MAF) of OSCC samples was downloaded from TCGA database. These specimens were equally separated into high- and low-risk subgroups. The waterfall plots of two subgroups were depicted for illustrating the different mutation events using the Maftools package (Mayakonda et al., 2018).

Univariate cox analyses were applied to evaluate the prognostic value of each gene in the prognosis-related gene signature for OSCC patients from the TCGA dataset. Moreover, their expression was visualized in OSCC and normal tissue specimens. Spearman correlation between genes in this signature was evaluated in OSCC samples. The protein expression of genes from the prognostic gene signature in OSCC tissues was assessed through the Human Protein Atlas (https://www.proteinatlas.org/) (Colwill and Gräslund, 2011).

Totally, 40 paired OSCC and adjacent normal frozen tissues were collected from patients who experienced operation in the Taihe Hospital, Hubei University of Medicine between January 2020 and January 2021. All subjects did not receive chemoradiotherapy before operation. Diagnosis and staging were performed by experienced pathologists according to the American Joint Committee on Cancer Staging System. The study protocol gained the approval of the Ethics Committee of Taihe Hospital, Hubei University of Medicine (KY 2020-024), with written informed consent acquired from each subject. In addition, the research followed the Declaration of Helsinki.

Each gene in this prognostic model was verified in OSCC by RT-qPCR. Total RNA was isolated from OSCC tissues utilizing TRIzol reagent (Beyotime, China), which was reverse transcribed into cDNA utilizing primers and SuperScriptIII reverse transcriptase. RT-qPCR was carried out through Prime Script RT Reagent Kit and 7500 Real-Time PCR System (Applied Biosystems, United States). GAPDH served as the reference control. These amplification procedures included: denaturation at 95 °C lasting 5 min, 40 cycles of denaturation at 95 °C lasting 15 s, annealing at 55 °C lasting 30 s, and extension at 60 °C lasting 1 min. Table 1 listed the primer sequences of target genes. The relative expression was determined with 2^−ΔΔCT method.

From The Human Protein Atlas (https://www.proteinatlas.org/), immunohistochemistry staining of genes in the EMT gene signature in OSCC and normal oral tissues from 10 OSCC patients. Staining, intensity, quantity and location were also obtained. The used antibodies are as follows: AREG (HPA008720), COL5A3 (HPA048256), GAS1 (HPA066902), PLOD2 (CAB025898) and SFRP1 (CAB008116).

Statistical analysis was carried out by R packages (version 3.5.2) and GraphPad Prism software (version 8.0.1). Comparisons between two subgroups were presented via Student’s t-test or Wilcoxon rank-sum test. Kaplan-Meier survival curves were conducted, and survival difference was analyzed through log-rank test. The predictive efficacy was estimated with ROC curves. Spearman’s test was executed for correlation analysis. p < 0.05 was considered significant.

To screen prognosis-related EMT genes in OSCC, univariate cox regression analysis was employed in the TCGA dataset. As a result, 11 genes including ANPEP, AREG, COL5A3, DKK1, FMOD, GAS1, GPX7, PLOD2, SFRP1, TNFRSF11B, VEGFA were significantly associated with OSCC patients’ prognosis (Table 2). These prognostic genes were further assessed by LASSO analysis. Totally, 9 genes were screened for establishing the LASSO model (Figure 1A). The regression coefficient of each gene was calculated, as shown in Figure 1B. The risk score of each patient was then determined, as follows: AREG expression * 0.0471791751815675 + COL5A3 expression * (−0.0543433112313582) + DKK1 expression * 0.0619741505236927 + GAS1 expression * (−0.118881742831249) + GPX7 expression * (−0.142215922759831) + PLOD2 expression * 0.259497264609766 + SFRP1 expression * (−0.0579599175758466) + TNFRSF11B expression * (−0.197512366515929) + VEGFA expression * (0.053308700735239). We separated all patients in the discovery dataset into two subgroups according to the median value of risk score (Figure 1C). Survival status was further compared in the two subgroups. There were more patients with dead status for high-risk in comparison to low-risk subgroups (Figure 1D). Heat map depicted the different expressions of the genes (SFRP1, TNFRSF11B, PLOD2, GPX7, COL5A3, GAS1, VEGFA, AREG and DKK1) in this prognostic model between high- and low-risk subgroups (Figure 1E). Our further analysis demonstrated that high-risk patients exhibited worse survival time in comparison to low-risk subjects (p = 6.615e-05; Figure 1F). These data indicated that the risk score could be employed for predicting OSCC patients’ clinical outcomes. We further assessed the predictive efficacy of the risk score for OSCC prognosis by ROCs. The AUCs under 1-year, 3-year and 5-year OS were separately 0.669, 0.715 and 0.622, confirming the well predictive performance for clinical outcomes (Figure 1G). Also, we compared the predictive value of the risk score with other clinical features. Our data demonstrated that the risk score displayed the highest AUC value (0.715) for OS time among age (AUC = 0.575), gender (AUC = 0.487), grade (AUC = 0.557) and stage (AUC = 0.625; Figure 1H), indicating that this signature was more advantageous in comparison to other clinical features regarding prediction of survival time.

The GSE41613 cohort was employed to externally validate the predictive efficacy of this EMT gene signature in OSCC patients’ prognosis. With the same formula, we calculated the risk scores of OSCC patients. Consistently, we separated OSCC subjects into two subgroups according to the median value of risk score (Figure 2A). Low-risk patients exhibited more optimistic survival status than high-risk individuals (Figure 2B). The expression of the genes in this model was visualized in each OSCC sample (Figure 2C). The survival difference between subgroups was further validated in the GSE41613 dataset. As expected, high-risk patients displayed shorter OS time than low-risk patients (p = 7.869e-04; Figure 2D). ROCs were conducted for evaluation of the predictive efficacy of the risk score. In Figure 2E, AUCs under 1-, 3- and 5-year OS were separately 0.800, 0.778 and 0.729, confirming that this risk score could predict OS time of OSCC patients. Univariate analysis revealed that age, stage, risk score displayed distinct associations to OSCC patients’ prognosis in the TCGA dataset (Figure 2F). As confirmed by multivariate analysis, age, stage as well as risk score were independently related to prognosis (Figure 2G).

The predictive potency of this EMT prognostic model was further evaluated among different subgroups from the discovery cohort. Kaplan-Meier OS analysis showed that high-risk patients were predicted to have worse clinical outcomes compared with low-risk patients in age >65 (p = 0.012; Figure 3A), age <65 (p = 0.001; Figure 3B), female (p = 0.064; Figure 3C), male (p < 0.001; Figure 3D), grade I-II (p = 0.001; Figure 3E), grade III-IV (p = 0.010; Figure 3F), stage I-II (p = 0.112; Figure 3G) as well as stage III-IV (p = 0.002; Figure 3H) subgroups.

ESTIMATE algorithm was employed to assess stromal score, immune score, and tumor purity of OSCC samples from the TCGA dataset. Our data showed that higher stromal (p < 0.001) and immune scores (p < 0.001) were found in high-risk samples than low-risk samples (Figure 4A). Also, there was lower tumor purity in high-risk samples compared with low-risk samples (p < 0.001). The infiltration levels of 22 kinds of immune cells of OSCC specimens were determined by applying CIBERSORT algorithm. There were lowered infiltration levels of naïve B cells (p < 0.001), T follicular helper cells (p < 0.01), Tregs (p < 0.001), T gamma delta cells (p < 0.05) and resting mast cells (p < 0.001) in high-risk specimens compared with low-risk specimens (Figure 4B). Meanwhile, higher infiltration levels of CD4 memory activated T cells (p < 0.05), resting NK cells (p < 0.05), activated dendritic cells (p < 0.05), activated mast cells (p < 0.001) and eosinophils (p < 0.01) were examined in high-risk compared with low-risk subgroups.

GSEA was applied to explore signaling pathways associated with the EMT gene signature. As a result, basal transcription factors (NES = 2.04, NOM p = 0.002 and FDR = 0.008), base excision repair (NES = 1.82, NOM p = 0.006 and FDR = 0.047), cell cycle (NES = 1.96, NOM p = 0.006 and FDR = 0.017), nucleotide excision repair (NES = 2.08, NOM p = 0.002 and FDR = 0.006) and spliceosome (NES = 2.13, NOM p < 0.001 and FDR = 0.004) were distinctly upregulated in high-risk OSCC samples (Figure 5A). Moreover, calcium signaling pathway (NES = −2.21, NOM p < 0.001 and FDR <0.001), cytokine-cytokine receptor interaction (NES = −1.99, NOM p = 0.002 and FDR = 0.008), ECM receptor interaction (NES = −1.98, NOM p = 0.006 and FDR = 0.009), MAPK signaling pathway (NES = −1.94, NOM p < 0.001 and FDR = 0.011) and VEGF signaling pathway (NES = −2.18, NOM p < 0.001 and FDR <0.001) were activated in low-risk samples (Figure 5B). The somatic mutation was further assessed in high- and low-risk OSCC samples. Our data showed the first 20 mutated genes across OSCC samples. We found that higher frequent genetic mutations occurred in high-risk subgroup (Figure 5C) than low-risk subgroup (Figure 5D), especially TP53, FAT1, and CDKN2A.

Prognostic value of each gene in this EMT gene model was evaluated for OSCC samples from the TCGA dataset. Our univariate cox regression analysis demonstrated that high expression of AREG (p < 0.001, HR = 1.85, 95%CI: 1.31-2.61; Figure 6A), DKK1 (p < 0.001, HR = 2.05, 95%CI: 1.47-2.87; Figure 6B), PLOD2 (p = 0.01, HR = 1.58, 95%CI: 1.13-2.22; Figure 6C) and VEGFA (p = 0.019, HR = 1.51, 95%CI: 1.03-2.21; Figure 6D) was indicative of poorer prognosis of OSCC patients than their low expression. Furthermore, high expression of COL5A3 (p = 0.036, HR = 0.66, 95%CI: 0.46-0.94; Figure 6E), GAS1 (p = 0.006, HR = 0.61, 95%CI: 0.44-0.86; Figure 6F), GPX7 (p = 0.004, HR = 0.58, 95%CI: 0.41-0.82; Figure 6G), SFRP1 (p = 0.002, HR = 0.58, 95%CI: 0.42-0.82; Figure 6H) and TNFRSF11B (p < 0.001, HR = 0.55, 95%CI: 0.39-0.77; Figure 6I) displayed significant associations with prolonged survival time compared with their low expression. The prognostic implications of above genes were also confirmed in the GSE41613 cohort (Figures 7A–I).

The expression of genes in the EMT gene signature was compared between OSCC and normal tissues. Our data showed that AREG (Figure 8A), COL5A3 (Figure 8B), DKK1 (Figure 8C), GAS1 (Figure 8D), GPX7 (Figure 8E) and PLOD2 (Figure 8F) were significantly upregulated in OSCC than normal tissues (all p < 0.05). Furthermore, lower SFRP1 expression was found in OSCC compared to normal specimens (p < 0.05; Figure 8G). Our correlation analyses demonstrated that COL5A3 exhibited significant correlations to GAS1, GPX7, PLOD2, SFRP1 and TNFRSF11B in OSCC samples (Figure 8H). GAS1 exhibited significant correlations to GPX7, PLOD2, SFRP1 and TNFRSF11B. GPX7 was distinctly associated with PLOD2 and TNFRSF11B. These data indicated that there were distinct correlations between the genes in the EMT gene signature.

This study further confirmed gene expression in the EMT gene signature between 40 paired OSCC and normal specimens by RT-qPCR. Consistently, our data confirmed that AREG (Figure 9A), COL5A3 (Figure 9B), DKK1 (Figure 9C), GAS1 (Figure 9D), GPX7 (Figure 9E) and PLOD2 (Figure 9F) were distinctly highly expressed in OSCC compared with normal tissues (all p < 0.0001). Also, SFRP1 exhibited lower expression in OSCC than normal specimens (p < 0.0001; Figure 9G). Abnormal expression of AREG, COL5A3, GAS1, PLOD2 and SFRP1 was also confirmed in OSCC tissues by immunohistochemistry (Figure 9H). We also evaluated the difference in genes from the EMT gene signature across distinct pathological stages, as shown in Figures 10A–G. Among them, COL5A3, PLOD2 and SFRP1 were differentially expressed among pathological stages, indicative of their potential relationships with disease progression.