A prognostic long non-coding RNA-associated competing endogenous RNA network in head and neck squamous cell carcinoma

Introduction

Head and neck squamous cell carcinoma (HNSCC) has the ninth highest global incidence of all malignancies, and accounted for approximately 53,260 new cases and approximately 10,750 deaths in 2019 in the United States of America (Ward et al., 2019; Siegel, Miller & Jemal, 2020; Thrift & Gudenkauf, 2019). Despite advances in surgical procedures, chemo-radiotherapy, and targeted therapy, clinical outcomes have remained unchanged for decades, with the five-year survival rate ranging from 40% to 50% over the past three decades (Braakhuis, Leemans & Visser, 2014). HNSCC is characterized by a remarkably aggressive biological behavior with a high incidence of lymph node and distant metastasis that leads to poor prognosis (Ding et al., 2019; Sun et al., 2017). The 8th edition of American Joint Commission Cancer (AJCC) introduced depth of invasion and extranodal extension as criteria for T and N stages, respectively (Amin et al., 2017). TNM staging are useful for determination of prognosis and development of strategic treatment plans based on the anatomical situation (Asare et al., 2015). However, even patients with the same TNM stage require different therapeutic approaches. Therefore, identification of additional molecular biomarkers is required to improve prognosis by accounting for the biological heterogeneity of HNSCC.

Differential expression of long non-coding RNAs (lncRNAs) might affect the corresponding function and play an important role in cancer pathogenesis. LncRNAs and microRNAs (miRNAs) are important subclasses of competing endogenous RNAs (ceRNAs), which regulate the protein levels of target genes, thus, participating in all aspects of various cellular biological processes (Lin & Gregory, 2015; Schmitt & Chang, 2016). Accumulating evidence suggests that the ceRNA network could be beneficial for prognosis prediction in patients with HNSCC (Xu et al., 2019; Yin et al., 2019; Pan et al., 2019; Fang et al., 2018; Zhao et al., 2018). Zhao et al. (2018) established a ceRNA regulation network related to HNSCC that included eight differently expressed mRNAs (DEmRNAs), 53 differentially expressed lncRNAs (DElncRNAs), and 16 differentially expressed miRNAs (DEmiRNAs). However, that study and others did not construct a prognostic model for HNSCC that assessed the relationship between ceRNA network RNAs and overall survival (OS). In addition, several important characteristics are closely related to the clinicopathological features of HNSCC, such as TP53 mutation status, angiolymphatic invasion (ALI), human papillomavirus (HPV) status, perineural invasion (PNI), and extracapsular spread (ECS). However, no previous studies have, to the best of our knowledge, assessed the relationship between ceRNA network RNAs and the above clinicopathological characteristics.

In the present study, we comprehensively evaluated the prognostic value of differentially expressed RNAs and developed multi-RNA-based models that predicted OS via the risk signature of HNSCC. In addition, we systematically investigated the relationship between prognostic RNAs and HNSCC clinicopathological features.

Materials & Methods

Differentially expressed RNAs related to HNSCC

The design and analytical approach of the present study are shown in the flow chart (Figs. 1A and 1B). lncRNA-seq and mRNA-seq data were obtained from 546 samples consisting of 502 HNSCC and 44 adjacent normal tissues. MiRNA-seq information was obtained from 569 samples consisting of 525 HNSCC and 44 adjacent normal tissue samples. According to the cutoff criteria of —log2 (fold-change [FC]) — > 2 and adjusted p < 0.05, DEmRNAs and DEmiRNAs in the HNSCC and normal tissues were identified using the “edgeR” package in R. Next, DElncRNAs were screened out using the same method above in R with the thresholds of —log₂FC— > 2.0 and adjusted p-value < 0.05. Then, the DElncRNAs were annotated and defined using the Encyclopedia of DNA Elements (ENCODE). The statistical significance of multiple comparisons was corrected using FDR (false discovery rate) and FDR < 0.05 was considered statistically significant.

Results

Identification of DElncRNAs, DEmRNAs and DEmiRNAs

Each sample from patients with HNSCC was analyzed to obtain the expression profiles of lncRNA, mRNA, and miRNA for comprehensive integrated analysis. A total of 1052 DElncRNAs, 82 DEmiRNAs, and 2001 DEmRNAs were differentially expressed between HNSCC and normal tissue samples. Of these, 766 lncRNAs (72.81%) were upregulated and 286 (27.19%) were downregulated; 873 mRNAs (43.63%) were upregulated and 1128 (56.37%) downregulated; and 44 miRNAs (53.66%) were upregulated and 38 (46.34%) were downregulated. The “gplots” package was used to plot the volcanos and heat maps of differentially expressed RNAs (Figs. 2A–2C, Figs. S1A–S1C).

