Identifying prognostic lncRNAs based on a ceRNA regulatory network in laryngeal squamous cell carcinoma

Growing evidence demonstrates that long non-coding RNAs (lncRNAs) play a crucial role as competing endogenous RNAs (ceRNAs) in tumor occurrence. The lncRNAs’ functions and clinical significance in laryngeal squamous cell carcinoma (LSCC) remain unclear. The study aims to reveal the lncRNA-associated ceRNA regulatory network of LSCC and clarify its clinical relevance. Here, we obtained LSCC transcriptome data from The Cancer Genome Atlas (TCGA) database and identified the differential expression profile of lncRNAs, miRNAs, and mRNAs by the EdgeR R package. The function enrichment analysis of mRNAs was performed using clusterProfiler R package and GSEA3.0. Then, we constructed a ceRNA network and prognosis model based on lncRNAs through bioinformatic methods. Moreover, we explored the functions of prognosis-related lncRNA in LSCC by CCK-8 and transwell assay. 1961 lncRNAs, 69 miRNAs, and 2224 mRNAs were identified as differentially expressed genes in LSCC tissues. According to the transcriptome differential expression profile, a ceRNA network containing 61 lncRNAs, 21 miRNAs, and 77 mRNAs was established. Then, four lncRNAs (AC011933.2, FAM30A, LINC02086, LINC02575) were identified from the ceRNA network to build a prognosis model for LSCC patients. And we found that LINC02086 and LINC02575 promoted the proliferation, migration, and invasion of LSCC cells while AC011933.2 and FAM30A inhibited these biological functions in vitro. Furthermore, we validated that LINc02086/miR-770-5p/SLC26A2 axis promoted migration in LSCC. Four lncRNAs of the ceRNA network were abnormally expressed and related to patient prognosis in LSCC. They played a significant role in the progress of LSCC via affecting the proliferation and metastasis of tumor cells.


Introduction
Laryngeal squamous cell carcinoma (LSCC) is a typical head and neck squamous cell carcinoma [1]. The conventional treatment for LSCC, including chemotherapy, radiation therapy, and surgical resection, have a satisfied curative effect on early-stage patients [2]. However, patients with advanced LSCC continuously have a low overall five-year survival rate [3]. Therefore, it is necessary to reveal the pathogenesis of LSCC and improve the level of diagnosis and treatment.
Non-coding RNAs (ncRNAs) are a type of RNA that cannot encode proteins, including small noncoding RNA (sncRNAs) and long non-coding RNA (lncRNAs) [4,5]. MicroRNAs (miRNAs) that belong to sncRNAs and have been widely studied decreases the expression level of mRNAs by binding to its 3′untranslated regions (3′-UTRs) to degrade mRNAs [6,7]. Long non-coding RNAs (lncRNAs), the length of over 200 nucleotides (nt), can play significant regulation functions in several biological processes, including transcription, pretranscription, chromatin modification, translation, post-translation [8][9][10]. Currently, researchers found that the abnormal expression of lncRNAs is involved in the occurrence and development of malignant tumors [11]. For example, lncRNA HOTAIR is upregulated in many cancers and promotes the proliferation, migration, and invasion of tumor cells, including breast cancer, rectal cancer, pancreatic cancer, and kidney cancer [12][13][14][15]. Upregulated-LINC02410 may serve a diagnostic marker for rectal cancer [16].
The ceRNA hypothesis, a hot topic of lncRNA research, is that some lncRNAs sharing the same miRNA response elements (MREs) with mRNAs bind miRNAs competitively to block the interaction between miRNA and mRNA, thereby the expression level of mRNAs has been elevated [17,18]. It has been reported that the ceRNA regulatory patterns extensively are present in many types of cancers [19]. For example, the lncRNA FAM225A upregulates ITGB3 by adsorbing miR-590-3p/miR-1275 to promote NPC cells' proliferation and invasion [20]. lncRNA PVT1 acts as a ceRNA to adsorb miR-143 and upregulates the expression of HK2 to encourage the proliferation and metastasis of GBC cells [21]. However, the lncRNA-associated ceRNA network in LSCC remains unclear. So, this study aims to build a ceRNA regulatory network for a better understanding of lncRNA's molecular mechanism in LSCC.
In this study, the LSCC transcriptome data was downloaded from the TCGA database, including the RNA expression profile of 111 tumor samples and 12 normal samples. Using bioinformatics tools, we constructed a ceRNA network containing 61 lncRNAs, 21 miRNAs, and 77 mRNAs. Meanwhile, a four-lncRNA prognosis model based on the ceRNA network was also established. Furthermore, we found that the four lncRNAs had apparent influences on the proliferation, migration, and invasion of LSCC cells in vitro. The above analysis and experiment results show that the four lncRNAs may serve prognosis biomarkers and become Therapeutic targets of LSCC in the future.

