The ESCC RNA expression profiles and corresponding clinical information were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database including GSE53624 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53624), GSE53622 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53622), and GSE53625 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53625) datasets. Samples with complete survival information are retained, while those patients without survival information are eliminated. To develop prognostic prediction lncRNA models, ESCC samples from GSE53624 were treated as a training set. GSE53622 and GSE53625 sets were test and validation datasets. The aforementioned datasets were generated with Agilent-038314 (GPL18109). Through re-annotating microarray probes (see details in the Supplementary Material) (Harrow et al., 2012; White et al., 2014; Guo et al., 2018), we gained the expression values of lncRNAs from ESCC cohorts (Supplementary Table S1). Probes with missing expression values in more than 20% of patients were discarded.

To single out those lncRNAs which were significantly associated with the prognosis of ESCC patients, both univariable Cox regression and KM survival analysis (the median lncRNA expression value as the cutoff value) were used in the training dataset. Those with Cox p < 0.05 and log rank p < 0.05 were considered OS-associated candidates. The LASSO regression method was then applied to obtain the strongest survival-related lncRNAs in the training set. Subsequently, the selected prognostic lncRNAs by KM, Cox, and LASSO regression were performed to develop combination models for estimating the ESCC prognosis risk as follows: risk score (RS) = ∑ Ni = 1 (Exp * coefficient), where N is the number of selected lncRNAs, Exp is the corresponding lncRNAs’ expression level, and the coefficient is calculated by the univariable Cox analysis. Based on the above formula, the RS of each combination model for each ESCC patient was calculated and ROC curve analysis was applied to make comparison of the survival prediction ability among those constructed multi-lncRNA signatures in the training set.

Human ESCC cell lines KYSE410 and TE5 were cultured in RPMI 1640 (Gibco) medium with 10% fetal bovine serum (TransSerum) and 1% streptomycin–penicillin solution (Gibco). All cells were cultured in a 5% CO₂ constant temperature incubator. Small interfering RNAs (siRNAs) targeting LINC01273 (siLINC01273-1: 5′-GAC​ACA​GAA​GGA​CAA​UGU​UTT-3′; siLINC01273-2: 5′-GAC​ACA​AAG​UGA​CAG​AAU​GTT-3′) were synthesized by GenePharma Co. (Suzhou, China). Following the instructions, siLINC01273 was transfected at a concentration of 40 nM using Lipofectamine RNAiMAX Transfection Reagent (Invitrogen) with Opti-MEM (Gibco). After transfection for 48 h, the RNAs were harvested.

Total RNA was reverse transcribed into cDNA by HiScript Q RT SuperMix for qPCR (Vazyme) after extracting by RNA-easy Isolation Reagent (Vazyme). The real-time RT-qPCR assay was conducted with an ABI 7500 system (Corbett Life Science) using ChamQ SYBR Color qPCR Master Mix (Vazyme) with the guide of its manufacturer’s instructions. The primers for RT-qPCR of LINC01273 were 5′-TGT​TGC​GGT​GTT​CAG​GGG​TTT-3′ (forward) and 5′-GTC​TGG​CTT​CTT​TCA​CTG​AGC-3′ (reverse). The primers for beta-actin were 5′-CAA​CTG​GGA​CGA​CAT​GGA​GAA​A-3′ (forward) and 5′-GAT​AGC​AAC​GTA​CAT​GGC​TGG​G-3′ (reverse). The relative mRNA expression was normalized to beta-actin as reference.

For the MTS assay, after transfection for 36 h, 5,000 cells/well were seeded into 96-well plates. After adding MTS solution (Promega) and incubating for 2 h, the absorbance was recorded at 490 nm using an ELISA plate reader. For the colony formation assay, 500 cells/well were planted in 12-well plates and continuously grown for 2 weeks until a single colony was formed. After fixing with methanol, these colonies were stained with 0.1% crystal violet.

ESCC cells were transfected with siRNAs for 36 h, and then serum starvation was performed for 12 h. For invasion assays, upper transwell chambers (Falcon) should be pre-coated with Matrigel (BD Biosciences) and then left in the incubator for 1 h. 5×10⁴ cells in 200 μL serum-free cell suspensions were added in the upper transwell chambers, while 500 μL medium containing 10% FBS was added in the bottom chamber. 36 h later, pictures were taken with a microscope magnifying ×200 after fixing and then staining the migrated or invasive cells from upper chambers.

