An Epithelial–Mesenchymal Transition Hallmark Gene-Based Risk Score System in Head and Neck Squamous-Cell Carcinoma

Introduction

Head and neck cancer is the seventh most common malignancy worldwide, with the greater being head and neck squamous-cell carcinoma (HNSC).¹ Despite advances in treatments, the prognosis of HNSC remains poor with mortality rates of 40–50%.² HNSC is a heterogeneous type of disease at phenotypic and genetic levels.³ The current clinical decision-making system in HNSC is mainly based on phenotypic heterogeneity, such as the American Joint Committee on Cancer TNM staging system⁴ and tumor grade. It is imperative to identify individuals at high-risk with the same phenotypes by revealing the genetic heterogeneity.

Epithelial-to-mesenchymal transition (EMT) is considered one of the hallmarks of cancer.⁵ The emerging evidence has shown that EMT program contributes to the induction of cancer stem cells, immune escape during cancer progression, and drug resistance in various types of cancers,^6–8 including HNSC.⁹ Thus, the expression pattern of different EMT hallmark genes may be one of the critical genetic heterogeneity of cancers. We assumed that there is an EMT hallmark gene-based risk score system as a prognostic indicator in HNSC. To investigate the hypothesis, we used the datasets of HNSC from TCGA, including clinical information and gene expression profiles, to generate an EMT hallmark gene-based risk score system for predicting prognosis for patients with HNSC, and validated it on another independent dataset.

Materials and Methods

Data Processing