Screening of LSCC prognostic signatures
Survival analysis and univariate cox regression analysis were performed to explore the correlation between the lncRNAs from the ceRNA network and OS of LSCC patients by the survival R package [30]. Then prognosis model was built by multivariate regression analysis. And the following formula was used to calculate the prognostic risk score of LSCC patients: Risk score = βlncRNA1 * explncRNA1 + βlncRNA2 * explncRNA2 + … + βlncRNAn * explncRNAn ('β' is the regression coefficient of lncRNAs and 'exp' is the expression of corresponding lncRNAs). Through the survi-valROC R package, a receiver operating characteristic (ROC) curve was drawn to assess the prognosis model's accuracy.  The cells were transfected with prognosis related-lncRNA using Lipofectamine™ 3000 (Thermo Fisher Scientific, USA).

Quantitative real-time PCR
According to the manufacturer's protocol, the total RNA was extracted from cells by TRIzol reagent (Invitrogen, CA, USA). qRT-PCR analysis on lncRNA was performed using HiScript® II Q Select RT Super-Mix for qPCR Kit (vazyme, Nanjing, China) and ChamQ Universal SYBR qPCR Master Mix (vazyme, Nanjing, China). 18S was used as the endogenous control.

CCK-8 assay
Cell proliferation was examined by Cell Counting Kit-8 (CCK-8; GLPBIO, USA). The cells were incubated in the 96-well plates (1 × 10 3 cells per well) for 24 h, and then 10 μl CCK-8 reagents were added to each well. Cell viability was determined by detecting the absorbance at 450 nm.

Cell migration and invasion assay
Transwell chambers with 8.0 μM pore polycarbonate membrane insert (Corning, USA) assessed cell migration and invasion abilities. For cell invasion assay, 40 μL matrigel solution (Matrigel: medium = 1: 4) was added to transwell inserts and solidified for 3 h at 37°C. Then 500 μL DMEM with 10% FBS was added Fig. 4 Kaplan-Meier survival curves of six lncRNAs (log-rank p < 0.05). One hundred eleven LSCC patients were divided into two groups (high expression and low expression) according to the median of relative expression value to the lower chamber. Cells were resuspended with serum-free medium and plated into transwell inserts at 1 × 10 5 cells/well. The cells on the filter's upper surface were removed after they were cultured at 37°C for 24-72 h. For cell migration assay, the matrigel solution was not needed on the inserts, and other steps were as same as invasion assay. Cells on the lower surface of the filter were fixed with 4% paraformaldehyde for 20 mins and then washed twice with PBS before they were stained with 1% crystal violet solution for 10 mins. The stained cells were counted by a microscope.

Statistical analysis
The statistical analysis software included R 3.4.3 and SPASS 22.0. All measurement data were presented as mean ± SD. Two-tailed Student's t-test was used to compare two groups. P-value < 0.05 was considered statistically significant.

DEmRNAs functional enrichment analysis
For exploring the function and molecular mechanism of DEmRNAs in the occurrence and development of laryngeal squamous cell carcinoma, we performed GO and KEGG enrichment analysis on DEmRNAs (Supplementary Material 2). The top 20 KEGG pathways and top 20 GO biological process terms (p < 0.05), sorted by the p-values, were chosen for Bar graph and bubble graph ( Fig. 2A, B). Among these pathways and biological process terms, some terms containing ECM-receptor interactions, Tight junction, and PPAR signaling pathway are correlated to the proliferation and metastasis of tumors. In

