Metastasis-related long non-coding RNAs AL359220.1, SH3BP5-AS1 and ZF-AS1 are significant for prognostic assessment of lung adenocarcinoma

Acquisition of MRLNRs

Transcriptome RNA-sequencing data and clinical data were downloaded from TCGA database. As a standard of P< 0.05 and |log2 fold change| > 1, we screened 3278 DElncRs, of which 2426 were upregulated and 852 were downregulated, and the top 50 upregulated and downregulated DElncRs were illustrated in Figure 1A. Following that, lncRNAs and mRNAs data were extracted from transcriptome data. We screened 44 MRGs in TOMIDA_METASTASIS_UP M17830 and TOMIDA_METASTASIS_DN M2583 of Molecular Signatures Database, of which 339 lncRNAs were identified to be the MRLNRs through correlation analysis.

The correlation between prognosis of LUAD patients and MRLNRs

Based on COX regression analysis, we then verified 6 MRLNRs which were associated with prognosis of LUAD patients, including AL359220.1, AL021368.2, SH3BP5-AS1, AL590666.2, AL109811.2 and ZF-AS1 (p<0.01). The relationships between these sMRLNRs and LUAD prognoses were clearly illustrated in the forest plot (Figure 2).

Prognostic features of the high-risk group and the low-risk group

The three sMRLNRs (AL359220.1, SH3BP5-AS1 and ZF-AS1) among the 6 sMRLNRs were selected by multivariate COX analysis (P<0.05)., and were used to establish the MRRS, by which the LUAD patients were divided into the high-risk group and the low-risk group (Figure 3A). As illustrated in the Figure 3B, the mortality rate constantly increased with the higher risk score. With the increase of the risk score, the expression levels of ZF-AS1 were increased, while SH3BP5-AS1 and AL359220.1 were decreased (Figure 3C). Furthermore, the survival curve of the high-risk group was significantly lower than that of the low-risk group (Figure 4).

The relationships between MRRS and clinical features and the regulation network of sMRLNRs

To further explore the relevance of the sMRLNRs and clinical features of LUAD, we analyzed the correlation of MRRS and the clinical and demographic characteristics including stage, T-stage and N-stage and M-stage (Figure 5A–5H). We found that the MRRS was significantly correlated with N-stage and M-stage (Figure 5C, 5D). We found the expression levels of SH3BP5-AS1 and AL359220.1 were lower in LUAD patients with the more advanced stage and T-stage (Figure 5E, 5F). Furthermore, the expression of SH3BP5-AS1 was correlated with early N-stage. Besides, the expression of ZF-AS1 was higher with the more advanced stage, T-stage, N-stage and M-stage (Figure 5E–5H). We further displayed the independent prognostic factor analysis, the results showed stage, T-stage, N-stage, M-stage and risk score were prominently correlated with OS in univariate analysis (P<0.05). However, only risk score illustrated as an independent risk factor in the multivariate analysis (Table 1). The ROC curves could represent the accuracy of the MRRS.

The AUCs of risk score was 0.661, sensitivity:0.215 and specificity:0.912 (Figure 6). These results suggested that the MRRS can be regarded as an accuracy predict model for LUAD. Next, the MRRS was integrated into a nomogram for predicting 1-, 3-, and 5-year survival rates of LUAD patients (Figure 7). Based on the nomogram scores, the 1-, 3-, and 5-year survival rates of LUAD patients could be well predicted by their nomogram scores. To further find the regulatory relationships of sMRLNRs and their target microRNAs, we selected the microRNAs of Starbase database and the results were demonstrated in regulatory network (Figure 8).

LncR-ZFAS1 showed significant relevance to metastasis

Then, we examined the levels of these sMRLNRs which were used to establish the MRRS in various LUAD tissues and cell lines, and found that ZFAS1 expressed highly in tumor tissues and showed the inverse results compared to SH3BP5-AS1 and AL359220.1 (Figure 9A–9C). Consistently, the expression of ZFAS1 in LUAD cell lines (A549, NCI-H1299 and HCC827) was dramatically higher than that in normal bronchial epithelial cell line, while SH3BP5-AS1 and AL359220.1 expressed increasingly in HBE cell line (Figure 9D–9F).

In addition, we further investigated the clinical relevance of ZFAS1. All the LUAD patients were divided into the ZFAS1-high-expression group and ZFAS1-low-expression group, followed by the analyses of their associations with various clinicopathologic features. As illustrated in Table 2, ZFAS1 showed prominent correlation to T-stage and metastatic status, which was consistent with the finding in Figure 5C, 5D.

LncR-ZFAS1 strikingly promoted migration and invasion abilities of LUAD cell lines

To further validate the potency of ZFAS1 in regulating the metastatic process of LUAD, ZFAS1 was silenced using siRNA in A549 cell line with the highest expression (Figure 10A). Interestingly, we found that siRNA-ZFAS1 considerably decreased the invasion and migration abilities of A549 cell line, suggesting the pro-metastatic effect of ZFAS1 (Figure 10B). Consistently, siRNA-ZFAS1 provoked the remarkable reductions of VIMENTIN and NCADHERIN in the epithelial-mesenchymal transition pathway, but showing the opposite role to the E-CADHERIN (Figure 10C). These findings indicate that lncR-ZFAS1 plays a pivotal role in promoting LUAD metastasis.

