Identification and Validation of an EMT-Related LncRNA Signature for HNSCC to Predict Survival and Immune Landscapes

Abstract

Long noncoding RNAs (lncRNAs) are increasingly recognized as decisive factors in the progression of head and neck squamous cell carcinoma (HNSCC), and they participate in the epithelial–mesenchymal transformation (EMT) of HNSCC. LncRNAs are closely related to the prognosis of patients with HNSCC; thus, it is essential to identify EMT-related lncRNAs with prognostic value for HNSCC. The coexpression network of EMT-related lncRNAs was constructed using The Cancer Genome Atlas (TCGA). An EMT-related eight-lncRNA-based prognostic signature was constructed using LASSO Cox regression and Cox proportional hazards analyses. Univariate and multivariate analyses and stratified prognosis confirmed that the prognostic signature was an independent predictive factor. Subsequently, we performed immune cell infiltration analysis, gene set enrichment analysis (GSEA), and single-sample GSEA (ssGSEA) pathway enrichment analysis to uncover the potential molecular mechanisms of prognostic differences in the high- and low-risk groups. Next, we discussed the relationship between the prognostic signature and immune checkpoint-related genes, their TIDE scores, and the sensitivity of common chemotherapeutics. Finally, we further verified the expression differences in lncRNAs that were included in our signature via RT–qPCR in eighteen paired tissues. In summary, this prognostic signature provides powerful prognostic biomarkers for HNSCC and could serve as a predictor for the sensitivity of common chemotherapeutics and immunotherapy responses as well as providing a reference for further personalized treatment.

Introduction

Human head and neck squamous cell carcinoma (HNSCC) is one of the most common malignant tumors and the most crucial pathological type of head and neck cancer, with an incidence of approximately 600,000 new cases yearly (Bray et al., 2018). The prognosis of HNSCC is poor due to high proliferation, regional lymph node metastasis, and a high recurrence rate; the 5 years survival rate of HNSCC is nearly 50% (Burusapat et al., 2015; Lambert et al., 2017). The prognosis of HNSCC patients is still unsatisfactory, although treatment methods have been improved in recent decades (Yi et al., 2020).

The first step in metastasis of cancer cells is achieving local invasiveness. To that end, cancer cells must abandon epithelial phenotypes and acquire mesenchymal morphological and transcriptional characteristics. This process is usually called the epithelial–mesenchymal transformation (EMT). EMT is crucial for physiological and pathological processes, such as embryonic development, wound healing, and metastasis of malignant tumors (Lamouille et al., 2014; Nieto et al., 2016). During EMT, there is separation between epithelial cells and between epithelial cells and the basement membrane (Nieto et al., 2016). EMT is essential for metastasis and the acquisition of drug resistance. Nearly 90% of cancer patients die from metastasis rather than primary tumors. Therefore, EMT has an important effect on the prognosis of HNSCC patients. Nevertheless, at present, the mechanism of EMT has not been fully elucidated.

Long noncoding RNAs (lncRNAs) are defined as transcripts of more than 200 nucleotides that are not translated into proteins (Xing et al., 2019; Guo et al., 2020). LncRNAs are crucial for cellular processes. An increasing number of studies have shown that lncRNAs have a regulatory effect on EMT (Lu et al., 2017; Jia et al., 2019). The function of lncRNAs in regulating cell fate has received widespread attention. The role of lncRNAs has been reported in multiple cancers, such as lung cancer, gastric cancer, bladder cancer, and colorectal cancer (Lv et al., 2017; Zhang et al., 2018; Wang et al., 2019a; Bure et al., 2019). Nevertheless, the mechanism of lncRNAs and EMT and the role of lncRNAs in the prognosis of HNSCC have still not been elucidated.

EMT is indispensable for immunosuppression, immune tolerance, evasion, and metastasis (Dongre et al., 2017; Mittal, 2018). EMT influences the prognosis of patients, simultaneously providing an approach to discovering novel therapeutic targets and predictive biomarkers for tumor treatment (Jiang and Zhan, 2020). Immune cells occupy a large proportion of the tumor microenvironment of HNSCC. Adaptive immune resistance induces an immunosuppressive tumor environment that enables immune evasion (Srinivasan et al., 2018). This capability likewise enhances tumor development and metastasis. Immunotherapy, a neoteric treatment strategy for tumors, has been widely studied in recent years, and the results pertaining to immune checkpoint inhibitors (ICIs) are particularly satisfactory. Tumor cell-derived PD-L1 contributes to EMT and tumor invasion in multiple tumor types (Dong et al., 2018). Anti-PD-1 and anti-CTLA4 strategies occupy the forefront of the immunotherapy field. Tumor stroma composition and immune cell infiltration have a profound impact on the prognosis of patients (Galon et al., 2006; Dongre et al., 2017). Reduced stromal composition indicates an increase in cancer cells, and less immune cell infiltration indicates a reduced response to immunotherapy, which would predict a worse prognosis.

The treatment of HNSCC is mainly surgery combined with radiotherapy and chemotherapy. Clinical doctors generally predict patient prognosis based on the TNM stage or clinical stage; however, this method is not sufficiently accurate. A novel prediction method urgently needs to be established for the clinic. Therefore, in this study, we established a new method to diagnose early disease and predict prognosis accurately. To this end, we further explored the regulatory mechanism of EMT in HNSCC and found potential effective biomarkers for early diagnosis and accurate prognosis. HNSCC RNA sequencing (RNA-seq) data and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database, and EMT-related genes were obtained from the study of Choi et al. (Choi et al., 2010; Yang et al., 2018). Univariate Cox regression analyses were used to identify key EMT-related lncRNAs. These lncRNAs were used to construct a prognostic signature with LASSO regression analysis (Vrieze, 2012). We divided the patients into high-risk and low-risk groups on the basis of our signature (Jiang et al., 2020). To evaluate the prognostic accuracy of this signature, we plotted Kaplan–Meier survival curves and receiver operating characteristic (ROC) curves and performed principal component analysis (PCA). We built a coexpression network to predict the function and mechanism of these lncRNAs. Univariate and multivariate Cox regressions were employed to appraise the independence of the signature. MultiROC curves were used to assess accuracy among multiple factors. A nomogram was used to predict the survival rate via multiple independent factors. Then, we compared differences between the two groups in terms of immune infiltration, sensitivity to frequently used chemotherapeutic drugs, and gene set enrichment. Finally, we verified lncRNA expression differences among eight paired paracarcinoma tissues and cancer tissues.

