LSCC tissues and paired ANM tissues were collected from patients who had undergone surgery at the Department of Otolaryngology Head and Neck Surgery, First Hospital of Shanxi Medical University. None of the patients received chemotherapy or radiotherapy before surgery. The clinical samples were collected from LSCC patients after obtaining their informed consent. Studies were performed with the approval of the medical ethics committee of First Hospital of Shanxi Medical University.

Total RNA was isolated using TRIzol reagent (Ambion, USA) following the manufacturer’s protocol. After evaluation of purity, concentration, and integrity, 3 μg of RNA per sample was used in library preparation for RNA sequencing. The rRNA was removed using the Ribo-zero rRNA Removal Kit (Epicentre, USA), and the rRNA-depleted sample was subsequently used to generate sequencing libraries using the NEB Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s protocols. In brief, the rRNA-depleted RNA samples were subjected to RNA fragmentation, first strand cDNA synthesis with random hexamer primers, second strand cDNA synthesis with dUTP, adenylation of 3’ ends, barcoded adapter ligation, purification, and amplification using PCR to generate sequencing libraries. After cluster generation, libraries were sequenced on an Illumina HiSeq 4000 platform by Novogene (Beijing, China).

Clean reads were obtained by removing reads adapter, reads containing poly-N, and low-quality reads from raw data. The reads mapping and transcript assembly was done using TopHat (v2.0.9) and Cufflinks (v2.1.1) based on the human reference genome (GRCh37) from the Ensembl database (https://www.ensembl.org/). Then, the coding potential of transcripts was predicted using CPC2 (beta) and CNCI (v2). Finally, the read counts in fragments per kilobase million (FPKM) of transcripts were used as input data for a more in-depth bioinformatic analysis. The TCGA expression data of LSCC were obtained from the TCGA data portal (https://cancergenome.nih.gov/). The associated clinical features were obtained from the University of California, Santa Cruz Xena (UCSC Xena; https://xenabrowser.net/) database. The significantly differentially expressed lncRNAs in the Shanxi and TCGA cohorts were screened using the R package DEseq2 (v1.26.0) (25) with the cut-off p-value< 0.05 and | log₂(Fold change) | ≥ 1.

Complete-linkage clustering, which is one of the hierarchical clustering algorithms, was used to analyze the differential expression of lncRNAs in relation to various clinical features. The heat map to visualize the expression data was constructed using the R package pheatmap (v1.0.12, https://CRAN.R-project.org/package=pheatmap). Principal component analysis (PCA) was used to reduce the dimensionality of high dimensional expression datasets and evaluate the degree of difference between LSCC and normal tissue. The 3D visualization of PCA was realized using the R package pca3d (v0.10.1, https://CRAN.R-project.org/package=pca3d).

Functional annotation of dysregulated lncRNAs in LSCC was performed using the GSEA (v3.0) software with Molecular Signatures Database (v7.0, https://www.gsea-msigdb.org/gsea/msigdb/). Significantly upregulated and downregulated lncRNAs in both the Shanxi and TCGA cohorts were obtained by Venn analysis. Pearson correlation analysis with Student’s t-test was used to obtain protein-coding RNAs that are correlated with lncRNAs based on their expression levels, with a cut-off set at | Pearson’s r | ≥0.4 and p-value< 0.05. Based on the ranking of Pearson’s r, the five most positively and five most negatively correlated protein-coding RNAs to each of the lncRNAs were selected as the input data for GSEA. Mutual functional pathways in two LSCC cohorts were mined using a cut-off set at FDR q-value ≤ 0.01 after GSEA.

The previously identified protein-coding RNAs that were correlated with the expression of the seven lncRNAs were used to determine the different associated biological signaling pathways using Metascape (https://metascape.org/).

Univariate Cox proportional hazards regression analysis was used to investigate the relationship between the continuous expression levels of each significantly dysregulated lncRNA and the overall survival (OS) of patients in the TCGA cohort. A total of 111 LSCC patients’ data from TCGA database were used for survival analysis, in which 50 patients were reached the survival endpoint. Based on the expression fold-change of each lncRNA, the output of the Cox regression analysis was divided into an upregulation set and a downregulation set. Under the criterion of p-value< 0.05, 20 lncRNAs in the upregulation set and 4 lncRNAs in the downregulation set were defined as significant high-risk lncRNAs for LSCC patients. The least absolute shrinkage and selection operator (LASSO) algorithm with a penalty parameter of 1000 and a 10-fold cross-validation was used to shrink sparse high-dimensional data and develop a prognostic prediction model from the output of the Cox regression analysis. The LASSO algorithm was implemented using the R package glmnet (v2.0_16, https://CRAN.R-project.org/package=glmnet). To determine the risk levels of the prognostic prediction model for OS, the risk score for each LSCC patient was calculated using the formula: risk score =  ∑ i = 1 n ( c lncRNAi   ×  e lncRNAi ) ; where c l n c R N A i is the relation coefficient of the risk factor lncRNAi; e l n c R N A i is the expressional value of the risk factor lncRNAi; and n is the number of lncRNAs included in the prognostic prediction model. The LSCC patients were classified into high or low risk groups according to the median value of the computed risk scores. The Kaplan-Meier (K-M) estimator model with the log-rank test was used to estimate the OS, disease-specific survival (DSS), and progression-free interval (PFI) for LSCC patients between the high and low risk groups. The survival curve was realized using the R packages survival (v2.44-1.1, https://cran.r-project.org/package=survival) and survminer (v0.4.4, https://CRAN.R-project.org/package=survminer). The sensitivity and specificity of the lncRNAs prognostic prediction model were analyzed using a time-dependent receiver operating characteristic (ROC) curve. The area under the curve (AUC) was compared among the 1-, 3- and 5-year survival rates. The calculation and visualization of the ROC curves were done using the R package timeROC (v0.3, https://CRAN.R-project.org/package=timeROC). Moreover, a total of 86 LSCC samples with follow-up data was used to validate the efficiency of prognostic model. The expression levels of 7 lncRNAs in the 86 LSCC samples were determined by qPCR analysis. A nomogram for predicting OS was built using the R package rms (v5.1-3, https://CRAN.R-project.org/package=rms). The survival probabilities were predicted for the 1-, 3- and 5-year marks. The validation of the nomogram prediction model was accessed using bootstrapped calibration curves.

The correlation among seven lncRNAs was analyzed using the R package stats (v4.1.0) and visualized using the R package corrplot (v 0.92, https://cran.r-project.org/package=corrplot). An lncRNAs-miRNAs-mRNAs network was constructed based on competitive endogenous RNA (ceRNAs). For target miRNA prediction of seven lncRNAs, miRanda (v3.3a, https://anaconda.org/bioconda/miranda), LncTar (v1.0, http://www.cuilab.cn/lnctar) and RNA22 (v2.0, https://cm.jefferson.edu/rna22/) were used to independently predict the binding between lncRNAs and miRNAs. Next, the top 30 miRNAs with strongest binding ability to each lncRNA in three softwares were intersected to obtained potential target miRNAs. For target mRNA prediction of miRNAs, Targetscan (v7.1, https://www.targetscan.org/), miRdb (https://mirdb.org/mirdb/) and mirWalk (http://mirwalk.umm.uni-heidelberg.de/) were utilized. Target Score >=20 in miRDB, energy<= -5 in miRWalk and Total context++ score<= -0.1 in Targetscan were used as the threshold, and target mRNAs for each miRNA were obtained from intersection of three databases. Finally, Cytoscape (v3.7.1. https://cytoscape.org/) was used to construct the lncRNAs-miRNAs-mRNAs network.

Human LSCC TU-177 cell lines (Bioleaf Biotech, China) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM; Life Technologies, USA) supplemented with 10% fetal bovine serum (FBS; BI, USA). Human LSCC FD-LSC-1 (26) cell lines were cultured with Bronchial Epithelial Cell Growth Medium (BEGM; Lonza, Switzerland) supplemented with 10% FBS. All cells were maintained in a humid atmosphere of 5% CO₂ at 37°C. The siRNAs used to target lncRNAs were synthesized by GenePharma Co., Ltd (Shanghai China). The siRNA sequences are listed in Table S1 . Cells (2.0 × 10⁵ cells per well) were seeded in 6-well plates 24 h before transfection. The siRNAs were transfected using Lipofectamine 3000 (ThermoFisher Scientific, USA) according to the manufacturer’s protocol. A qPCR analysis was performed to determine the relative expression levels of the target lncRNAs. Expression levels of 18S rRNA were used as internal controls. The primer sequences used are listed in Table S2 .

Cell proliferation was assessed using the EdU staining kit (RiboBio, Guangzhou, China) according to the manufacturer’s protocol. Trans-well migration and invasion assays were performed as previously described with modification (27). Briefly, for the migration assay, cells were seeded in 24-well trans-well plates (Falcon, USA) with a density of 1 × 10⁶ cells per well in 100 µL of serum-free media. Meanwhile, for the invasion assay, 1 × 10⁵ cells were seeded in each well of the 24-well trans-well plates coated with Matrigel (Corning, USA). After 24 h of incubation, the cells in the upper chamber of the membranes were removed using cotton swabs, then cells in the lower chamber were stained with crystal violet (Amresco, USA) for 10 min and washed with PBS.

A total of 107 LSCC and paired ANM tissues were collected from 107 patients who were clinically diagnosed with LSCC to constitute the Shanxi cohort. The RNA samples from the Shanxi cohort were subjected to transcriptome sequencing ( Figure 1A ). Furthermore, the RNA-seq data of 111 LSCC patients were obtained from the TCGA database. The analysis workflow used in this study is shown in Figure S1 . The PCA results showed that the lncRNA expression patterns in LSCC tissues are different from that of the paired ANM tissues ( Figure 1B ). This trend was further confirmed by the TCGA cohort ( Figure 1C ). Next, an unsupervised hierarchical clustering based on the lncRNA expression data and clinical features was performed. The results showed that numerous lncRNAs were differentially expressed between the LSCC and ANM/normal groups ( Figures 1D, E ).

Based on a threshold of p< 0.05 and |fold-change| ≥ 2, a total of 6,060 differentially expressed lncRNAs were identified in the Shanxi cohort. Out of these, 2,949 were upregulated and 3,111 were downregulated in LSCC tissues compared with ANM tissues ( Figure 2A and Table S3 ). Meanwhile, in the TCGA cohort, 4,381 upregulated and 1,524 downregulated lncRNAs were identified in LSCC samples compared to normal samples ( Figure 2B and Table S4 ). Using Venn analysis, we obtained 594 upregulated and 197 downregulated lncRNAs from both cohorts ( Figure 2C ), which were named as significantly dysregulated lncRNAs (SDLs). Next, we predicted the biological processes and pathways associated with the SDLs using GSEA. Results showed 36 mutual biological processes and pathways (FDR q-value ≤ 0.01) in both cohorts, including terms such as epithelial to mesenchymal transition (EMT), nasopharyngeal carcinoma, cell cycle checkpoints, tumorigenesis, thyroid carcinoma poorly differentiated, and malignant mesothelioma ( Figure 2D ). These suggested that the SDLs were involved in the regulation of multiple tumorigenesis-related processes and pathways.

EMT plays an important role in the development and metastasis of cancer (28, 29). We observed that the SDLs in both LSCC cohorts were significantly enriched in EMT ( Figures 2E, F ), suggesting that the dysregulated lncRNAs may affect LSCC invasion and metastasis by regulating the EMT process.

To screen for lncRNAs that are significantly associated with the OS of LSCC patients, all SDLs and prognosis-related clinical information of TCGA patients were subjected to Cox proportional hazards regression analysis. Using the criterion of p-value ≤ 0.05, 20 upregulated (TEX41, LINC00392, LINC02575, ENSG00000233397 (named as ENSG233397 thereafter), BARX1-DT, TSPEAR-AS1, ENSG00000236841, LSAMP-AS1, HOXB-AS4, MNX1-AS1, LINC01391, LINC02086, LINC02315, LINC01385, ENSG00000253477, CABP1-DT, ENSG00000260978, ENSG00000267919, LINC02893 (also known as AL513318.2), and ENSG00000272763) and 4 downregulated (LINC00278, HOXA-AS2, LINC02568, and LINC02188) lncRNAs were obtained from the univariate Cox analysis ( Figures 3A, B ). The association between the expression levels of these 24 lncRNAs and OS of LSCC patients was verified using K-M survival analysis ( Figures S2, S3 ). Results showed that LSCC patients with high expression levels (≥ median FPKM) of each of the 20 upregulated lncRNAs had a shorter median OS time than that of the low expression group (< median FPKM). In contrast, LSCC patients with high expression levels of each of the 4 downregulated lncRNAs had a longer median OS time than that of the low expression group ( Figure 3C ).

To construct a model for predicting the prognosis of LSCC patients, LASSO regression analysis was used to narrow down the candidate lncRNAs based on the output of the Cox analysis. This led to the identification of seven lncRNAs including ENSG233397, BARX1-DT, LSAMP-AS1, HOXB-AS4, MNX1-AS1, LINC01385, and LINC02893 ( Figure 4A ). These seven lncRNAs were significantly upregulated in LSCC tissues compared to normal counterparts in both the Shanxi and TCGA cohorts ( Figures S4, S5 ). Combining the relative expression of the lncRNAs and their corresponding coefficients, the risk score for the seven-lncRNAs of each patient was calculated as follows: risk score = 0.1331 × e ENSG 233397 + 0.0940 × e BARX 1 - DT + 0.1254 × e LSAMP - AS 1 + 0.1400 × e HOXB - AS 4 + 0.1560 × e MNX 1 - AS 1 + 0.1611 × e LINC 01385 + 0.0959 × e LINC 02893 . The risk scores of the seven-lncRNAs panel were significantly associated with poor OS (HR:6.21 [3.27-11.81], p<0.0001), DSS (HR:4.34 [1.83-10.26], p=0.0008), and PFI (HR:3.78 [1.92-7.43], p=0.0001) in TCGA cohort. Based on the median value of the risk scores, all patients in the TCGA cohort were divided into high and low risk groups. K-M survival analysis revealed that LSCC patients in the high risk group had shorter median OS (23.17 months), DSS (48.87 months), and PFI (19.23 months) than those in the low risk group (OS: 91.10 months, DSS: 213.90 months, and PFI: 213.90 months; Figures 4B , S6A, C ).

Time-dependent ROC curves were drawn to assess the prognostic ability of the seven-lncRNAs risk scores. In the time frame of 1-, 3-, and 5-years, the AUC index of the seven-lncRNAs risk score for OS were 65.77%, 78.97%, and 80.76%; for DSS were 58.47%, 75.99%, and 76.55%; and for PFI were 63.20%, 73.41%, and 74.70%, respectively ( Figures 4C , S6B, D ). To validate the seven-lncRNAs prognostic model, the expression levels of these lncRNAs in 86 LSCC samples were determined by qPCR, then the association between seven-lncRNAs risk score and overall survival of the 86 LSCC patients was evaluated by Kaplan-Meier analysis. As expected, results further confirmed that the seven-lncRNAs model was able to predict prognosis of LSCC patients ( Figure 4D ). Furthermore, time-dependent ROC curves analysis showed that the AUC index of the seven-lncRNAs risk score of the 86 LSCC patients were 54.90%, 62.90%, and 69.10% in the time frame of 1-, 3-, and 5-years for OS ( Figure 4E ).

We also noted that LSCC patients in the low risk group had a lower mortality rate (20/55 patients, 36.36%) than those in the high risk group (30/56 patients, 53.57%; Figures 5A, B ). Furthermore, the identified seven lncRNAs had higher expression levels in patients of the high risk group than those in the low risk group ( Figure 5C ). Further analysis of the seven-lncRNAs data of the 86 LSCC patients showed that patients in high risk group had a higher mortality rate (13/43 patients, 30.23%) than those in the low risk group (6/43 patients, 13.95%), and patients in high risk group had higher expression levels of seven lncRNAs ( Figures 5D–F ).

To analyze the interconnections among seven lncRNAs, a total of 18 miRNAs was predicted based on ceRNA hypothesis ( Figure S7A ). LINC01385 and BARX1-DT had a common targeted miRNA named as miR-4739 ( Figure S7B ). ceRNA network analysis revealed that direct or indirect regulatory relationships were existed among these seven lncRNAs ( Figure S7C ). The Pearson’s correlation analysis of the seven lncRNAs versus each other in the Shanxi and TCGA cohorts is seen in Figure 6A . LNM is a predominant independent factor associated with the mortality rate in LSCC (30). We found that the seven lncRNAs had relatively higher expression levels in patients with LNM than those of non-metastatic patients in both cohorts ( Figures 6B, C , S8 ). To further improve the predictive accuracy of our model, a nomogram for OS prediction was constructed by integrating the risk score and clinical features of LNM. As shown in Figure 6D , the upward vertical line can be respectively drawn from the line of the risk score and LNM to the line where their respective points can be calculated. Based on the sum of points of the seven-lncRNAs risk score and LNM, a downward vertical line can be drawn from the line of the total points to the line of the 1-, 3-, or 5-year OS probability to respectively calculate the OS probability. Additionally, the calibration curves for the prognostic nomogram estimating the probability of 1-, 3-, and 5-years OS performed well with the ideal model ( Figure 6E ). Decision curve analysis of the prognostic nomogram for 1-, 3-, and 5-years OS risk showed that the model has a net benefit of risk stratification in the gradient prediction time ( Figure S9 ). Collectively, these findings demonstrated that a seven-lncRNAs panel for predicting the prognosis of LSCC patients was constructed using an integrated analysis of our large-scale RNA-seq data and the TCGA data.

We sought to determine whether the lncRNAs are involved in regulating LSCC progression. We used the protein-coding genes whose expression levels significantly correlated with the seven lncRNAs for functional annotation through Metascape ( Table S5 ). Results showed that the seven lncRNAs were mainly involved in the Wnt signaling pathway and pathways regulating the pluripotency of stem cells ( Figure 7A ). Furthermore, protein-protein interaction enrichment analysis revealed that the seven lncRNAs were also involved in cell cycle checkpoints and the JAK-STAT signaling pathway ( Figures 7B, C ). Since these pathways are known to contribute to cellular growth and differentiation, we hypothesized that these lncRNAs may regulate the progression of LSCC.

To investigate the functional roles of each of the seven lncRNAs in our model, we performed siRNA-mediated knockdown targeting the identified lncRNAs in LSCC cell models. Our qPCR analysis confirmed that the expression levels of the seven lncRNAs decreased significantly in the LSCC cell lines FD-LSC-1 and TU-177 after siRNAs transfection ( Figure 8 ). The EdU staining results showed that knocking down ENSG233397, BARX1-DT, LSAMP-AS1, HOXB-AS4, MNX1-AS1, LINC01385, and LINC02893, separately, inhibits the proliferation of LSCC cells to different degrees ( Figures 9A , S10A ). In addition, we found that knockdown of each of the seven lncRNAs resulted in significantly decreased migration and invasion abilities of LSCC cells ( Figures 9B, C ; S10B, C ). Collectively, our findings suggest that the lncRNAs used in our model promote the proliferation, migration, and invasion abilities of LSCC cells in vitro.