Construction of the lncRNA-miRNA-mRNA ceRNA network
To explore the ceRNA Molecular regulatory network centered on lncRNA in LSCC, we constructed the ceRNA network according to the workflow (Fig. 3A). First, the interactions between lncRNAs and miRNAs were predicted by the starbase database among the DElncRNAs and DEmiRNAs. We obtained 168 lncRNA-miRNA pairs (Supplementary Material 4). Next, the interactions between miRNAs and mRNA were screened by TargetScan, miRDB, and miRWalk among the DEmiRNAs and DEmRNAs. We obtained 118 miRNA-mRNA pairs that all interactions exist in three databases (Supplementary Material 4). Then we used Cytoscape software to construct a lncRNA-miRNA-mRNA ceRNA network consisting of 61 lncRNAs, 21 miRNAs, and 77 mRNAs with a total of 179 interactions, based on the lncRNA-miRNA pairs and the miRNA-mRNA pairs (Fig. 3B).

Screening lncRNA-associated prognosis factors based on ceRNA network
Kaplan-Meier analysis results showed that six lncRNAs, including AC011933.2, LINC00689, LINC02570, FAM225B, LINC02086, and LINC02575 were closely   (Fig. 5) was a forest plot for the four lncRNAs consisted of the prognosis model. The survival risk score of all LSCC patients was presented in (Fig. 6A, C). (Fig. 6B) was a survival state graph for all LSCC patients. For detecting the specificity and sensitivity of the model. The Kaplan-Meier analysis result showed that the high-risk group's survival rate tended to lower than the low-risk group (Fig. 7A; p < 0.05). The ROC analysis results showed that the area under curve (AUC) value was 0.609 in the 1st year, 0.793 in the 3rd year, 0.752 in the 5th year, 0.876 in the 7th year, 0.915 in the 9th year, which indicated that the cox risk model has high specificity and sensitivity (Fig. 7B-F).

Prognosis-related lncRNAs affected the proliferation, migration, and invasion of LSCC cells
To understand the role of prognosis-related lncRNAs in the oncogenesis of LSCC, we elevated the expression level of prognosis-related lncRNAs by transfected pcDNA3.1(+)-prognosis-related lncRNAs into LSCC cells. The results showed that all four lncRNAs expression increased in the experimental groups compared to the control groups (Fig. 8A, B). The CCK-8 assay demonstrated that both LINC02086 and LINC02575 promoted Hep-2 and TU 177 proliferation while both AC011933.2 and FAM30A inhibited Hep-2 and TU 177 proliferation (Fig. 8C, D; p < 0.05). Then, the migration and invasion were detected by transwell assay. The results revealed that both LINC02086 and LINC02575 promoted Hep-2 and TU 177 cell migration and invasion while both AC011933.2 and FAM30A inhibited Hep-2 and TU 177 migration and invasion (Fig. 8E, F; p < 0.05).

LINC02086 promotes migration by miR-770-5p/SLC26A2 axis in LSCC
In order to validate the ceRNA model, LINC02086/miR-770-5p/SLC26A2 axis was extracted from the ceRNA network for experimental verification. Firstly, SLC26A2 mRNA level was significantly increased after the Hep-2 and TU 177 cells were transfected with LINC02086 overexpression vector compared with the control group, while SLC26A2 mRNA level was decreased after the Hep-2 and TU 177 cells were transfected with LINC02086 and minic miR-770-5p (Fig. 9A-B). Dualluciferase report analysis showed that the activity of Dual-luciferase was significantly decreased in the wildtype LINC02086 and SLC26A2 groups after transfection with minic miR-770-5p (Fig. 9C-D). The migration experiments showed that miR-770-5p could partially reverse the effect of LINC02086 on the migration of larynx carcinoma cells, while SLC26A2 could weaken the effect of miR-770-5p on the migration of larynx carcinoma cells (Fig. 9E). These results suggested that LINC02086 promoted the up-regulation of SLC26A2 by adsorbing miR-770-5p, thereby promoting the migration of larynx carcinoma cells.

