RNA sequencing (fragments per kilobase of transcript per million mapped reads [FPKM]) and clinical information, including age, sex, and TNM stage, were downloaded from publicly available data from TCGA head and neck squamous cell carcinomas and GEO databases. The primary site in these samples was the tongue. One TCGA dataset and three eligible GEO datasets (GSE65858, GSE41613, and GSE31056) provided data for background adjustment and quantified normalization.

The Kaplan–Meier survival package of R was used to perform Cox proportional hazard regression analysis to identify differentially expressed genes (DEGs) linked to overall survival (OS). The Limma R language package was applied to screen DEGs according to the criterion of false discovery rate (FDR) < 0.05 compared with the differences between normal and cancer tissues. Genes with prognostic value and similar cancer association were included in co-expression analysis.

Screened DEGs were prepared for further scale-free network construction and were inputted to test their availability based on the WGCNA R package (9). First, the appropriate soft threshold power β was determined via power calculation. Six clinical characteristics were incorporated: race, status, M grade, T grade, N grade, and stage. A sample clustering tree that was constructed combined clinical characteristics and used the dynamic tree cut algorithm to detect the modules. Finally, Pearson’s correlation coefficient was used to construct an adjacency matrix to describe the relationship between co-expression modules and clinical factors.

Identification of the co-expression module containing similar functional genes may reveal the major signaling pathways in TSCC. We also used the Metascape (http://metascape.org) portal to perform Gene Oncology (GO) analysis with representative enriched terms and to construct a network colored according to identified clusters. Furthermore, the clusterProfiler R package was used to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

Cox analysis was performed to calculate significant prognostic genes in TCGA. The genes were validated from the co-expression modules. LASSO Cox regression analyses were then applied to filter core genes and construct the prognostic models according to the risk score as follows:

Patients were assigned to high-risk (above the median cutoff of risk score) and low-risk groups (below the median cutoff of risk score). Kaplan–Meier survival and receiver operating characteristic (ROC) curve analyses were performed to further evaluate the prediction accuracy. To enhance prediction accuracy and interpretability, we downloaded three GEO cohorts (GSE65858, GSE41613, and GSE31056) as validation sets to select key prognostic modulators.

To further understand the correlations among these key prognostic modulators, consensus clustering analysis was performed to classify TSCC samples using the consensus cluster plus package. To explore the immune landscape of TSCC, the “estimate” component of R software was used to calculate tumor purity, stromal score, immune score, and ESTIMATE score of each tongue tumor sample. The CIBERSORT package was then used to accurately determine the composition of immune cells in the large tumor sample data. Significant results (p < 0.05) were used in further analysis.

The Cal-27 human TSCC cell line was obtained from Procell Company (Wuhan, China). Cal-27 cells were maintained in Dulbecco’s Modified Eagle Medium (DMEM; Procell Company) supplemented with 10% fetal bovine serum (FBS; Gibco, Carlsbad, CA, USA) and 1% penicillin-streptomycin liquid (Procell Company). Cal-27 cells were cultured in a humidified atmosphere containing 5% CO₂ at 37°C. Small interfering RNAs (siRNAs) targeting plasminogen activator urokinase (PLAU) and PPFIA binding protein 2 (PPFIBP2) were synthesized by GenePharma (Shanghai, China). Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) was used to transfect the siRNAs into Cal-27 cells according to the manufacturer’s protocol. The transfection efficiency was assessed by quantitative reverse transcription polymerase chain reaction (qRT-PCR).

TRIzol reagent (TaKaRa Bio, Dalian, China) was used to extract total RNAs following the manufacturer’s instructions. RNA (1 μg) was used to synthesize cDNA by a reverse transcription kit (RR037A; TaKaRa Bio). qRT-PCR was performed using the TB Green Premix Ex Taq Kit (RR820A; TaKaRa Bio). The 2^−ΔΔCt method was used to calculate the silencing efficiency of the siRNAs. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control. The primer sequences are summarized in Supplementary Table S1 .

The proliferation of Cal-27 cells was examined using the CCK-8 kit (Biosharp, Hefei, Anhui, China). Cal-27 cells were transfected with negative control or siRNAs targeting PLAU and PPFIBP2. CCK-8 (10 μl) was added to cells and incubated at 37°C for the indicated times. The optical density at 450 nm (OD₄₅₀) was measured using a microplate spectrophotometer.