The 50th percentile of the risk score is defined as the threshold to classify the high-risk group and the low-risk group. KM analysis was applied to evaluate and validate the survival prediction performance of the lncRNA signature in different ESCC cohorts. The time-dependent ROC curve was used to compare the prediction ability of the lncRNA signature with that of other clinical features at different survival times. And univariable and multivariable Cox regression and stratification analysis were used to test whether the multi-lncRNA risk score model was independent of other clinical characters. The R program (3.5.1) including R packages named survival, survminer, glmnet, pROC, and timeROC was used to perform the above analyses.

To explore the potential biological functions of lncRNAs, the Pearson correlation test was used to construct co-expressed networks of lncRNAs and the protein-coding genes (PCGs) in the GSE53625 dataset, and the PCGs that were highly correlated with lncRNAs (correlation coefficient >0.60/< -0.6, p < 0.001) were selected for GO and KEGG pathway enrichment analysis by the Cluego plugin in Cytoscape (Guo et al., 2018). SubpathwayMiner was also used to identify related pathways of the co-expressed PCGs in the KEGG database including entire pathways and sub-pathways.

All experiments were repeated for at least three times. The values are shown as mean ± SD. Prism 8 software was used to perform statistical analyses. Student’s t-test was employed for comparisons between two groups, and one-way ANOVA was performed for multiple-group comparisons. The differences with *p < 0.05, **p < 0.01, ***p < 0.001 were considered statistically significant.

There were a total of 179 ESCC samples used in this study, including 119 from GSE53624 and 60 from GSE53622, respectively. GSE53625 is the union of GSE53624 and GSE53622. The median survival age was 60 years. There were more male patients with ESCC than females (146 vs. 33), and most of the patients were dead (survival time, 3 days to 60 months). Other clinical characters are shown in Table 1. In addition, through re-annotating microarray probes, a total of 6,253 expressed lncRNAs and 17,434 expressed PCGs were obtained from GSE53624 and GSE53622.

ESCC samples from GSE53624 (n = 119) were treated as the training dataset to evaluate the relationship between ESCC OS and lncRNAs. After univariate Cox and KM analysis of lncRNAs’ expression level with clinical survival information, we identified a total of 209 lncRNAs (Figure 1A) related to ESCC patients’ OS significantly (Cox p < 0.05 and log rank p < 0.05), which could be used as prognostic candidates. Then, the LASSO regression algorithm via regression coefficient shrinkage based on a penalty that is proportional to size was utilized to screen out lncRNAs which were mostly correlated with ESCC survival among the 209-lncRNA set. As shown in Figure 1B, we found that the value of independent coefficients tended to zero with the increase of lambda value. Finally, we used threefold cross-validation and selected seven lncRNA candidates to construct the multi-lncRNA classifiers (Figure 1C).

To select a better predictive multi-lncRNA model with fewer lncRNAs, ROC curve analysis was performed to compare the prognostic prediction performance of the 2⁷-1 = 127 risk score combinations in the training dataset (Supplementary Table S2). All risk scores for each ESCC based on the corresponding lncRNA signature were calculated as the method described. Then, the six-lncRNA combination with the largest AUC value composed by AL445524.1, AC109439.2, LINC01273, AC015922.3, LINC00547, and PSPC1-AS2 was obtained (Figure 1D; Table 2). The RS of the six-lncRNA signature is as follows: RS = (-0.5460037×AL445524.1) + (-0.2473264× AC109439.2) + (0.4223392× LINC01273) + (-0.81843 × AC015922.3) + (0.7987309× LINC00547) + (0.8210199× PSPC1-AS2). The AUC of the six-lncRNA signature was 0.863 (95% CI: 0.798–0.928), higher than that of the seven-lncRNA model (0.855, 95% CI: 0.787–0.924, Figures 1E,F) and other lncRNA combinations. Therefore, we chose the six-lncRNA signature with fewer nodes and better survival prediction ability as the candidate classifier.

In the GSE53624 set, on the basis of the median risk score calculated by the six-lncRNA signature, patients were distinguished into two groups with different OS. Unfortunately, patients with ESCC from the high-risk group suffered a worst survival outcome than those from the low-risk group (log rank p < 0.001, Figure 2A). The five-year survival rate of patients in the low-risk group was 63.3%, which was significantly more than 15.25% of patients in the high-risk group.

For verifying the survival classification power of the lncRNA model, each patient from the validation GSE53622 set obtained their risk score values. Figure 2B shows the KM curves for patients with ESCC from the low/high-risk group in the GSE53622 dataset. We found that the median survival time in the high-risk group was 39.17 months less than 50.6 months in the low-risk group (five-year survival rate: 30% vs. 60%, log rank test p = 0.021). As for the entire dataset (GSE53625), patients with high risk scores suffered more undesirable outcomes than those with low risk scores (median survival time: 23.13 months vs. 51.3 months; log rank test p < 0.001, Figure 2C).

