Differentially expressed lncRNAs, miRNAs and mRNAs
To construct the metastasis-associated ceRNA network in lung cancer, we downloaded RNA expression profiles from TCGA database, which included 773 lung cancer samples. The differentially expressed (DE) mRNAs, lncRNAs and miRNAs between 32 metastatic tumor tissues and 741 nonmetastatic tumor tissues were explored using the “edgeR” package in R software. With the criteria of |logFC| > 1 and FDR < 0.05, we identified 1249 DE mRNAs (569 upregulated and 680 downregulated), 440 DE lncRNAs (221 upregulated and 219 downregulated) and 26 miRNAs (21 up- regulated and 5 downregulated) between nonmetastatic and metastatic tumor tissues. The results indicate that aberrant expression of these genes might be involved in lung cancer migration.
GO and KEGG pathway analysis of DE genes
To better understand the function of the identified DE genes, DE mRNAs were divided into upregulated and downregulated groups for analysis by DAVID and KOBAS bioinformatics resources. The top 15 GO biological processes and KEGG pathways of dysregulated genes are shown in Fig. 2 Among the top 15 GO and KEGG pathways of upregulated mRNAs (Fig. 2a, b), “GO:0006810~transport”, “GO:0003341~cilium movement”, “MAPK signaling pathway”, “ECM-receptor interaction” and “focal adhesion” are reported to promote invasion and metastasis in cancer [20, 21]. Downregulated genes participated in “GO:0031424~keratinization”, “GO:0008544~epidermis development” and “drug metabolism - cytochrome P450” (Fig. 2c, d), which reportedly inhibit invasion and metastasis of cancer [22, 23]. These results suggest that these above pathways may play important roles in lung cancer metastasis.
Construction of the ceRNA network
To construct the metastasis-associated ceRNA network, we predicted interactions among DE mRNAs, miRNAs and lncRNAs using bioinformatics tools. Target mRNAs of DE miRNAs were predicted using TargetScan and miRDB, and the intersection of these two databases was selected. After we discarded target mRNAs that were not among DE mRNAs, 117 mRNAs (95 upregulated and 22 downregulated) were used to construct the ceRNA network (Table S1). Next, we utilized DIANA to predict interactions between DE lncRNAs and DE miRNAs. We found the 22 miRNAs that are predicted to interact with 23 lncRNAs (Table S2). Ultimately, 23 lncRNAs (20 downregulated and 3 upregulated) (Table S3) and 22 miRNAs (4 downregulated and 18 upregulated) (Table S4) were included in the ceRNA network.
Based on the above findings (Table S1, 2), we constructed the ceRNA network using Cytoscape 3.6. Fig. 3 shows that 4 downregulated miRNAs, 22 upregulated mRNAs and 3 upregulated lncRNAs are involved in one ceRNA network. Meanwhile, 18 upregulated miRNAs, 95 downregulated mRNAs and 20 downregulated lncRNAs were involved in another ceRNA network.
Identification of a prognostic signature of lung cancer
Increasing studies have shown that lncRNAs effectively predict overall survival (OS) in cancer patients [24, 25]. To validate the association between the lncRNAs in the ceRNA network and OS of lung cancer, we randomly divided lung cancer patients into the training set (n = 387) and testing set (n = 386), and there was no significant difference in the clinical covariates between the two sets (P > 0.05) (Table S5). We analyzed the 23 DE lncRNAs in the training set by multivariate Cox regression analysis. The results demonstrated that a 6 lncRNA signal was a significant prognostic factor for lung cancer (Table S6). LINC01287, SNAP25-AS1 and LINC00470 were protective, whereas AC104809.2, LINC00645 and LINC01010 were associated with increased risk (Fig. 4a). Next, lung cancer patients in the training set were divided into high‐risk and low‐risk groups based on the following risk score formula: Risk Score= (0.0662) * LINC01287+ (0.0978) * SNAP25-AS1+ (0.0567) * LINC00470+ (-0.0745) * AC104809.2+ (-0.2018) * LINC00645+ (-0.0726) * LINC01010. Kaplan‐Meier analysis revealed that OS in the high-risk group was significantly lower than in the low-risk group (P=0.0016, Fig. 4b). The overall ten-year relative survival rates of the high-risk and the low-risk groups were 14.9% and 33%, respectively. Furthermore, the area under ROC curve (AUC) at 10 years for OS was 0.641 (Fig. 4c).
In order to validate the prognostic ability of the 6- lncRNA model, we were further performed Kaplan-Meier analysis and ROC curve analysis in the testing set. Lung cancer patients in the testing set were divided into high-risk and low-risk group (Fig. 5a) with statistically significant different overall survival (P=0.046, Fig. 5b) and ROC curve (AUC=0.61, Fig. 5c).
Next, all lung cancer patients were stratified by gender. The six-lncRNA signature could classify 486 male patients into high-risk group and low-risk group and there was the statistically significant difference between the high-risk and the low-risk group in Kaplan-Meier analysis (p = 0.042, Fig. 5d). Similarly, 287 female patients were divided into high-risk group and low-risk group with statistically significant different overall survival (p = 0.027, Fig. 5e). These results indicate that the 6 lncRNA signal is indeed an independent prognostic factor for predicting OS in lung cancer patients.
Clinical feature analysis of LINC01010
To further analyze whether these 6 lncRNAs affect clinical features in lung cancer from TCGA database, we divided lung cancer patients into high-expression and low-expression groups based on lncRNA expression. Kaplan-Meier survival curves revealed that only LINC01010 was positively correlated with OS (P = 0.02779) (Fig. 6a). OS at 5 and 10 years in the LINC01010 high-expression group was 49.9% and 26.8%, respectively, compared with 40.6% and 18.2% in the low-expression group. Therefore, LINC01010 was chosen for further study. Although there was no significant difference between normal tissues and lung cancer tissues (Fig. 6b), expression levels of LINC01010 were downregulated in M1 (metastatic samples) (Fig. 6c) and in stage IV tumors (Fig. 6d) compared to M0 (nonmetastatic samples) and stage I, respectively.
To further investigate the function of LINC01010, we performed co-expression analysis of LINC01010 using the Multi-Experiment Matrix (MEM) resource and selected the top 200 mRNAs co-expressed with LINC01010 for GO and KEGG pathway enrichment analysis (Table S7). The top 15 biological processes identified in GO and KEGG pathways revealed that mRNAs positively co-expressed with LINC01010 participated primarily in “integrin-mediated signaling pathway”, “extracellular matrix organization”, “leukocyte migration”, “regulation of actin cytoskeleton” and “proteoglycans in cancer” (Fig. 6e, f), suggesting that LINC01010 may regulate tumorigenesis and metastasis in lung cancer [26, 27].
LINC01010 represses lung cancer cell migration
To further determine whether LINC01010 affects lung cancer cell migration ability, we knocked down its expression by transfecting siRNA for LINC01010 into A549 or SPC-A-1 cells. qRT-PCR results showed that expression levels of LINC01010 were significantly decreased in three siRNA-transfected cells compared with controls (Fig. 7a, b). SiLINC01010-3 was selected for subsequent functional experiments. Transwell assays showed that suppression of LINC01010 promoted both A549 (Fig. 7c, e) and SPC-A-1 cell (Fig. 7d, f) migration and invasion compared to control cells, indicating that LINC01010 might play an important role in lung cancer metastasis.
LINC01010 associated ceRNA network
The ceRNA networks illustrated that LINC01010 may target hsa-mir-372, hsa-mir-373, hsa-mir-488 and hsa-mir-541. By Pearson correlation analysis, we found that LINC01010 was negatively correlated with hsa-mir-372 (Fig. 8a). By GSEA analysis, the gene set "KRAS.DF.V1_UP” (P < 0.0001) was significantly enriched in high levels of hsa-mir-372 (Fig. 8b), which was also verified by the literature [28, 29]. Meanwhile, the gene set "KRAS.DF.V1_UP” (P = 0.011) was significantly enriched in low levels of LINC01010 (Fig. 8c) and "KRAS.LUNG_UP.V1_DN" (P < 0.0001) was significantly enriched in high levels of LINC01010 (Fig. 8d), which suggested that LINC01010 might inhibit the function of hsa-mir-372 in the MAPK pathway.
We further investigated whether expression of LINC01010 was positively correlated with hsa-mir-372 target mRNAs. Pearson correlation analysis showed that LINC01010 was positively correlated with GRIA2, HS3ST4, SLC5A7 and TRDN (Fig. 8e). Therefore, LINC01010 may play an important role in tumorigenesis in lung cancer by competing with miR-372.