The BeyoClick™ EdU-488 Cell Proliferation Kit (Beyotime, Shanghai, China) was used to stain cells to determine proliferation. Briefly, Cal-27 cells (1×10⁵ cells/well) were seeded in 6-well plates, transfected with different siRNAs, and cultured in an incubator. Cells were incubated with EdU for 2 h after 72 h of transfection, fixed with 1 ml of 4% paraformaldehyde for 20 min, and permeabilized with 0.3% Triton X-100 (Beyotime) for 20 min. Subsequently, the cells were incubated with 500 µl of the click reaction mixture for 30 min at room temperature in the dark, washed three times, and incubated with Hoechst stain for 10 min.

Cal-27 cells (4×10⁵ per well) in 6-well plates were transfected with different siRNAs. The cells were collected by trypsinization after 48 h of transfection. Cells (n = 1,000) were seeded in 6-well plates and incubated at 37°C for 2–3 weeks. Cell colonies were fixed with 4% paraformaldehyde for 30 min and stained with crystal violet solution (G1073; Solarbio, Beijing, China) for another 30 min. Fixed and stained cells that added up to 50 was counted as a clone.

Transwell chambers (3422; Costar, Corning, NY, USA) that were uncoated or coated with Matrigel Matrix (356234, Corning, USA) were used to detect the migration and invasion of Cal-27 cells. After 24 h of transfection, Cal-27 cells were collected by trypsinization. Cells (at a density of 5×10⁵/ml) were diluted in serum-free DMEM. The upper Transwell chamber was seeded with 200 μl of cell suspension (1×10⁵ cells/well). The lower chamber contained 600 μl of DMEM containing 20% FBS. The Cal-27 cells were incubated for 24 h. The cells that invaded across the chamber membrane were fixed and stained with crystal violet stain solution (G1073; Solarbio) for 30 min. Five randomly photographed fields were selected and the invading cells were counted using an inverted microscope.

Analyses were performed using R language software. Cox proportional hazards regression analysis was used to select the independent prognostic genes linked to OS. Kaplan–Meier curves were used to compare the clinical outcomes of the subgroups. The Wilcoxon rank-sum test was used to compare the ESTIMATE algorithm. In all analyses, statistical p-values were bilateral, and statistical significance was set at p < 0.05.

Cox analysis identified 1,391 DEGs as protective biomarkers and 1,346 DEGs as risk factors associated with poor OS ( Supplementary Table 1 ). Analysis using the Limma package revealed that 6,620 molecules prevented cancer from developing and 7,306 molecules might induce the formation of carcinoma ( Supplementary Table 2 ). We further identified 1,007 DEGs associated with TSCC. These were used to construct a weighted co-expression network ( Supplementary Table 3 ). To filter the DEGs, clinical traits, including status, race, and M, N, and T stages, were merged. The dendrogram and trait heatmap for these patients are shown in Figure 1A . The data show that the samples can be divided into multiple groups according to differences in clinical features.

The appropriate soft threshold power β was determined as six for subsequent adjacency calculations ( Figure 1B ). We identified gene co-expression modules based on the dynamic tree cut algorithm, as shown in Figure 1C . Finally, we evaluated and visualized the relationship between co-expression modules and clinical features, which indicated that the turquoise, green, and gray modules were significantly associated with the survival status ( Figure 1D ). Owing to the limited predictive value of TNM stage, we chose the turquoise module, which contained 244 genes, for further investigation. Applying Metascape, we conducted the GO function with representative enriched terms ( Figure 2A ) and a network colored by cluster ID ( Figure 2B ), which showed that the pathways of cell junction organization, response to topologically incorrect proteins, and regulation of cell adhesion may play an essential role in the turquoise module. KEGG analysis showed that focal adhesion pathways were enriched during the development of TSCC ( Figure 2C ).

To accurately filter crucial genes, data from the GSE65858 database were downloaded and prognostic biomarkers were identified using Cox analysis ( Supplementary Table 4 ). The results showed that 747 DEGs were correlated with better survival and 1,244 DEGs were associated with poor OS. The collective findings from TCGA, WGCNA, and GSE65858 identified 44 genes as crucial prognostic biomarkers. These were examined as the first step in establishing the risk model ( Supplementary Table 5 ). Based on the LASSO Cox regression analysis of the corresponding 44 mRNAs, 11 genes (PLXNB1, N4BP3, KDELR2, INTS8, PLAU, PPFIBP2, OAF, LMF1, IL34, ZFP3, and MAP7D3) were selected to build the risk model ( Figures 3A, B , Supplementary Table 6 ). Coefficient_mRNA of the risk score is shown in Supplementary Table 6 . Subsequently, patients were separated into low- or high-risk groups based on the median cutoff of risk score. The high-risk patients had worse OS than those in the low-risk group (p = 0.0002) ( Figures 3C, D ). Univariate and multivariate analyses showed that the risk score was an independent prognostic biomarker ( Figure 3E ). The ROC curve showed that the risk score had better prognostic effectiveness than other clinical traits, with an area under the ROC curve of 0.657 and 0.664 at 1 and 3 years, respectively ( Figure 3F ).