Discussion
LSCC, a type of common head and neck squamous cell carcinomas, causes a severe threat to peoples' health all over the world [31]. As the incidence of laryngeal squamous cell carcinoma is increasing year by year, it is urgent to clarify the pathogenesis and identify significant prognosis biomarkers for improving the current treatment to LSCC. Our study demonstrated that a four-lncRNA prognosis model based on the ceRNA network had a sufficient ability to predict the prognosis of LSCC. Moreover, the four lncRNAs played a significant role in the proliferation and metastasis of LSCC. Therefore, these prognosis-related lncRNAs may serve as a new therapeutic target in the future.
LncRNA plays a vital role in the occurrence and development of tumors. Recently, more and more studies have shown that the ceRNA molecular regulation model is a common way for lncRNA to promote or inhibit tumor growth [32]. For example, ABHD11-AS1 serves as a competitive endogenous RNA to upregulate STAT3 by sponging miR-1301-3p in PTC [33]. Linc-DYNC2H1-4 promotes EMT and CSC phenotypes by acting as a sponge of miR-145 in pancreatic cancer cells [34]. However, most previous studies are limited to a single molecular regulation axis, and there is a lack of comprehensive and in-depth research on molecular regulation networks. As we know, the interaction between lncRNA and miRNA is not the only correspondence relationship. One lncRNA may have many different MREs, and other lncRNAs may also have the same MREs. For example, lnc HOTAIR can act as a molecular sponge for both miR-206 and miR-124 [15,35]. Both lncRNA TUG1 and lncRNA PVT1 can interact with miR-145 [36,37]. Similarly, the interaction between miRNA and mRNA is not the only correspondence relationship. Therefore, it is helpful for a comprehensive understanding of the molecular mechanism of lncRNA to elucidate the regulatory network of lncRNA-miRNA-mRNA.
Our study determined 1961 lncRNAs, 69 miRNAs, and 2224 mRNAs with abnormal expression genes in LSCC according to the transcriptome sequencing data derived from the TCGA. Then, functional enrichment analysis, including GO, KEGG, and GSEA, were performed by bioinformatics tools. We found that DEmRNAs were enriched in some tumor-related biological processes or pathways such as ECM-receptor interactions, Tight junction, PPAR signaling pathway, Human papillomavirus infection, PI3K-Akt signaling pathway. Next, a lncRNA-miRNA-mRNA ceRNA network was constructed according to DElncRNAs, DEmRNAs, and DEmiRNAs, and it included 61 lncRNAs, 21 miRNAs, and 77 mRNAs. Furthermore, we extracted a four-lncRNA (AC011933.2, FAM30A, LINC02086, LINC02575) prognosis model from the ceRNA network. These results indicate that the ceRNA network is useful for uncovering the molecular mechanism of lncRNA in LSCC and may become prognosis biomarkers of LSCC.
Several risk score systems were established to predict LSCC patients' prognosis in previous studies [38,39]. However, all these prognosis biomarkers based on the predicting model have not been verified by experiments in vitro. Our study constructed a four-lncRNA prognosis model and confirmed the biological functions of the four lncRNAs in LSCC.
The four lncRNAs' functions in the progress of LSCC remain unclear because of the few associated studies. In our research, we found that the expression of AC011933.2 and FAM30A is abnormal in LSCC, and they have a negative correlation with the OS of LSCC patients. On the contrary, LINC02086 and LINC02575 are positively correlated with the overall survival of LSCC patients. Furthermore, the four lncRNAs were found to promote or inhibit LSCC cell proliferation, migration, and invasion in vitro. Therefore, the four lncRNA may be considered as new tumor suppressor genes and oncogenes that deserve further research.
Inevitably, our research has several innate limitations that need to be addressed. Although the bioinformatics process of constructing the ceRNA network was designed reasonably, and experiments were performed in vitro to verify the functions of prognosis-associated lncRNAs, our study's main disadvantage was the lack of validation of the lncRNA-miRNA-mRNA molecular regulatory axis extracted from the ceRNA network. Besides, our cox risk model was based on data derived from the TCGA database, and clinical validation of larger samples is necessary. Despite the above drawbacks, our study results show that our cox risk model with 4 lncRNAs play a significant role in the occurrence and development of LSCC.

Conclusion
In conclusion, we have identified a lncRNA-related prognosis model that can effectively predict OS in LSCC. Moreover, we verified these lncRNAs could influence the biological function of LSCC cells.