Functional analysis

GO enrichment analysis showed that on the BP, the differential expression genes were associated with “keratinization”, “anterior/posterior pattern specification”, etc. For MF, the genes were related to “receptor ligand activity”, “receptor regulator activity”, etc. CC were closely correlated with “intermediate filament”, “intermediate filament cytoskeleton”, etc (Fig. 3A). KEGG enrichment analysis revealed that mRNAs were dramatically enriched in pathways related to “transcriptional misregulation in cancer” and “neuroactive ligand–receptor interaction”, which are significantly correlated with carcinogenesis and metastasis (Fig. 3B). Further, correlations of the top 185 most significant mRNAs (p < 0.05) were analyzed by constructing PPI networks (Fig. 3C). The results showed that these mRNAs play key roles during the initiation and progression of HNSCC.

Predictive model for OS

In total, 19 lncRNAs, 18 miRNAs, and 25 mRNAs identified in univariate Cox regression analysis were selected for the LASSO regression analysis. Of these, 13 lncRNAs (LINC00460, AL136987.1, MYOSLID, MIR9-3HG, AC073130.1, AC079160.1, LINC01305, AP002478.1, LINC02434, HOTTIP, ATP6V1B1-AS1, AC023310.4, and AL158209.1), eight miRNAs (hsa-miR-411, hsa-miR-4510, hsa-miR-410, hsa-miR-99a, hsa-miR-499a, hsa-miR-4652, hsa-miR-206 and hsa-miR-520e), and 17 mRNAs (CELSR3, STC2, ADGRD2, ZNF541, GRB14, NOSTRIN, TIMP4, HOXB9, ADPRHL1, SPINK1, PTX3, FRZB, ODF4, PLAU, OLR1, GNG7, and DTHD1) were independent prognostic factors related to OS in patients with HNSCC. Associations between each of the above-mentioned RNAs and clinicopathological features and molecular characteristics were assessed (Tables S2–S4). Based on the powerful prognostic RNAs, three risk signatures were constructed (Figs. 4A–4F). LASSO regression analysis of the 13 lncRNAs, eight miRNAs, and 17 mRNAs was conducted, and the RS was calculated according to the coefficients shown in Table 1. To evaluate the prognostic performance of the risk signature, patients with HNSCC were divided into low- and high-risk groups based on the median RS. The results suggested that patients in the high-risk group had worse prognoses than those in the low-risk group (Figs. 4G–4I). In addition, specificity and sensitivity were assessed using a time-dependent ROC curve, and the AUC values were 0.62, 0.619, and 0.688, respectively, suggesting better prediction performance (Figs. 4J–4L).

Table 1:

The coefficients estimated by LASSO regression analysis.

DElncRNA	Coef		DEmiRNA	Coef		DEmRNA	Coef
LINC00460	2.1735E−04		has-miR-411	2.0123E−04		CELSR3	−8.2570E−05
AL136987.1	2.3890E−03		has-miR-4510	−3.0047E−02		STC2	3.8044E−05
MYOSLID	1.6926E−04		has-miR-410	1.9585E−04		ADGRD2	−6.0742E−04
MIR9-3HG	−2.8041E−05		has-miR-99a	−6.1385E−05		ZNF541	−1.2459E−04
AC073130.1	4.1471E−04		has-miR-499a	1.7026E−03		GRB14	2.7241E−04
AC079160.1	1.9663E−03		has-miR-4652	2.6325E−03		NOSTRIN	−3.0643E−03
LINC01305	−7.1734E−07		has-miR-206	1.8113E−06		TIMP4	1.6223E−03
AP002478.1	2.1455E−03		has-miR-520e	1.1725E−02		HOXB9	9.9531E−05
LINC02434	1.6812E−03					ADPRHL1	3.1718E−04
HOTTIP	4.2863E−04					SPINK1	3.3189E−04
ATP6V1B1-AS1	1.1889E−03					PTX3	1.1693E−04
AC023310.4	6.1627E−04					FRZB	−3.4147E−05
AL158209.1	2.8547E−02					ODF4	−1.9280E−02
						PLAU	5.3363E−06
						OLR1	1.1427E−05
						GNG7	−3.3920E−04
						DTHD1	−1.2034E−03

