Head and neck cancer is the sixth most common types of cancer, with about half a million new cases annually diagnosed worldwide [1], among which 350,000 individuals die of it [2]. Although local control rate and quality of life have improved for head and neck squamous cell cancer (HNSCC) owing to advances in surgical techniques and comprehensive treatment techniques, overall survival (OS) has not increased significantly in recent decades. Moreover, the 5-years survival rate of patients with this disease is only 40–50% [3]. In recent years, age, clinical stage, and smoking status which are characteristics emerging as important contributors to the clinical outcome might help improve survival prediction for patients afflicted by HNSCC [4–6]. However, due to the complex molecular mechanisms underlying cancer regulation, the conventional clinical information has limited predictive ability.

Increasing evidence has revealed the important clinical significance of mRNA expression in various pathological and physiological processes of multiple tumor histotypes, including HNSCC. For instance, decreased calpain six expression has been known to be associated with tumorigenesis and poor prognosis of HNSCC [7]. Furthermore, FcGBP expression was upregulated by HPV infection and correlated with longer survival time of HNSCC patients [8]. However, compared with single-gene biomarkers, tumor signatures including several genes have been identified, which might be better choices to facilitate clinical application, provide insights into cancer progression, as well as reveal potentially new therapeutic targets [9, 10]. Thus, it is necessary to establish an expression-based gene signature for predicting survival of HNSCC patients for effective clinical decisions making with respect to optimal treatment regimen.

Metastasis is a complex, highly inefficient, but deadly process that has been under intense investigation in hopes to eliminate distant spread of metastasis and reduce cancer-associated mortality. A key event in promoting stationary tumor cells to migrate and invade is the epithelial mesenchymal transition (EMT) [11]. Specific tumor cell populations with EMT are associated with more aggressive tumor phenotypes and worse outcomes [12, 13]. Therefore, understanding the relationship between EMT and tumor is crucial for elucidating the underlying mechanism of tumorigenesis. Gene Set Enrichment Analysis (GSEA) is not concerned with a limited number of different genes that have significantly changed, but with whether the expression of these detected genes exhibit a common expression pattern in the defined functional groups, and interprets biological information to elaborate its biological significance from another perspective. In the present study, we tried to identify the relationship between the metastatic cascade system and HNSCC by GSEA analysis [14]. As expected, EMT-related risk signature could independently identify patients afflicted by HNSCC with high-risk and poor prognosis.

Recent studies have shown that clinicopathological features such as gender, age, tumor margin status, and metastasis are not sufficient for accurately predict the prognosis of patients. As a result, more and more mRNAs might be examined as molecular markers at an increasing rate to predict cancer development and prognosis, indicating that its clinical significance needs to be explored [18]. For example, analysis of the prognostic value of a single lncRNA from qRT-PCR array of 84 gastric cancer patient samples found that higher level of BANCR could predict a poor prognosis for GC patients [19]. The high expression of FAM83H-AS1 is involved in the progression of bladder cancer and serves as a prognostic biomarker and potential therapeutic target for patients with bladder cancer [20]. Compared to single biomarker, integrating multiple biomarkers into the aggregation model could improve the prognostic value [21]. In particular, the expression of a single gene could be controlled by a variety of factors, and hence might not provide a strong predictive effect. Therefore, a statistical model was constructed gene signature contains multiple genes, combined with the effect of each component gene prediction, improve forecasting efficiency. This model has been used widely, and is superior to single biomarkers in predicting disease prognosis [22–24].