Moreover, Figure 3 shows the lncRNAs’ expression pattern of ESCC patients, the distribution of survival status, and their risk scores. For ESCC patients with high risk scores from the training set, the expression values of four lncRNAs (LINC01273, AC015922.3, LINC00547, PSPC1-AS2) were high, while the expression values of protective lncRNAs (AL445524.1, AC109439.2) were low. In contrast, the expression of prognostic lncRNAs showed the opposite pattern in patients with low risk scores in the training set (Figure 3A). Subsequently, we confirmed the similar survival distribution and risky or protective lncRNAs’ expression pattern in GSE53622 and GSE53625 sets (Figures 3B,C).

To evaluate the independence of the signature in survival prediction with other clinical characters including age, gender, and TNM stage, Cox regression analysis in GSE53624, GSE53622, and GSE53625 datasets was performed, and the multivariate Cox results of the multiple ESCC datasets showed that the six-lncRNA signature in OS prediction was independent of age and gender (high vs. low risk, HR = 4.97, p < 0.001, n = 119; HR = 2.26, p = 0.025, n = 60; HR = 2.11, p < 0.001, n = 179, Table 3). In addition, TNM stage affected the OS of patients with ESCC in GSE53624, GSE53622, and GSE53625 datasets (III vs. I + II: HR = 1.8, p < 0.001, n = 119; HR = 2.37, p = 0.009, n = 60; HR = 1.95, p < 0.001, n = 179, Table 3).

Time-dependent ROC curve analysis from 1 year to 5 years was applied to compare the survival prediction ability of the lncRNA signature with that of tumor grade, TNM stage, T stage, and N stage in the entire ESCC group (GSE53625, n = 179). The AUC values showed the predictive ability of the lncRNA signature (AUC from 1 year to 5 years: 0.698–0.909) was better than that of TNM stage (AUC from 1 year to 5 years: 0.486–0.67) and other features, especially at 5 years (Figure 4A). And the AUC of the combined model was the largest one compared to that of TNM stage or signature alone (AUC = 0.712, 95% CI = 0.645–0.779, Figure 4B), which further suggested the signature has potential to become a novel prognostic biomarker.

To evaluate whether the signature can further subgroup ESCC patients at high (III)/low (I, II) TNM stage, we performed stratification analysis in the entire dataset (GSE53625, n = 179). According to the TNM stage information of all the 179 patients, we found 87 patients at TNM low stage and 92 at TNM high stage. For patients at low TNM stage, the six-lncRNA signature could separate them into low- and high-risk groups with significantly different survival (five-year survival rate 59.1% vs. 18.6%, log rank test p < 0.001, Figure 4C). The signature can further classify patients at the high TNM stage into two groups with different prognostic outcomes (median survival: 28.7 months vs. 58.2 months; log rank test p < 0.001, Figure 4D). This result showed the potential ability of the six-lncRNA signature as a clinical auxiliary marker for TNM stage to subgroup patients with ESCC more accurately.

The Pearson test observed that the expression of 491 PCGs was significantly related to at least one of the six prognostic lncRNAs (coefficient >0.60/< −0.6, p < 0.001). GO and KEGG function analysis was then performed by Cluego and SubpathwayMiner. The results showed the 491 PCGs correlated with lncRNAs were significantly enriched in 37 GO terms and 36 KEGG pathways (p < 0.05, Supplementary Table S3). All these vital GO terms were organized into an interaction network based on similar functions in Cytoscape, and several clusters of functionally related GO terms were found such as ncRNA metabolic process, RNA process via interacting with those PCGs that affect cell cycle, regulation of actin cytoskeleton, MAPK signaling pathway, cell cycle, and TGF−beta signaling pathway (Supplementary Figure S1B).

We next investigated the biological roles of LINC01273 in maintaining the malignant phenotypes of ESCC cells. LINC01273 expression was examined in ESCC cell lines which our lab owned using qRT-PCR, and the results showed that LINC01273 was highly expressed in KYSE410 and TE5 cells (Figure 5A). Therefore, KYSE410 and TE5 cell lines were selected for further experiments. Firstly, we, respectively, transfected two individual siRNAs and confirmed LINC01273 was successfully knocked down by qRT-PCR (Figure 5B). We found that, by using the MTT assay and cell colony formation assay, silencing LINC01273 remarkably attenuated both the proliferation and colony formation capability of ESCC cells (Figures 5C,D). Transwell assays showed a significant suppression of the migration and invasive abilities of the two ESCC cell lines due to LINC01273 downregulation (Figures 5E,F). These results suggested that LINC01273 might enhance the ability of proliferation, migration, and invasion of KYSE410 and TE5 cells, demonstrating that LINC01273 may play oncogenic roles in ESCC.