Materials and Methods

Data Collection and Preparation

RNA-seq data and clinical data were downloaded from the TCGA database (https://portal.gdc.cancer.gov/) (Yang et al., 2018). Fragments per kilobase million (FPKM) values were employed in the next analysis. The samples in this study, mostly from the tongue, oral floor, and cheek, included 44 normal structures and 525 malignant tumors from TCGA. The RNA-seq data and clinical data were combined by using R (R version 4.0.4). Protein-coding genes and lncRNAs were annotated and classified using the Ensemble human genome browser GRCh38.p13 (http://asia.ensembl.org/index.html) (Sun et al., 2020). This study excluded repeated samples and samples for which clinical data were not known.

Identification of Prognostic EMT-Related Hub LncRNAs

We used TCGA data to identify lncRNAs associated with EMT-related genes via the R software package “limma.” EMT genes were obtained from the study of Choi et al. (Choi et al., 2010). Pearson’s correlational coefficients assessed the association of lncRNAs and EMT-related genes. The threshold value of Pearson’s correlational coefficient was >0.3 or < −0.3, and the p value was <0.01, in the identification of the differentially expressed EMT-related lncRNAs of tumor tissue and normal tissue. The association of EMT-related lncRNA expression and prognosis was assessed by univariate Cox regression, and the overall survival rate was calculated using the R software package “survival.” A p value <0.001 indicated a significant difference.

Estimation of the Prognosis Signature

The data from TCGA were divided randomly and equally into two groups: a training set (n = 251) and a validation set (n = 248). We adopted LASSO regression via the R software package “glmnet” to filter prognosis-related lncRNAs. ROC curves were applied to verify this signature (Zhang et al., 2018). Then, the validation set further validated the signature. Risk scores of every patient from total TCGA were calculated (Yoshihara et al., 2013). Then, all HNSCC patients were divided into a high-risk group (n = 254) and a low-risk group (n = 245) based on the median risk score. Kaplan–Meier survival curves were used for eight EMT-related lncRNAs to explore the relationship between prognosis and lncRNAs. PCA was used to evaluate the expression profile between the two groups.

Validation of the Independence and Forecast Efficiency of the Prognostic Signature

Univariate and multivariate Cox analyses were utilized to discover independent risk factors among age, sex, grade stage, clinical stage, TNM stage, and risk score. MultiROC curves were employed to compare efficiency among selected risk factors.

CIBERSORT Analysis of Immune Condition and Function

We employed the CIBERSORT database to predict the differences among 22 immune cells in infiltration conditions of two groups (https://cibersort.stanford.edu/) (Newman et al., 2015). The algorithm of 100 permutations was adopted. Only samples with a p value <0.05 were included to perform the subsequent analysis comparing differential immune infiltration levels between two groups grouped by clustering subtypes and risk scores.

ssGSEA Assessment of Tumor-Infiltrating Cells

Bindea et al. found numerous marker genes of tumor-infiltrating immune cells (TIICs) (Angelova et al., 2018). Using ssGSEA, the enrichment scores of 16 immune cells and 13 immune functions for each HNSCC sample were quantified based on gene signatures using the R software package “gene set variation analysis” (GSVA) (Hänzelmann et al., 2013). The immune infiltration levels and functions between the two groups were compared by the Kruskal–Wallis test. ESTIMATE was used to evaluate the tumor microenvironment and predict the tumor purity and abundances of tumoral stromal cells based on the gene expression data of the high-risk group and low-risk group (Jiang et al., 2020).

Evaluation of Immunotherapy and Drug Sensitivity Prediction

We further assessed the response to immune checkpoint inhibitors (ICIs), such as CTLA4 and PD-1. The response to ICI was predicted by the Tumor Immune Dysfunction and Exclusion (TIDE; http://tide.dfci.harvard.edu/). The differences between the two groups were discovered by the Wilcoxon test. Furthermore, the therapeutic effects of anti-PD1 and anti-CTLA4 were predicted for HNSCC patients in each of the two groups. The sensitivity of each patient to commonly used chemotherapy drugs for HNSCC was estimated using the Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) database (Yang et al., 2013). The half-maximal inhibitory concentration (IC50) was quantified through the R software package “pRRophetic” (Geeleher et al., 2014).

Gene Set Enrichment Analysis

In terms of the signaling pathway analysis, differential expression analysis was first conducted on all genes of the samples with high- and low-risk scores by applying the R software package “limma.” Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses with GSEA software were performed on the differentially expressed genes (DEGs) between the high-risk and low-risk groups to research the signaling pathways in which DEGs are involved. (p < 0.05 and false discovery rate, FDR <0.25).

Deciphering Mutational Landscape in the Genome

Mutation data of patients with HNSCC were downloaded from the TCGA database. The top 20 driver genes with the highest frequency of change were further analyzed. Oncoplots were sketched between high- and low-risk groups by R software “maftools.”

Verification of Risk Signature Stability and Prediction of Survival Rate

Kaplan–Meier survival curves were employed to analyze the stability of the signature among imparity clinical features, including age and sex. We ultimately structured a nomogram to predict the survival rate accounting for the independent features.

Verification of Expression Differences

LncRNAs incorporated into signature expression differences of eighteen paired paracarcinoma tissues and cancer tissues were verified by qRT–PCR. The primers used for qRT–PCR are summarized in Table 1. Clinical data of all patients recruited were showed in Supplementary Table S1. According to the manufacturer's instructions, total RNA was extracted and reverse transcribed from tissue from eighteen patients by a Takara Reagent kit. The SYBR Master Mixture Kit (Takara) was utilized to amplify the resulting cDNA, and qRT–PCR was performed using a QuantStudio(TM) 6 Flex System according to the manufacturer's instructions. Each experiment was conducted at least three times. The expression of target lncRNAs was calculated using the 2-∆Ct method relative to the internal reference (ACTB).

Statistical Analysis

All data were analyzed by R software. A p value <0.05 was considered a significant difference. The independent Student’s t test for continuous data and the χ² test for categorical data were utilized for pairwise comparisons between groups. The Wilcoxon test was used to compare non-normally distributed variables between two groups. The log-rank test was used to compare two groups when Kaplan–Meier survival curves were generated.

Results

Identification of EMT-Related LncRNAs

A total of 14,142 lncRNAs were identified by analyzing RNA-seq data from HNSCC patient tissue samples in the TCGA database and 200 EMT-related genes. A total of 1,598 lncRNAs were screened by calculating Pearson’s correlation coefficients between the lncRNAs and EMT-related genes. The identified lncRNAs were defined as EMT-related lncRNAs. We conducted our study as shown in Figure 1.

FIGURE 1

Establishment of an EMT-Related Eight-lncRNA Prognostic Signature

Out of 1,598 EMT-related lncRNAs, 22 lncRNAs were selected through univariate Cox regression analysis and were associated with the prognosis of HNSCC patients (Figures 2A,B), and the p value was <0.001. Patients were randomly and equally divided into training and validation sets. LASSO regression analysis was conducted to assess the correlations between these 22 lncRNAs and survival status in the training set (Figures 2C,D). Then, we included eight lncRNAs as a predictive prognostic signature. Of the 8 lncRNAs, 6 lncRNAs were protective factors (hazard ratio, HR < 1), and 2 lncRNAs were risk factors (HR > 1) (Figure 2E). The network of 8 EMT-related lncRNAs and EMT-related genes was revealed by a network diagram (Figure 2F). Only AL138902.1 negatively regulated the downstream genes LGALS1, NNMT, and TNFRSF12A in 18 pairs of interactions.

FIGURE 2

Evaluation of the Prognostic Signature

Risk scores of each HNSCC patient were calculated, and the patients were divided into high-risk and low-risk groups based on the median risk score in the training set. A correlation diagram between the survival rate and risk score (Figure 3A), a scatterplot of the patients (Figure 3B), a heatmap (Figure 3C), Kaplan–Meier survival curves (Figure 3D), and 1, 3, and 5 years survival ROC curves (Figure 3E) were used to evaluate this signature. The validation set and all sets were used to further verify the accuracy of this signature (Figures 3F–O). HNSCC patients were sorted according to their risk scores, and the survival rate was significantly associated with the risk score. The survival time of patients in the high-risk group was lower than that of patients in the low-risk group. This result suggested that the higher the score was, the worse the prognosis. In the heatmap, the expression of 8 EMT-related lncRNAs was remarkably different between groups. The expression of AL138902.1, AL356481.3, AC024267.3, LINC00996, AC114730.3, and AP003774.4, as risk factors, was higher in the high-risk group than in the low-risk group. On the other hand, the protective factors AC245041.2 and AL596223.1 were expressed in the opposite direction. According to the Kaplan–Meier survival curves, patients with low-risk scores survived notably longer than those with high-risk scores. ROC curves were plotted to assess the predictive value of the prognostic signature. The training set AUC values of the 1, 3, and 5 years ROC curves were 0.723, 0.735, and 0.683, respectively, and the AUC values of the 1, 3, and 5 years ROC curves of all patients were 0.688, 0.705, and 0.610. These results confirmed that our signature had medium accuracy for evaluating the prognosis of HNSCC patients. PCA was conducted to investigate the distributional differences based on different patterns (Figures 3P–S). We found that eight EMT-related lncRNAs were better prognostic genes by which high- and low-risk patients could be more easily separated than risk genes, EMT genes, and all risk genes.

FIGURE 3

Univariate and multivariate Cox regression analyses demonstrated that the risk score was an independent predictor of prognosis in HNSCC patients, as were age, M stage, and N stage (Figures 4A,B). Furthermore, the efficacy of the risk score was higher than that of other clinical factors at 1 year (Figure 4C), 3 years (Figure 4D), and 5 years (Figure 4E) ROC curves.

FIGURE 4

Tumor Immune Microenvironment Characteristics of Different Groups

We estimated the tumor immune microenvironment using the CIBERSORT database, ssGSEA, and ESTIMATE. As shown in Figures 5A,B, we mapped the immune landscape for each patient. Sixteen immune cells were remarkably different between the two groups, as shown in the violin plot (Figure 5C). Of the 16 immune cells, 11 kinds of cells were reduced, and five kinds were increased. The infiltration of naïve B cells, plasma cells, CD8 T cells, activated memory CD4 T cells, follicular helper T cells, Tregs, gamma delta T cells, resting NK cells, activated NK cells, and M2 macrophages was lower in the high-risk group. However, CD4 naïve T cells, M0 macrophages, resting dendritic cells, activated dendritic cells, activated mast cells, and neutrophils were more abundant in the high-risk group. We used the CIBERSORT database to explore correlations between immune cells in all HNSCC patients (Figure 5D). T cell memory activation had a high correlation with resting NK cells (correlation coefficient, CC = 0.44) and CD8 T cells (CC = 0.65). Naïve B cells had a high correlation with plasma cells (CC = 0.57). Resting memory T cells had a high correlation with CD8 T cells (CC = −0.44). B cell memory was correlated with CD4-naïve T cells (CC = 0.72). Moreover, we evaluated the correlation between the risk score and immune cells (Table 2) and discovered that naive B cells, CD8 T cells, CD4 memory activated T cells, and resting NK cells had a negative correlation, and activated NK cells, resting dendritic cells, activated dendritic cells, activated mast cells, and neutrophils had a positive correlation.

FIGURE 5

In the ssGSEA, we found that 15 of 16 tumor-infiltrating immune cells (except iDCs) (Figure 5E) and 10 of 13 immune functions (except APC costimulation, parainflammation and type I IFN response) (Figure 5F) were clearly reduced in the high-risk group compared with the low-risk group.

In the integration of CIBERSORT and ssGSEA, most immune infiltrating cells were reduced, such as B cells, neutrophils, and CD8⁺ T cells. This result implied that as the tumor progressed, fewer responsive immune cells infiltrated the tumor, which might lead to the insensitivity of advanced tumors to immunotherapy.

The risk index was conspicuously negatively correlated with the immune, stromal, and ESTIMATE scores, indicating that the infiltration levels of immune and stromal cells decreased with the elevation of the risk score of HNSCC. Higher immune or stromal scores in tumors indicated higher abundances of immune cells or stromal cells. We discovered that patients in the high-risk group had lower stromal, immune, and ESTIMATE scores. Moreover, the higher the score, the better the immune or stromal condition.

Assessment of Immunotherapy and Drug Sensitivity

We discovered that 37 kinds of immune checkpoint genes were notably different (Figure 6A). CD44, CD276, and TNFSF9 were more highly expressed in the high-risk group than in the low-risk group. The expression of 34 other types of genes was the opposite. This result indicated that treatment targets in the high-risk group should differ from those in the low-risk group. We compared anti-PD-1 and anti-CTLA4 therapeutic differences in the two groups via TIDE data (Figure 6B) and found that patients in the low-risk group had better responses.

FIGURE 6

We compared the differences in the estimated IC50 levels of five chemotherapy drugs, including paclitaxel (Figure 6C), docetaxel (Figure 6D), AKT inhibitor VIII (Figure 6E), docetaxel (Figure 6F), and cisplatin (Figure 6G). Our data revealed higher estimated IC50 levels of paclitaxel (p = 0.00029) and docetaxel (p = 1.7e-9) in the low-risk group than in the high-risk group. In contrast, the estimated IC50 levels of the AKT inhibitor VIII (p = 1.1e-5), docetaxel (p < 2.22e-16), and cisplatin (p = 3.4e-5) in the low-risk score group were significantly lower than those in the high-risk score group. This result implied that patients with high-risk scores were more sensitive to paclitaxel and docetaxel; patients with low-risk scores were more sensitive to the AKT inhibitor VIII, docetaxel, and cisplatin.

GSEA

GSEA was conducted to compare the high- and low-risk groups. Most of the KEGG enrichment was concentrated in “cell adhesion” and immune-related pathways, such as “chemokine signaling pathway,” “natural killer cell mediated cytotoxicity,” and “T cell receptor signaling pathway,” in the low-risk group (Figure 7A). Moreover, hallmark enrichment was associated with “epithelial-mesenchymal transformation” and “hypoxia” in the high-risk group and “G2M checkpoint,” “IL2 STAT5 signaling,” “KRAS signaling dn,” and “peroxisome” in the low-risk group (Figure 7B).

FIGURE 7

The Correlation Between the Prognostic Signature and Somatic Variants

We evaluated the distribution of somatic variants in patients with HNSCC driver genes between the low- and high-risk groups and applied the R package “maftools” to visualize the landscape of mutation profiles in patients with HNSCC(Figure 8). Among all types of mutations, missense mutations accounted for the highest proportion. Our research shows that the alteration frequency of TP53, CDKN2A, and NSD1 was significantly different between the low-risk group and the high-risk group (Table 3).

FIGURE 8

Prediction Stability and Prediction of Survival Rate

We found that our risk signature had high stability among different clinical subgroups, such as age, sex, grade, T stage, and N stage, via Kaplan–Meier survival curves (Figures 9A–E). Demographic data from TCGA patients were presented in Table 4. We established a nomogram including age, M stage, N stage, and risk score to predict 1, 3, and 5 years survival rates (Figure 6F). We found that the nomogram corresponded fairly with the actual survival rate.

FIGURE 9

Verification of Expression Differences

Using RT–qPCR, we found that eight lncRNAs included in our signature had distinctive expression differences between eighteen paired paracarcinoma tissues and cancer tissues. AC024267.3 (p = 0.021) (Figure 10A), AL596223.1 (p = 0.018) (Figure 10C), LINC00996 (p = 0.030) (Figure 10D), AL356481.3 (p = 0.037) (Figure 10F), and AL138902.1 (p = 0.029) (Figure 10G) were upregulated in tumors compared with normal tissues. In addition, AC114730.3 (p = 0.048) (Figure 10B), AP003774.4 (p = 0.018) (Figure 10E), and AC245041.2 (p = 0.047) (Figure 10H) were downregulated.

FIGURE 10

Discussion

EMT is a process in which epithelial cells acquire mesenchymal cell characteristics and is closely related to tumor invasion, metastasis, and chemotherapy resistance in tumor progression (Lamouille et al., 2014; Pastushenko and Blanpain, 2019; Xing et al., 2019; Guo et al., 2020). The molecular mechanisms involved in the transition from an epithelial phenotype to a mesenchymal state are complex and controlled by multiple signaling pathways (Gugnoni and Ciarrocchi, 2019). Recently, more and more evidence has emphasized the regulatory role of lncRNA in the tumor EMT process (Grelet et al., 2017; Chen and Shen, 2020; Hu et al., 2020; Peng et al., 2020).

In this paper, EMT-related lncRNAs were comprehensively identified in HNSCC for the first time, and an EMT-related lncRNAs prognostic signature was established and validated. Then, we evaluated the relationships of the prognostic signature with tumor immune microenvironment characteristics, somatic variants, and chemotherapy. Finally, we verified lncRNA expression differences among eighteen paired paracarcinoma tissues and cancer tissues.

The Pearson's correlational coefficients were used to identify EMT-related lncRNAs. Univariate Cox regression was used to search for EMT-related lncRNAs with prognostic values. LASSO regression analysis was used to construct a predictive model to identify patients as a high-risk and low-risk group (Figure 2). The AUC values of training set and test set confirm that our signature had medium accuracy for evaluating the prognosis of HNSCC patients (Figures 3E–O). Univariate and multivariate Cox regression analyses demonstrated that the risk score was an independent predictor of prognosis in HNSCC patients (Figures 4A,B). Additionally, PCA showed that the high-risk and low-risk groups were distributed in different directions, suggesting significant differences in immune status between the high-risk and low-risk groups (Figures 3P–S).

ESTIMATE, a tool based on gene expression signatures, was used to predict the proportions of immune cells and stromal cells in tumors (Yoshihara et al., 2013). The high-risk score group had higher tumor purity and lower stromal and immune scores than the low-risk group. Moreover, we found that patients with high stromal and immune scores had shorter survival times, indicating that high-risk scores could predict a poor prognosis for HNSCC patients. This finding ultimately indicated that EMT could be related to the reconstruction of the tumor microenvironment while affecting tumor growth and progression (Li et al., 2020; Xu et al., 2021a; Guo et al., 2021; Lin et al., 2021). EMT is a critical event in the process of tumor metastasis, and our results support this conclusion. There were also differences between the two groups in terms of the tumor immune microenvironment. The presence of inhibitory immune cells was higher in the high-risk group, and the proportion of effector cells was higher in the low-risk group, suggesting that the outcome of immunotherapy might be better in the low-risk group. Many factors impact the effect of immunotherapy in HNSCC. Despite the favorable clinical benefits of immunotherapy, only a portion of HNSCC patients show a response to immunotherapy. Using our signature, we found that the sensitivity of PD-1 and CTLA4 was clearly different in the high- and low-risk groups. Patients in the low-risk group responded to anti-PD-1 and anti-CTLA4 drugs. It is worth noting that TIDE is built based on data from melanoma and NSCLC cancers. Although it has been proved that TIDE score can be used to predict HNSCC in a lot of literature, its accuracy is still worth considering. We found that patients with high risk scores were more sensitive to paclitaxel and docetaxel. Meanwhile, patients with low risk scores were more sensitive to the AKT inhibitor VIII, docetaxel, and cisplatin. Based on these findings, clinicians can implement comprehensive and individualized treatment plans for patients.

The GSEA results showed a significant connection between the EMT-related LncRNA signature and the tumor immune microenvironment, including chemokines, natural cells, and T cells. Hallmarks of EMT and the “hypoxia” pathway enriched in the high-risk group promote tumor development and metastasis, and “G2M checkpoint”-related gene overexpression and deficiency of KRAS inhibition-related genes promote abnormal cell division in different tumors (Mittal, 2018; Shirazi et al., 2020; Smith et al., 2020). The peroxisome generally exhibits dysfunction in tumors, especially advanced tumors (Islinger et al., 2018).

Growing evidence indicates that TP53 mutations are potential predictors to immunotherapy (Caponio et al., 2020; Zhang et al., 2020). Meanwhile, the overall survival of TP53 mutant-type (TP53 MT) patients is worse than that of TP53 wild-type (TP53 WT) patients (Wu et al., 2020a). Our research shows that the alteration frequency of TP53, CDKN2A, and NSD1 was significantly different between the low-risk group and the high-risk group (Figure 8; Table 3).

Some of these genes have been reported in a variety of cancers. Our results are consistent with the conclusions of the existing articles. AC114730.3 is a protective factor for patients with HNSCC and bladder cancer. The higher the AC114730.3 expression, the better the prognosis (Wu et al., 2020b; Xu et al., 2021b). AC245041.2 could influence transcription factors via LAMA3 and then alter the expression of other protein molecules. Furthermore, the high expression of AC245041.2 supports the survival of patients with clear cell renal cell carcinoma (Tian et al., 2021). Wang et al. also similarly reported that AC245041.2 was a protective factor in clear cell renal cell carcinoma (Wang et al., 2019b). Lina et al. reported that LINC00996 was mainly involved in biological processes such as vascular development, positive metabolic processes, and growth regulation (Lina, 2021). Ge H et al. reported that low expression of LINC00996 was related to metastasis and a low survival rate. In this process, the JAK-STAT, NF-κB, HIF-1, TLR and PI3K-AKT signaling pathways might have a crucial effect (Ge et al., 2018).

Our study lays the foundation for further research into EMT-related lncRNAs and has profound implications for screening criteria for patients receiving immunotherapy. The original EMT-based classification could provide prognostic predictors for HNSCC and might guide clinicians in selecting potential patients who respond to immunotherapy for the preferential use of immunotherapy. This study still has a few limitations. First, our model may have some deviation due to the use of TCGA data and not external data. Although we only used TCGA data, univariate and multivariate Cox analyses, grouping modeling, stratified prognosis analyses, and ROC curves all indicated that our model was stable and accurate. Second, further prospective external verification with a large sample of HNSCC patients from our hospital is required. Future functional studies are essential to clarify the underlying mechanism of EMT-related lncRNAs in HNSCC.

Conclusion

In summary, we constructed and verified the independent predictive value of an EMT-related lncRNA-based prognostic signature for patients with HNSCC for the first time. More importantly, these results provide important references for personalized chemotherapy and immunotherapy.

References

• 1

  AngelovaM.MlecnikB.VasaturoA.BindeaG.FredriksenT.LafontaineL.et al (2018). Evolution of Metastases in Space and Time under Immune Selection. Cell175, 751–e16. 10.1016/j.cell.2018.09.018

  • CrossRef
  • Google Scholar

• 2

  BrayF.FerlayJ.SoerjomataramI.SiegelR. L.TorreL. A.JemalA. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinicians68, 394–424. 10.3322/caac.21492

  • CrossRef
  • Google Scholar

• 3

  BureI. V.NemtsovaM. V.ZaletaevD. V. (2019). Roles of E-Cadherin and Noncoding RNAs in the Epithelial-Mesenchymal Transition and Progression in Gastric Cancer. Int. J. Mol. Sci.20. 10.3390/ijms20122870

  • CrossRef
  • Google Scholar

• 4

  BurusapatC.JarungroongruangchaiW.CharoenpitakchaiM. (2015). Prognostic Factors of Cervical Node Status in Head and Neck Squamous Cell Carcinoma. World J. Surg. Onc13, 51. 10.1186/s12957-015-0460-6

  • CrossRef
  • Google Scholar

• 5

  CaponioV. C. A.TroianoG.AdipietroI.ZhurakivskaK.ArenaC.MangieriD.et al (2020). Computational Analysis of TP53 Mutational Landscape Unveils Key Prognostic Signatures and Distinct Pathobiological Pathways in Head and Neck Squamous Cell Cancer. Br. J. Cancer123, 1302–1314. 10.1038/s41416-020-0984-6

  • CrossRef
  • Google Scholar

• 6

  ChenS.ShenX. (2020). Long Noncoding RNAs: Functions and Mechanisms in colon Cancer. Mol. Cancer19, 167. 10.1186/s12943-020-01287-2

  • CrossRef
  • Google Scholar

• 7

  ChoiY.-L.BocanegraM.KwonM. J.ShinY. K.NamS. J.YangJ.-H.et al (2010). LYN Is a Mediator of Epithelial-Mesenchymal Transition and a Target of Dasatinib in Breast Cancer. Cancer Res.70, 2296–2306. 10.1158/0008-5472.can-09-3141

  • CrossRef
  • Google Scholar

• 8

  DongP.XiongY.YueJ.HanleyS. J. B.WatariH. (2018). Tumor-Intrinsic PD-L1 Signaling in Cancer Initiation, Development and Treatment: Beyond Immune Evasion. Front. Oncol.8, 386. 10.3389/fonc.2018.00386

  • CrossRef
  • Google Scholar

• 9

  DongreA.RashidianM.ReinhardtF.BagnatoA.KeckesovaZ.PloeghH. L.et al (2017). Epithelial-to-Mesenchymal Transition Contributes to Immunosuppression in Breast Carcinomas. Cancer Res.77, 3982–3989. 10.1158/0008-5472.can-16-3292

  • CrossRef
  • Google Scholar

• 10

  GalonJ.CostesA.Sanchez-CaboF.KirilovskyA.MlecnikB.Lagorce-PagèsC.et al (2006). Type, Density, and Location of Immune Cells within Human Colorectal Tumors Predict Clinical Outcome. Science313, 1960–1964. 10.1126/science.1129139

  • CrossRef
  • Google Scholar

• 11

  GeH.YanY.WuD.HuangY.TianF. (2018). Potential Role of LINC00996 in Colorectal Cancer: a Study Based on Data Mining and Bioinformatics. Ott11, 4845–4855. 10.2147/ott.s173225

  • CrossRef
  • Google Scholar

• 12

  GeeleherP.CoxN.HuangR. S. (2014). pRRophetic: an R Package for Prediction of Clinical Chemotherapeutic Response from Tumor Gene Expression Levels. PLoS One9, e107468. 10.1371/journal.pone.0107468

  • CrossRef
  • Google Scholar

• 13

  GreletS.LinkL. A.HowleyB.ObellianneC.PalanisamyV.GangarajuV. K.et al (2017). A Regulated PNUTS mRNA to lncRNA Splice Switch Mediates EMT and Tumour Progression. Nat. Cel Biol19, 1105–1115. 10.1038/ncb3595

  • CrossRef
  • Google Scholar

• 14

  GugnoniM.CiarrocchiA. (2019). Long Noncoding RNA and Epithelial Mesenchymal Transition in Cancer. Int. J. Mol. Sci.20. 10.3390/ijms20081924

  • CrossRef
  • Google Scholar

• 15

  GuoK.FengY.ZhengX.SunL.WasanH. S.RuanS.et al (2021). Resveratrol and its Analogs: Potent Agents to Reverse Epithelial-To-Mesenchymal Transition in Tumors. Front. Oncol.11, 644134. 10.3389/fonc.2021.644134

  • CrossRef
  • Google Scholar

• 16

  GuoY.YangP. T.WangZ. W.XuK.KouW. H.LuoH. (2020). Identification of Three Autophagy-Related Long Non-coding RNAs as a Novel Head and Neck Squamous Cell Carcinoma Prognostic Signature. Front. Oncol.10, 603864. 10.3389/fonc.2020.603864

  • CrossRef
  • Google Scholar

• 17

  HänzelmannS.CasteloR.GuinneyJ. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics14, 7.

  • Google Scholar

• 18

  HuX.-T.XingW.ZhaoR.-S.TanY.WuX.-F.AoL.-Q.et al (2020). HDAC2 Inhibits EMT-Mediated Cancer Metastasis by Downregulating the Long Noncoding RNA H19 in Colorectal Cancer. J. Exp. Clin. Cancer Res.39, 270. 10.1186/s13046-020-01783-9

  • CrossRef
  • Google Scholar

• 19

  IslingerM.VoelklA.FahimiH. D.SchraderM. (2018). The Peroxisome: an Update on Mysteries 2.0. Histochem. Cel Biol150, 443–471. 10.1007/s00418-018-1722-5

  • CrossRef
  • Google Scholar

• 20

  JiaY.DuanY.LiuT.WangX.LvW.WangM.et al (2019). LncRNA TTN-AS1 Promotes Migration, Invasion, and Epithelial Mesenchymal Transition of Lung Adenocarcinoma via Sponging miR-142-5p to Regulate CDK5. Cell Death Dis10, 573. 10.1038/s41419-019-1811-y

  • CrossRef
  • Google Scholar

• 21

  JiangH. Y.NingG.WangY. S.LvW. B. (2020). 14-CpG-Based Signature Improves the Prognosis Prediction of Hepatocellular Carcinoma Patients. Biomed. Res. Int.2020, 9762067. 10.1155/2020/9762067

  • CrossRef
  • Google Scholar

• 22

  JiangY.ZhanH. (2020). Communication between EMT and PD-L1 Signaling: New Insights into Tumor Immune Evasion. Cancer Lett.468, 72–81. 10.1016/j.canlet.2019.10.013

  • CrossRef
  • Google Scholar

• 23

  LambertA. W.PattabiramanD. R.WeinbergR. A. (2017). Emerging Biological Principles of Metastasis. Cell168, 670–691. 10.1016/j.cell.2016.11.037

  • CrossRef
  • Google Scholar

• 24

  LamouilleS.XuJ.DerynckR. (2014). Molecular Mechanisms of Epithelial-Mesenchymal Transition. Nat. Rev. Mol. Cel Biol15, 178–196. 10.1038/nrm3758

  • CrossRef
  • Google Scholar

• 25

  LiL.WangT.LiS.ChenZ.WuJ.CaoW.et al (2020). TDO2 Promotes the EMT of Hepatocellular Carcinoma through Kyn-AhR Pathway. Front. Oncol.10, 562823. 10.3389/fonc.2020.562823

  • CrossRef
  • Google Scholar

• 26

  LinW.-H.ChangY.-W.HongM.-X.HsuT.-C.LeeK.-C.LinC.et al (2021). STAT3 Phosphorylation at Ser727 and Tyr705 Differentially Regulates the EMT-MET Switch and Cancer Metastasis. Oncogene40, 791–805. 10.1038/s41388-020-01566-8

  • CrossRef
  • Google Scholar

• 27

  LinaS. (2021). Identification of Hub lncRNAs in Head and Neck Cancer Based on Weighted Gene Co‐expression Network Analysis and Experiments. FEBS Open Bio11, 2060–2073. 10.1002/2211-5463.13134

  • CrossRef
  • Google Scholar

• 28

  LuW.ZhangH.NiuY.WuY.SunW.LiH.et al (2017). Erratum to: Long Non-coding RNA Linc00673 Regulated Non-small Cell Lung Cancer Proliferation, Migration, Invasion and Epithelial Mesenchymal Transition by Sponging miR-150-5p. Mol. Cancer16, 144. 10.1186/s12943-017-0716-6

  • CrossRef
  • Google Scholar

• 29

  LvM.ZhongZ.HuangM.TianQ.JiangR.ChenJ. (2017). lncRNA H19 Regulates Epithelial-Mesenchymal Transition and Metastasis of Bladder Cancer by miR-29b-3p as Competing Endogenous RNA. Biochim. Biophys. Acta (Bba) - Mol. Cel Res.1864, 1887–1899. 10.1016/j.bbamcr.2017.08.001

  • CrossRef
  • Google Scholar

• 30

  MittalV. (2018). Epithelial Mesenchymal Transition in Tumor Metastasis. Annu. Rev. Pathol. Mech. Dis.13, 395–412. 10.1146/annurev-pathol-020117-043854

  • CrossRef
  • Google Scholar

• 31

  NewmanA. M.LiuC. L.GreenM. R.GentlesA. J.FengW.XuY.et al (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods12, 453–457. 10.1038/nmeth.3337

  • CrossRef
  • Google Scholar

• 32

  NietoM. A.HuangR. Y.-J.JacksonR. A.ThieryJ. P. (2016). EMT: 2016. Cell166, 21–45. 10.1016/j.cell.2016.06.028

  • CrossRef
  • Google Scholar

• 33

  PastushenkoI.BlanpainC. (2019). EMT Transition States during Tumor Progression and Metastasis. Trends Cel Biol.29, 212–226. 10.1016/j.tcb.2018.12.001

  • CrossRef
  • Google Scholar

• 34

  PengP.-H.Chieh-Yu LaiJ.HsuK.-W.WuK.-J. (2020). Hypoxia-induced lncRNA RP11-390F4.3 Promotes Epithelial-Mesenchymal Transition (EMT) and Metastasis through Upregulating EMT Regulators. Cancer Lett.483, 35–45. 10.1016/j.canlet.2020.04.014

  • CrossRef
  • Google Scholar

• 35

  ShiraziF.JonesR. J.SinghR. K.ZouJ.KuiatseI.BerkovaZ.et al (2020). ActivatingKRAS,NRAS, andBRAFmutants Enhance Proteasome Capacity and Reduce Endoplasmic Reticulum Stress in Multiple Myeloma. Proc. Natl. Acad. Sci. USA117, 20004–20014. 10.1073/pnas.2005052117

  • CrossRef
  • Google Scholar

• 36

  SmithH. L.SouthgateH.TweddleD. A.CurtinN. J. (2020). DNA Damage Checkpoint Kinases in Cancer. Expert Rev. Mol. Med.22, e2. 10.1017/erm.2020.3

  • CrossRef
  • Google Scholar

• 37

  SrinivasanP.WuX.BasuM.RossiC.SandlerA. D. (2018). PD-L1 Checkpoint Inhibition and Anti-CTLA-4 Whole Tumor Cell Vaccination Counter Adaptive Immune Resistance: A Mouse Neuroblastoma Model that Mimics Human Disease. Plos Med.15, e1002497. 10.1371/journal.pmed.1002497

  • CrossRef
  • Google Scholar

• 38

  SunZ.JingC.XiaoC.LiT. (2020). An Autophagy-Related Long Non-coding RNA Prognostic Signature Accurately Predicts Survival Outcomes in Bladder Urothelial Carcinoma Patients. Aging12, 15624–15637. 10.18632/aging.103718

  • CrossRef
  • Google Scholar

• 39

  TianC.LiX.GeC. (2021). High Expression of LAMA3/AC245041.2 Gene Pair Associated with KRAS Mutation and Poor Survival in Pancreatic Adenocarcinoma: a Comprehensive TCGA Analysis. Mol. Med.27, 62. 10.1186/s10020-021-00322-2

  • CrossRef
  • Google Scholar

• 40

  VriezeS. I. (2012). Model Selection and Psychological Theory: a Discussion of the Differences between the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC). Psychol. Methods17, 228–243. 10.1037/a0027127

  • CrossRef
  • Google Scholar

• 41

  WangL.ChoK. B.LiY.TaoG.XieZ.GuoB. (2019). Long Noncoding RNA (lncRNA)-Mediated Competing Endogenous RNA Networks Provide Novel Potential Biomarkers and Therapeutic Targets for Colorectal Cancer. Int. J. Mol. Sci.20. 10.3390/ijms20225758

  • CrossRef
  • Google Scholar

• 42

  WangS.ChaiK.ChenJ. (2019). A Novel Prognostic Nomogram Based on 5 Long Non-coding RNAs in clear Cell Renal Cell Carcinoma. Oncol. Lett.18, 6605–6613. 10.3892/ol.2019.11009

  • CrossRef
  • Google Scholar

• 43

  WuX.LvD.CaiC.ZhaoZ.WangM.ChenW.et al (2020). A TP53-Associated Immune Prognostic Signature for the Prediction of Overall Survival and Therapeutic Responses in Muscle-Invasive Bladder Cancer. Front. Immunol.11, 590618. 10.3389/fimmu.2020.590618

  • CrossRef
  • Google Scholar

• 44

  WuY.ZhangL.HeS.GuanB.HeA.YangK.et al (2020). Identification of Immune-Related LncRNA for Predicting Prognosis and Immunotherapeutic Response in Bladder Cancer. Aging (Albany NY)12, 23306–23325. 10.18632/aging.104115

  • CrossRef
  • Google Scholar

• 45

  XingL.ZhangX.ChenA. (2019). Prognostic 4-lncRNA-Based Risk Model Predicts Survival Time of Patients with Head and Neck Squamous Cell Carcinoma. Oncol. Lett.18, 3304–3316. 10.3892/ol.2019.10670

  • CrossRef
  • Google Scholar

• 46

  XuF.ShangguanX.PanJ.YueZ.ShenK.JiY.et al (2021). HOXD13 Suppresses Prostate Cancer Metastasis and BMP4 ‐induced Epithelial‐mesenchymal Transition by Inhibiting SMAD1. Int. J. Cancer148, 3060–3070. 10.1002/ijc.33494

  • CrossRef
  • Google Scholar

• 47

  XuY.XuF.LvY.WangS.LiJ.ZhouC.et al (2021). A ceRNA-Associated Risk Model Predicts the Poor Prognosis for Head and Neck Squamous Cell Carcinoma Patients. Sci. Rep.11, 6374. 10.1038/s41598-021-86048-x

  • CrossRef
  • Google Scholar

• 48

  YangC.HanC.WangX.LiaoX.LiuX.QinW.et al (2018). Clinical Implication and the Hereditary Factors of NM23 in Hepatocellular Carcinoma Based on Bioinformatics Analysis and Genome-wide Association Study. J. Oncol.2018, 6594169. 10.1155/2018/6594169

  • CrossRef
  • Google Scholar

• 49

  YangW.SoaresJ.GreningerP.EdelmanE. J.LightfootH.ForbesS.et al (2013). Genomics of Drug Sensitivity in Cancer (GDSC): a Resource for Therapeutic Biomarker Discovery in Cancer Cells. Nucleic Acids Res.41, D955–D961. 10.1093/nar/gks1111

  • CrossRef
  • Google Scholar

• 50

  YiL.WuG.GuoL.ZouX.HuangP. (2020). Comprehensive Analysis of the PD-L1 and Immune Infiltrates of m6A RNA Methylation Regulators in Head and Neck Squamous Cell Carcinoma. Mol. Ther. - Nucleic Acids21, 299–314. 10.1016/j.omtn.2020.06.001

  • CrossRef
  • Google Scholar

• 51

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun.4, 2612. 10.1038/ncomms3612

  • CrossRef
  • Google Scholar

• 52

  ZhangC.HaoY.WangY.XuJ.TengY.YangX. (2018). TGF-β/SMAD4-Regulated LncRNA-LINP1 Inhibits Epithelial-Mesenchymal Transition in Lung Cancer. Int. J. Biol. Sci.14, 1715–1723. 10.7150/ijbs.27197

  • CrossRef
  • Google Scholar

• 53

  ZhangY.LinA.LiY.DingW.MengH.LuoP.et al (2020). Age and Mutations as Predictors of the Response to Immunotherapy in Head and Neck Squamous Cell Cancer. Front. Cel Dev. Biol.8, 608969. 10.3389/fcell.2020.608969

  • CrossRef
  • Google Scholar