With the development of high-throughput genetic testing technology, we are entering a new era of big biological data [25]. Currently, RNA sequencing or micro array data for gene mutations and expression levels often use Cox proportional hazard regression models to construct new prognostic features [26, 27]. In this study, we attempted to apply GSEA using expression data of 57,072 genes of patients from TCGA-HNSCC cohort, and we found significant differences in five biological functions with p-values < 0.05. As mentioned above, we focused on EMT, using GSEA to select genes to predict patient survival rather than exploring a broad range of genes. Univariate and multivariate Cox regression analyses were performed to identify combinations of four genes with predictive values for HNSCC patients, rather than just one gene [21]. Compared to some known predictive biomarkers, the risk profile of this option may have a more targeted and prognostic ability to support positive clinical results and be an effective classification tool for HNSCC patients. In addition, this study used bioinformatics methods to explore the risk characteristics of mRNA and its clinical significance, providing a novel method for mining potential prognostic markers, which not only supports our previous understanding of HNSCC, but also lays the foundation for future studies. It should be noted that, Kaplan-Meier analysis showed that risk score on survival is small and appears even smaller than lymph vascular invasion or perineural invasion. Therefore, more comprehensive and homogeneous datasets are needed to assess the four-mRNA signature before clinical applications. Notably, the stratification analysis found that the four-mRNA signature could generate superior performances for predicting the survival benefit in patients with P16 negative or alcohol or age <61, indicating that the four-mRNA signature could be helpful in predicting HNSCC prognosis, especially in P16 negative or alcohol or age <61 patients. It is highly significant given that P16 positive tumors typically have a good prognosis regardless of stage or high-risk pathologic features, this prognostic mRNA signature classifier for P16 negative-related HNSCC may help clinicians to pinpoint those HNSCC patients at high risk of unfavorable OS. EMT is recognized to be important in cancer cell migration/metastasis. For example, to metastasize, disseminated mesenchyme-like tumor cells of well-differentiated carcinomas must regain their epithelial function (invasion and metastasis) [28, 29]. The disruption of epithelial-cell homeostasis leading to aggressive cancer progression is correlated with the loss of epithelial characteristics and the acquisition of a migratory phenotype, known to EMT, and is considered to be a crucial event in malignancy [30, 31]. EMT and mesenchymal-related gene expression are associated with aggressive breast cancer subtypes and poor clinical outcome in breast cancer patients [32, 33]. Recently, some studies have reported the occurrence and development of gene expression induced EMT and metastasis in HNSCC [34, 35]. Additionally, several studies have been conducted which have predicted that in HNSCC patients, survival is associated with EMT. For instance, TEAD4 might act as a putative oncogenic gene by enhancing the proliferation, migration, and invasion of HNSCC cancer cells [36]. However, the EMT gene signature used to predict HNSCC prognosis have not been established. In the present study, we reported an EMT gene signature of four genes (PPIB, BASP1, WIPF1 and PLOD2), identified using bioinformatics methods which demonstrate the prognostic value of HNSCC. Of the four candidate genes, three have positively correlated coefficients in the prognostic model and correlated with poor survival. PPIB, called cyclophilin B, regulates the protein conformation of its substrate through prolyl cis–trans-isomerization in endoplasmic reticulum (ER) lumen and nucleus, and possesses multiple functions, including chemotaxis and prolactin signaling [37–39]. Additionally, it is a novel wild-type p53 (p53WT)-inducible gene [40]. In cancer biology, PPIB is associated with malignant progression and regulation of genes involved in the pathogenesis of gastric cancer [41],hepatocellular carcinoma [42], pancreatic cancer [43] and is considered as a candidate biomarker for these cancers. However, the exact molecular mechanism by which PPIB leads to cancer cell survival is unclear. BASP1 (brain acid soluble protein 1), was initially isolated from 23,000 brain extracts. BASP1 comprises an effector domain which dynamically couples to the plasma membrane and is also involved in neuronal sprouting process [44]. Recently, some reports have identified BASP1 as a transcriptional co-suppressor for WT1 protein [45]. WT1 can promote tumor cell proliferation, inhibition of apoptosis, and interacts with cytoskeletal proteins to promote migration, invasion and angiogenesis [46]. WT1 gene is expressed highly abnormally in a variety of tumors, including breast [47], thyroid [48, 49], non-small cell lung cancers [50], and HNSCC [51], and is considered to have the characteristics of an oncogene. PLOD2, also called LH2, is the key enzyme mediating the formation of stabilized collagen cross-links [52]. Gain and loss of function studies show that LH2 hydroxylated telopeptide lysine residues on collagen, shifted the tumor stroma toward higher levels of hydroxylysine aldehyde–derived collagen cross-links (HLCCs), lower levels of lysine aldehyde–derived cross-links (LCCs), increased tumor stiffness, and enhanced tumor cell invasion and metastasis [53, 54]. PLOD2 is overexpressed in different cancers such as bladder cancer [55], renal cell carcinoma [56], and oral carcinoma [57], and is closely associated with poor prognosis [58]. All these three genes are associated with EMT in cancer. In contrast, WIPF1 confers an onco-protective effect. WIPF1, also known as WIP, the WASP-interacting protein (WIP) drives the oncogenic activity of mutant p53. WIPF1 knockdown in glioblastoma and breast cancer cells expressing mtp53 greatly reduced the proliferation and growth capacity of cancer stem cell-like cells and reduced the expression of cancer stem cell-like markers (CD44, CD133 or TAZ/YAP). WIPF1 knockdown inhibits the growth of glioblastoma tumor cells and breast cancer cells in vivo [59]. Since EMT is affected in HNSCC, the targeted genes might be useful in controlling cancer. Finally, the prognostic value of the four genes in HNSCC still needs further experimental studies to elucidate biological function.

In conclusion, we identified and verified a four-gene risk signature related to EMT that can predict the survival of HNSCC patients, wherein a high-risk score shows poor patient prognosis, suggesting that this gene signature could be a potential new biomarker for HNSCC prognosis. Further investigation of these genes will provide theoretical guidance for basic research and for the clinical treatment of HNSCC.