Pearson’s correlation analysis was performed to investigate the potential interactions involved in the selected genes and reveal relationships among these prognostic biomarkers ( Figure 4A ). Consensus clustering analysis assigned the TSCC samples into two clusters ( Figure 4B ). The expression pattern of genes in cluster 1 was remarkably different from that in cluster 2 ( Figure 4C ). Interestingly, Kaplan–Meier analysis revealed that the OS between the two clusters was significantly different (p = 0.0010; Figure 4D ). The findings suggest that TSCC is a heterogeneous cancer mediated by specific biomarkers.

TSCC is considered to be immunogenic, in part because of higher levels of tumor-infiltrating lymphocytes. Based on two clusters, a heatmap was drawn to comprehensively describe the immune landscape ( Figure 5A ). The ESTIMATE algorithm was then used to display immune-related scores. The stromal score in cluster 1 was higher than that in cluster 2 (p = 0.0002; Figure 5B ). Results of the CIBERSORT method showed that many immune cells contributed to the alteration of the immune environment, especially naïve B cells, plasma cells, CD8 T cells, CD4 memory resting and memory activated T cells, follicular helper T cells, and T regulatory cells (Tregs) ( Figure 5C ). These results indicated that functional clusters based on these prognostic biomarkers play an important role in deciphering the heterogeneous immune landscape of TSCC.

To further identify and verify the predictive ability of the 11 significant mRNAs, we downloaded one GEO (GSE65858) with OS and two GEO (GSE41613 and GSE31056) with progression-free survival (PFS). The risk model and two genes (PLAU and PPFIBP2) were successfully verified as effective prognostic biomarkers for OS and PFS ( Figures 6A–C , Supplementary Figures 1 - 3 ). N4BP3 presented a positive trend in all databases, yet the trend was opposite. Compared to those in normal tissues, the expression level of PLAU remarkably increased, and the expression level of PPFIBP2 noticeably declined. There was no significant change in the expression level of N4BP3 ( Figure 7A ). Therefore, we investigated the immune environment using these three biomarkers. The results showed that all three hub genes were intimately associated with the expression of PDCD1/PD-1, CD274/PD-L1, and multiple types of T cells ( Figures 7B, C ). Pan-cancer analysis explored differences of PLAU and PPFIBP2 in multiple cancers ( Supplementary Figures 6A and 7A ), as well as the correlation between immune cells and hub genes. The results indicated that PLAU and PPFIBP2 may influence the immune environment in multiple cancers ( Supplementary Figures 6B and 7B ).

We further investigated whether these mRNAs promoted the development and progression of TSCC through in vitro. We selected PLAU and PPFIBP2 after a careful literature review revealed that they have not been or rarely been studied in TSCC. The analysis sought to further verify the reliability and accuracy of our diagnostic model. Subsequently, siRNAs targeting PLAU and PPFIBP2 were designed. Compared with NC-transfected Cal-27 cells, the expression level of PLAU or PPFIBP2 was lower in cells transfected with si-PLAU-2 (named si-PLAU) or si-PPFIBP2-1 (named si-PPFIBP2) ( Figure 8A ). CCK8 and EdU assays showed that silencing PLAU significantly suppressed the proliferation of TSCC cells, whereas depletion of PPFIBP2 could promote the growth of TSCC cells ( Figures 8B, C ). Similarly, the proliferative capacities of TSCC cells were markedly suppressed after knockdown of PLAU and strikingly enhanced after depletion of PPFIBP2 by performing clone formation ( Figure 8D ), indicating that PLAU and PPFIBP2 play important roles in the growth. Finally, Transwell assay findings showed that the depletion of PLAU markedly attenuated the migration and invasion ability of TSCC cells, while knockdown of PPFIBP2 significantly enhanced the migration and invasion ability of TSCC cells ( Figure 8E ).