The RNA sequencing (RNA-seq) data (displayed as raw counts) and clinical information of HNSC in The Cancer Genome Atlas (https://www.cancer.gov/) were downloaded to generate an EMT hallmark gene-based risk score (RS) to predict prognosis. Another HNSC data set GSE65858¹⁰ based on platform of GPL10558 was downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) and used to validate the RS. The RNA-seq data was normalized using voom function from limma package¹¹ in R software (version 4.0.2) (https://www.r-project.org/). The gene expression profiles in GSE65858 were normalized by the contributor. If multiple probes correspond to a gene, then the average value of these probes is considered as the expression value of this gene. EMT hallmark gene set included 200 genes (Table S1) was obtained from the Molecular Signatures Database (version 7.2).^12,13 The present study only included 192 EMT hallmark genes because their expression values are available in both TCGA and GSE65858. The workflow of the present study is displayed in Figure 1.

Figure 1 The workflow of the present study. Abbreviations: TCGA, The Cancer Genome Atlas; HNSC, head and neck squamous-cell carcinoma; LASSO, least absolute shrinkage and selection operator; GSEA, gene set enrichment analysis.

Screening Differentially Expressed Genes (DEGs)

The expression profiles of the 192 EMT hallmark genes were extracted from the TCGA HNSC dataset, and subsequently used to screen the DEGs in HNSC compared to healthy tissues using limma package. Genes with P adjusted by false discovery rate < 0.05 and log2 (fold change) >1 were considered significant.

Cox Regression and Least Absolute Shrinkage and Selection Operator Analysis

The expression profiles of DEGs and overall survival data were used to perform univariable Cox regression to identify the prognosis-associated EMT hallmark genes. Subsequently, the expression profiles of prognosis-associated EMT hallmark genes were performed with least absolute shrinkage and selection operator (LASSO) analysis using glmnet (https://CRAN.R-project.org/package=glmnet) R package to select the optimal prognostic EMT hallmark genes. Thus, the EMT hallmark gene-based risk score (RS) was created as:

RS = Expr_gene1*Coef_gene1 + Expr_gene2*Coef_gene2+ Expr_gene2*Coef_gene2+ …

The “Coef” is the regression coefficient of gene and is derived from the LASSO Cox regression, and “Expr” indicates the expression values of the gene. Each patient with HNSC got a RS, and was divided into the high- or low-risk group according to the median RS. The OS between the two different risk groups were compared using log-rank method. In addition, the routine clinical factors were included in the multivariable Cox regression analysis to assess whether the RS is an independent prognostic factor.

Gene Set Enrichment Analysis (GSEA)

To explore the biological state of high-risk group patients, GSEA^12,14 was performed using the GSEA JAVA program (version 4.0.1) (http://www.broadinstitute.org/gsea). The hallmark gene set included 50 gene sets obtained from the Molecular Signatures Database (version 7.2)^12,13 used as the reference gene set. Gene sets with nom p < 0.05 after performing 1000 permutations were considered to be significantly enriched.

Validation of the EMT Hallmark Gene-Based RS

As it was in the TCGA HNSC data set, each patient in GSE65858 got an RS according to the above formula, and was divided into high- or low-risk groups. The OS between the two different risk groups were compared and the multivariable Cox regression analysis was carried out.

Statistical Analysis

All analyses were performed using R software (version 4.0.2). The unpaired t-test from limma package was used to screen DEGs. Kaplan–Meier survival analysis and Log rank test were used to compare survival between the two groups of patients. Time-receiver operating characteristic (tROC) curve analysis was performed using the timeROC package (https://CRAN.R-project.org/package=timeROC). All tests were two-sided and p < 0.05, unless otherwise stated, was considered to indicate statistical significance.

Results

Multiple EMT Hallmark Genes Upregulated in HNSC

A total of 93 EMT hallmark genes (Table S2) were differentially expressed in HNSC compared to healthy paracancer tissue, including 16 downregulated and 77 upregulated genes (Figure 2A). This indicates that EMT plays a crucial role in promoting HNSC due to the fact that most EMT hallmark genes are upregulated. The expression heat map of the DEGs shows that it has a promising effect of distinguishing tumor from paracancer tissue (Figure 2B).

Figure 2 Differentially expressed EMT hallmark genes. (A) The volcano plot of the differentially expressed genes. Red indicates up-regulated, and blue indicates down-regulated. (B) The expression heat map of the differentially expressed EMT hallmark genes.

EMT Hallmark Gene-Based RS as an Independent Prognostic Factor

After univariable Cox regression analysis, nineteen EMT hallmark genes were identified as the prognosis-associated genes (Table 1). Unsurprisingly, most of (17 from 19) the prognosis-associated EMT hallmark genes show an association with poor prognosis in HNSC. Subsequently, nine EMT hallmark genes (SFRP1, TGFBR3, DKK1, PCOLCE2, PTX3, CAP2, PLOD2, VEGFC, and IL6) were considered as optimal features through LASSO Cox analysis (Figure 3A). Thus, all patients got RS according to the coefficients. The RS is significantly associated with poor prognosis (Hazard Rate (HR) = 3.254, 95% CI = 2.367–4.473, p < 0.001). The RS showed promising prognostic value with AUC approximately 0.7 (Figure 3B), and the AUC of 5-year tROC was 0.660 (Figure 3C). The high-risk group HNSC patients showed significantly shorter OS than the low-risk group HNSC patients (Figure 3D). Furthermore, the RS remained independent compared to some routine clinical factors, including TNM staging system, tumor grade, and tumor primary subdivision (Figure 4).

Table 1 The Results of Univariable and Least Absolute Shrinkage and Selection Operator Cox Analysis

Figure 3 The nine-EMT hallmark gene-based risk score in the head and neck squamous-cell carcinoma data set from The Cancer Genome Atlas. (A) Nine genes were considered as the optimal features in the least absolute shrinkage and selection operator Cox analysis. (B) The time-dependent receiver operating characteristic (ROC) curve analysis for the risk score. (C) 5-year The time-dependent ROC curve analysis for the risk score. (D) The Kaplan–Meier curves with the Log rank test of the high- and low-risk groups. Abbreviation: AUC, area under the curve.

Figure 4 The results of the multivariable Cox analysis for the nine-EMT hallmark gene-based risk score (RS) and routine clinical factors. *P < 0.05, **P < 0.01, and ***P < 0.001.

Biological Phenotypes of High-Risk HNSC

The results of GSEA showed that the EMT hallmark gene set was significantly enriched in high-risk HNSC (Figure 5A). In addition to this, hallmark gene set of angiogenesis (Figure 5B), coagulation (Figure 5C), glycolysis (Figure 5D), hypoxia (Figure 5E), MTORC1 signaling (Figure 5F), unfold protein response (Figure 5G), and UV response up (Figure 5H) were also enriched in high-risk HNSC.

Figure 5 The results of gene set enrichment analysis for the high-risk head and neck squamous-cell carcinoma. Eight hallmark gene sets enriched in the samples of high-risk head and neck squamous-cell carcinoma, including (A) epithelial–mesenchymal transition, (B) angiogenesis, (C) coagulation, (D) glycolysis, (E) hypoxia, (F) MTORC1 signaling, (G) unfold protein response, and (H) UV response up.

The RS Was Validated in an Independent Data Set

As it was in the TCGA HNSC data set, the RS was generated for each individual in GSE65858 according to the formula. The RS remained significantly associated with poor prognosis (HR = 13.261, 95% CI = 2.136–82.352, p = 0.006). The high-risk group HNSC patients remained significantly shorter OS than the low-risk group HNSC patients in GSE65858 (Figure 6A). The RS also remained an independent prognostic factor compared to routine clinical factors (Figure 6B).

Figure 6 The nine- EMT hallmark gene-based risk score in GSE65858. (A) The Kaplan–Meier curves with the Log rank test of the high- and low-risk groups. (B) The results of the multivariable Cox analysis for the nine-EMT hallmark gene-based risk score (RS) and routine clinical factors. *P < 0.05, and **P < 0.01.

Discussion

Previous studies constructed a prognostic stratification system from multiple perspectives, such as immunity-related gene-based signature,¹⁵ microRNA-based signature,^16,17 and microenvironment-based system.¹⁸ Few studies focused on EMT-related gene-based signature. In our present study, a nine-EMT hallmark gene-based RS was generated to successfully identify the relatively high-risk HNSCs. It was an independent prognostic factor compared to routine clinical factors; moreover, it was validated in an external data set. This may provide more references for clinical decision-making.

Unsurprisingly, some of the nine EMT hallmark genes reported were associated with HNSC or other types of cancer. CAP2 was considered to be related to multistage hepatocarcinogenesis.¹⁹ Elevated DKK1 was reported as an independent unfavorable prognostic indicator of survival in HNSC,²⁰ which is consistent with our result. IL-6 plays an important role in HNSC tumor proliferation and metastasis, and IL-6/STAT3 signaling may be a potential target for treating HNSCC patients.²¹ PLOD2 was found to contribute to drug resistance in laryngeal cancer by promoting cancer stem cell-like characteristics.²² PTX3 is an extrinsic oncosuppressor regulating complement-dependent inflammation in cancer.²³ The loss of SFRP1 expression is associated with colorectal cancer, prostate cancer, and renal cell cancer.^24,25 Our analysis found that it may also play a tumor suppressor role in HNSC in association with a good prognosis. Gene TGFBR3 was found to play a dual role in bladder cancer, acting as both a tumor suppressor and as a tumor promoter.²⁶ However, a previous study showed it can block lymph node metastasis in head and neck cancer.²⁷ Our analysis shows that TGFBR3 is a protective gene in HNSC. VEGFC may contribute to HNSC growth and motility.²⁸ There are few reports of functional experiments regarding PCOLCE2 in HNSC. In the present study, we found that PCOLCE2 is associated with poor prognosis in HNSC.

Furthermore, we also conducted GSEA to explore the biological characteristics of patients identified as high-risk patients by this nine-EMT hallmark gene-based RS. Based on the results of GSEA, high-risk HNSC is characterized by active EMT program, angiogenesis, MTORC1 signaling, unfold protein response (UPR) program, and high hypoxia. Hypoxia is common in HNSC cells and contributes to malignant behaviors, such as tumor progression, invasion, metastasis, and resistance to chemotherapy and radiotherapy.²⁹ Recent studies suggest that the UPR may affect many hallmarks of cancer, including metastasis, genome stability, angiogenesis, inflammation, and drug resistance.³⁰ These partially explain the reasons for the poor prognosis of the high-risk group. There are numerous ongoing efforts to target the mTOR signaling pathway for cancer therapy,³¹ however, whether active MTORC1 signaling in high-risk HNSC indicates the response to therapies for targeting mTORC1 signaling still needs further exploration.

Although the present study may provide new insight into the prognostic systems in HNSC, it has several noticed limitations. First, not all EMT hallmark genes were included in the analysis due to the data sets coming from different centers. Thus, the RS may be improved in further study. Secondly, the role of some of these nine genes in HNSC is not yet clear, therefore, it is not clear whether these genes are causal or merely markers for predicting prognosis in HNSC. Thirdly, the present study lacks molecular experimental verification of candidate molecules.

In conclusion, we proposed and validated a nine-EMT hallmark gene-based risk score system for predicting prognosis for patients with HNSC, and also preliminarily explained the biological characteristics of high-risk patients.