Sample collection and cell culture

One hundred and fifty-six LUAD tissues and adjacent normal tissues were collected from patients admitted to the First Affiliated Hospital of Chongqing Medical University between May 2018 and September 2022. All specimens were frozen in liquid nitrogen immediately and transferred to –80° C refrigerator for further analyses.

LUAC cell lines (A549, NCI-H1299 and HCC827) and human bronchial epithelial cell line (HBE) were purchased from the American Type Culture Collection (Manassas, VA, USA). DMEM and 1640 basic medium, supplementing with 10% fetal bovine serum, 100 μ/mL penicillin and streptomycin (Gibco, Gaithersburg, MD, USA) was utilized for cell culture. Cells were incubated at 37° C with an atmosphere of 5% CO₂.

Real-time quantitative PCR

According to the manufacturer's instructions, total RNA from tissues and cell lines under various experimental conditions was extracted by Triazole (Invitrogen, Waltham, MA, USA). cDNA Synthesis Kit (TaKaRa, Osaka, Japan) combining with RNA (1μg) was utilized to reverse transcribed cDNA. The quantitative polymerase chain reaction (qPCR) was performed on an ABI 7500 real-time PCR system (Applied Biosystems, Waltham, MA, USA) using SYBR-Green method (TaKaRa). Relative expression levels of lncRNAs normalized to GAPDH was calculated by the 2−ΔCt method. The primer sequences are shown in Table 3. Three assays were performed per cDNA sample.

Transcriptome data download and preprocessing

RNA-sequencing data and clinical data of LUAD patients were downloaded from the TCGA data portal (https://portal.gdc.cancer.gov/), containing 402 LUAD and 54 non-tumor samples (Supplementary Table 1). Our data were obtained from the TCGA database. We conducted preprocessing of the clinical data, excluding patients with a survival period of less than 30 days from the analysis. This exclusion was done because this subgroup of patients may have experienced early mortality due to factors such as bleeding and infection. These data were currently updated in Dec.12, 2022. Raw data were collected to do further analyses. RNA-seq results and clinical results were combined into a matrix file by a merge script of the Perl (http://www.perl.org/).

MRLNRs extraction

The limma package of R software was used to screen differentially expressed lncRNAs (DElncRs) from TCGA dataset by comparing the LUAD and non-tumor samples. The screen criterion is as follows: P< 0.05 and |log2 fold change| > 1. MRGs were screened by The Molecular Signatures Database v4.0 (TOMIDA_METASTASIS_UP M17830, TOMIDA_METASTASIS_DN M2583, http://www.broadinstitute.org/gsea/msigdb/index.jsp). MRGs were used to establish the immune score of LUAD gene by GSEA. The correlations between immune score and the expression levels of lncRNAs in LUAD patients were calculated by Pearson correlation analysis. A standard of |r|>0.6 and P<0.01 was used to verify the MRLNRs.

Acquire the survival-related MRLNRs (sMRLNRs)

MRLNRs with correlation of OS were regarded as sMRLNRs in LUAD patients. The sMRLNRs were screened by univariate COX analysis using R software survival packages (P<0.01). Besides, the sMRLNRs were divided into protective and deleterious portion by Hazard ratio. These sMRLNRs were used for subsequent research.

Metastasis-related risk score model (MRRS) creation

The sMRLNRs have been analyzed by multivariate COX analysis, and the integrated sMRLNRs were utilized as an independent prognostic factor to develop the MRRS (P<0.05). We performed the MRRS to classify LUAD patients into the high-risk group and the low-risk group.

The creation of MRRS was based on the expression data multiplied by Cox regression coefficients. The formula was as follows, [Expression levels of AL359220.1 * (-0.27974)] + [Expression levels of SH3BP5-AS1 * (-0.16829)] + [Expression levels of ZF-AS1 * (0.142708)]. Patients were divided into high-risk group and low-risk group according to the median score of the MRRS.

Bioinformatics analysis

Kaplan-Meier curve was used to evaluate the OS between high-risk group and low-risk group. Univariate Cox regression analysis and Pearson correlation analysis were used to identify the interest MRLNRs. Univariate and multivariate Cox regression analysis were used for identify the independent prognostic factors of LUAD patients. Based on the MRRS, ROC curve was employed to assess the sensitivity and specificity of the prognosis. Gene set enrichment analysis (GSEA) was displayed to explore the different functional phenotypes between the high-risk group and the low-risk group. The nomogram was constructed by rms package of R software to provide a reference for clinical evaluation of LUAD patients’ prognosis.

Statistical analysis

SPSS21.0 software (SPSS Inc, Chicago, IL, USA) and GraphPad Prism5 (GraphPad Software Inc, La Jolla, CA, USA) were employed to analyze data. Data were expressed as means ± SD. The relevance between lnR-ZFAS1 and clinicopathological features of LUAD patients was detected by using Fisher’s exact probability method. The differential comparison between two or more groups were analysed through Student T-test or ANOVA and post-hoc test respectively. P<0.05 was considered a significantly statistical difference.

Availability of data and materials

Authors can provide all of the datasets on reasonable request.