DOI: 10.7717/peerj.9701/table-1

Notes:

DElncRNA

differentially expressed long noncoding RNA

DEmRNA

differentially expressed messenger RNA

DEmiRNA

differentially expressed microRNA

Coef

coefficient

The expression levels of differentially expressed RNAs related to clinicopathological characteristics, including the three risk signatures, are shown in the heat map (Figs. 5A–5C). Classification into the high-risk group was closely related to patient died status, TP53 mutation status, a high T stage, HPV negative, PNI, and disease recurrence/progression in the lncRNA-based model; to patient died status, TP53 mutation status, high M stage, and HPV negative in the miRNA-based model; and to patient died status, TP53 mutation status, high clinical stage, a high M stage, HPV negative, PNI, ECS, and disease recurrence/progression in the mRNA-based model (Figs. 5A–5C, Table 2). Thus, age, stage, T stage, N stage and the risk score were strongly related to OS in univariate analysis (Figs. 6A–6C); risk score and age, in the lncRNA and mRNA models; and risk score, age, and N stage, in the miRNA model (multivariate analysis) (Figs. 6D–6F) (p < 0.05). These results showed that the risk score could be regarded as an independent prognostic factor in HNSCC.

Verification of the prognostic values of ceRNA

The verified results of GEPIA showed that HOTTIP, LINC00460, and STC2 were significantly associated with OS, strengthening the reliability of our findings (Figs. 8A–8P). Notably, the expression of hsa-miR-206 (GSE45238) and STC2 (GSE2379, GSE13601, GSE25099, GSE12452 and GSE9844) in six microarrays was significantly differential expression in HNSCC tissues than in normal tissues (Figs. 8Q–8V). Furthermore, qRT-PCR results showed that HOTTIP and LINC00460 were upregulated (Figs. 8W and 8X), hsa-miR-206 was downregulated (Fig. 8Y) to varying degrees in many HNSCC cell lines (HN4, HN6, Cal27, HN30, SCC9, and SCC25).

Table 2:

Clinicopathological features are different between high-risk and low-risk in DElncRNA, DEmiRNA and DEmRNA.

Clinicopathological features		DElncRNA			DEmiRNA			DEmRNA
		High-risk	Low-risk	p	High-risk	Low-risk	p	High-risk	Low-risk	p
HPV p16	Negative	29	37	9.11e−3**	42	22	1.67e−5***	35	30	4.60e−5***
	Positive	4	25		4	24		2	27
HPV ish	Negative	23	30	8.19e−3**	34	18	6.87e−5***	27	25	2.95e−4***
	Positive	1	17		1	16		0	18
Perineural Invasion	Yes	93	62	0.01*	88	64	0.17	103	51	8.34e−6***
	No	80	95		87	88		73	102
Angiolymphatic Invasion	Yes	58	51	0.73	62	44	0.25	65	44	0.12
	No	106	104		107	103		104	105
Extracapsular spread pathologic	Gross Extension	20	12	0.34	15	16	0.22	19	13	6.12e−3**
	Microscopic Extension	35	35		43	26		47	22
	No Extranodal Extension	110	116		114	110		106	120
Lymph nodes counts	<18	43	34	0.63	41	37	1	46	31	0.20
	≥18	161	149		160	145		157	152
Disease Free Status	Disease Free	88	138	1.25e−3**	105	121	0.64	73	155	1.20e−8***
	Recurred/Progressed	78	59		67	68		86	50

DOI: 10.7717/peerj.9701/table-2

Notes:

DElncRNA

differentially expressed long noncoding RNA

DEmRNA

differentially expressed messenger RNA

DEmiRNA

differentially expressed microRNA

*p <  0.05.

**p <  0.01.

***p <  0.001.

Discussion

Determining the potential molecular mechanisms of the onset and development of HNSCC and new prognostic targets will help improve the survival rate of patients. To date, few studies have reported single or multiple biomarkers for HNSCC due to its multiple genetic alternations and mutations and their complicated interactions. Therefore, identification of specific and sensitive lncRNA biomarkers and prognostic targets for HNSCC is urgently needed (Guglas et al., 2017). ceRNA regulatory networks are composed of lncRNAs, miRNAs, and mRNAs and have attracted great interest in recent years with respect to their contribution to the molecular mechanisms underlying tumor development and subsequent malignant progression (Deng et al., 2015). The present study explored how the HNSCC-specific lncRNAs function as ceRNAs to regulate their target genes and constructed RNA-based prognostic models for HNSCC.

Many studies have implicated the ceRNA regulatory network comprising lncRNAs, miRNAs, and mRNAs in the development of HNSCC (Xu et al., 2019; Yin et al., 2019; Pan et al., 2019; Fang et al., 2018; Zhao et al., 2018). All the above studies are based on data from the TCGA database for the construction of a ceRNA regulatory network and to obtain survival-related RNAs of HNSCC. The present study screened HOTTIP, LINC00460, hsa-miR-206, and STC2 and, similar to the results of previous studies, a comparison with which indicated that our research methods were credible. However, these studies only assessed the relationship between signal lncRNA/miRNA/mRNA and survival and did not construct a prognostic model. Conventional methods are often directly used to screen for differentially expressed RNAs while constructing a ceRNA regulatory network. In contrast to other studies, in this study, we first analyzed differentially expressed RNAs to identify prognosis-related RNAs using univariate Cox regression analysis, and then selected the prognosis-related DElncRNAs, DEmiRNAs, and DEmRNAs for the construction of the ceRNA network. This approach ensured that all RNAs incorporated into the network were associated with survival.

Another study by Yin et al. (2019) showed that 71 lncRNAs, eight miRNAs, and 16 mRNAs were established an HNSCC-specific ceRNA network. HOTTIP was proven to be the most powerful prognostic factor and was markedly associated with clinical stage and histological grade in three patients with HNSCC. Moreover, our study used univariate Cox regression analysis and LASSO analysis to establish RNA-based risk signatures. Compared with multivariate Cox regression analysis, LASSO analysis may be more beneficial in reducing the occurrence of data model overfitting. However, a previous study identified only lncRNAs and mRNA prognostic signatures among the ceRNA network and could not establish a miRNA-based model. Forthermore, we analyzed the relationship between the prognostic model and different clinicopathological features, and thereby strictly explained the predictive role of risk signatures. More importantly, our study not only established three prognostic models consisting of prognostic RNAs, but also assessed the relationship between these prognostic models and conventional clinicopathological features (including TNM stage and histological grade) and specific HNSCC clinicopathological features (including TP53 mutation status, ALI, HPV status, PNI, and ECS).

Previous publications have focused in part on ceRNAs in laryngeal squamous cell carcinoma (LSCC) (Kong et al., 2019; Liu & Ye, 2019; Gao et al., 2019; Zhao et al., 2019; Sun et al., 2018; Feng et al., 2016), tongue squamous cell carcinoma (TSCC) (Sturgis & Cinciripini, 2007), and nasopharyngeal carcinoma (NPC) (Fullerton et al., 2020; Agrawal et al., 2011). These studies have identified differentially expressed RNAs by lncRNA and mRNA integrated microarrays of cancer tissue or by immunohistochemistry (Gao et al., 2019; Zhao et al., 2019; Sun et al., 2018; Feng et al., 2016). The different key RNAs screened may be attributed to the patients included from different databases and the different tumor types. The evaluation ability was 0.62 for the lncRNA-based model, 0.619 for the miRNA-based model, and 0.688 for the mRNA-based model.

Importantly, while establishing prognosis-related prediction models, we also analyzed the relationship between the three prognostic models and clinicopathological features, including the relationship between the expression of each independent RNA in the model and clinicopathological features. This has, to the best of our knowledge, never been explored. Our results showed that the three prognostic models significantly correlated with clinicopathological features, such as OS, DFS, TP53 mutation status, and HPV status. Furthermore, the RS of the three risk signatures as an independent predictor of prognosis plays an important role in prognostic prediction in HNSCC. Notably, HPV is regarded as a crucial risk factor for the development of oropharyngeal cancer (Sturgis & Cinciripini, 2007). New studies have shown that patients with HPV-positive and HPV-negative oropharyngeal cancer have dramatically different rates of both head and neck carcinoma mortality and competing-cause mortality. Compared with HPV-positive patients, HPV-negative patients have a higher risk of two-year cumulative incidence of all-cause mortality and a lower risk of both head and neck carcinoma-specific mortality and competing-cause mortality (p < 0.0001) (Fullerton et al., 2020). Therefore, we believe that the HPV status is a critical determinant in oropharyngeal cancer. TP53 mutation was observed in 60% to 80% of HPV-negative HNSCC patients and was associated with poor prognosis, increased resistance to standard therapy (mainly cisplatin and radiation), and the shortest time to distant metastases (Agrawal et al., 2011; Stransky et al., 2011; Neskey et al., 2015). The significance of TP53 mutations is in predicting the efficacy of chemotherapy, targeted therapy, and radiotherapy.

In our KEGG pathway analysis, the prognostic mRNAs were significantly enriched in “Transcriptional misregulation in cancer” and “Neuroactive ligand–receptor interaction”. The pathways related to “Transcriptional misregulation in cancer” included the P53 signaling pathway, renal cell carcinoma, thyroid cancer, and acute myeloid leukemia. Xu et al. (2019) reported the same KEGG pathways to be enriched in HNSCC.

We further validated the RNAs in the ceRNA regulatory network using the TCGA database, both in vivo and in vitro. HOTTIP and LINC00460 were significantly associated with OS validated by GTEx data and several gene microarrays. They also showed varying degrees of upregulation in six HNSCC cell lines. HOTTIP is a newly identified lncRNA encoded by chromosome 7p15.2 located at the 5 ′-end of the HOXA cluster (Wang et al., 2011). Many studies have reported that higher expression of HOTTIP is associated with positive lymph node metastasis and unfavorable prognosis in various cancer types (Lian et al., 2016; Chen et al., 2017). Our previous studies have confirmed that HOTTIP is also involved in the construction of a ceRNA regulatory network in renal clear-cell carcinoma (Wang et al., 2018). Zhang et al. (2015) reported that the overexpression of HOTTIP is an independent poor prognostic factor for patients with TSCC. Xue et al. (2019) showed that the downregulation of LINC00460 can facilitate cancer cell apoptosis and autophagy in HNSCC. Our previous studies have confirmed that the lncRNA signature and tumor grade are independent prognostic factors and that the 3-lncRNA (KTN1-AS1, LINC00460, and RP5-894A10.6) signature could be a novel prognostic biomarker for HNSCC (Cao et al., 2017). Jiang et al. (2019) proposed that LINC00460 could enhance the proliferation and metastasis of HNSCC cells and promote EMT by accelerating PRDX1 entry into the nucleus. LINC00460 contributes to the progression of NPC by sponging miR-149-5p to upregulate IL6 (Kong et al., 2018).

In this study, only one prognostic miRNA (hsa-miR-206) involved in ceRNAs network. In terms of hsa-miR-206, it has been reported to be downregulated and associated with tumor progression and worse prognosis in several malignancies such as gastric cancer, osteosarcoma and colorectal carcinoma (Arigami et al., 2013; Bao et al., 2013; Liu et al., 2017). In HNSCC, it was demonstrated that hsa-miR-206 plays a critical role progression by targeting HDAC6 via PTEN/AKT/mTOR pathway, which might be a potential target for diagnostic and therapeutic (Vickers et al., 2012). The target prognostic mRNA of hsa-miR-206 was STC2 in ceRNAs regulatory network in present study. STC2 (Stanniocalcin 2) is a secreted glycoprotein with important functions in gastric cancer, which could be a powerful marker of poor prognosis (Yang et al., 2013; Yang et al., 2017). Yang et al. concluded that STC2 may promote HNSCC metastasis via PI3K/AKT/Snail signaling axis. Targeted therapy against STC2 may become a novel strategy to effectively treat patients with metastatic HNSCC (Yokobori et al., 2010).

This study has various advantages and limitations. The advantage of bioinformatics technology is that it can provide novel insights into the potential molecular mechanisms involved in the development and progression of HNSCC (Yu et al., 2014). However, a well-designed and verified classifier model is needed to ensure the reliability of our findings. The highlights of this study are as follows: First, we established three multi-RNA-based prognostic models and a prognostic ceRNA regulatory network for HNSCC, all of which can increase the effectiveness of prognosis prediction. Second, we comprehensively analyzed the relationship between the prognosis model and each RNA included in the prediction model and conventional clinicopathological features (including TNM stage and histological grade) and specific HNSCC clinicopathological features (including TP53 mutation status, ALI, HPV status, PNI, and ECS) to further verify the significant correlation with survival and clinical therapeutic outcomes. lncRNAs can regulate gene expression at several levels and are involved in the evolution and progression of many cancers. However, the present study only identified one regulatory mechanism. Therefore, further clinical investigations and scientific research are required to validate these conclusions.

Conclusions

In conclusion, the present study analyzed an independent cohort of patients with HNSCC from the TCGA database, constructed a ceRNA regulatory network, and established three multi-RNA-based prognostic models with potential prognostic value, and the findings may provide a basis for further experimental and clinical research.

Supplemental Information