SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1, and FEZF1-AS1 are crucial long non-coding RNAs associated with the prognosis of lung squamous cell carcinoma.

Lung cancer has high incidence and mortality rates, and lung squamous cell carcinoma (LUSC) is a common form of non-small-cell lung carcinoma (NSCLC). The aim of our study was to discover long non-coding RNAs (lncRNAs) associated with LUSC prognosis. RNA-sequencing data obtained from LUSC samples were extracted from The Cancer Genome Atlas database. Using the limma package, differentially expressed genes (DEGs; including differentially expressed lncRNA genes (DELs), coding genes (DECs), and other genes (DEOs)) between LUSC and control samples were analyzed. Using Kaplan-Meier survival analysis, prognosis-associated lncRNAs were further selected. Following the calculation of Pearson's correlation coefficients between DELs and other DEGs, the DEL-DEG co-expression network was visualized using Cytoscape software. Using the clusterProfiler package, potential functions for DECs co-expressed with DELs were predicted. There were 1,305 DEGs in LUSC samples, including 153 DELs, 1,109 DECs, and 43 DEOs. Based on survival analysis, 22 prognosis-associated lncRNAs (including surfactant associated 1, pseudogene (SFTA1P), long intergenic non-protein coding RNA 968 (LINC00968), GATA6 antisense RNA 1, (GATA6-AS1) TBX5 antisense RNA 1 (TBX5-AS1) and FEZF1 antisense RNA 1 (FEZF1-AS1)) in LUSC were selected from these DELs, and the associated abnormal expression levels were also verified in LUSC clinical samples. A DEL-DEG co-expression network was constructed, which involved 93 DELs. Co-expressed DECs were enriched for only 8 prognosis-associated DELs, including LINC00968, SFTA1P, and TBX5-AS1. Specifically, mitogen-activated protein kinase (MAPK) signaling pathway-associated genes were enriched in DECs co-expressed with LINC00968, SFTA1P, GATA6-AS1, TBX5-AS1 and FEZF1-AS1, which may be prognosis-associated lncRNAs in LUSC. In addition, LINC00968 may affect the outcome of patients with LUSC via the MAPK signaling pathway. Copyright: © Xiong et al.

Introduction

Lung cancer is classified as small-cell lung carcinoma (SCLC) or non-small-cell lung carcinoma (NSCLC) (1). Patients with lung cancer are frequently effected by weight loss, shortness of breath, coughing and chest pains (2). Lung cancer has a high rate of incidence and mortality; in 2012, 1.8 million people were effected, resulting in 1.6 million deaths globally (3). NSCLC accounts for ~85 % of all lung cancer cases, and is less sensitive to chemotherapy than SCLC (4,5). As a type of NSCLC, lung squamous cell carcinoma (LUSC) has a higher incidence in males compared with females, and is closely associated with a history of smoking (6). It was therefore necessary to investigate the pathogenesis of LUSC in order to improve the clinical outcomes of patients, and approaches to treatment.

Long non-coding RNAs (lncRNAs) are non-coding RNAs longer than 200 nucleotides, which serve roles in various biological phenomena, including tumorigenesis and cancer progression (7,8). A number of lncRNAs associated with the outcome of patients with NSCLC have been reported; downregulation of HMlincRNA717 (also known as gastric cancer associated transcript 2) is associated with NSCLC progression, and HMlincRNA717 expression serves as an independent marker for predicting the survival of patients with NSCLC (9). Expression of TatD DNase domain containing 1 (TATDN1) is associated with the increased invasion and metastatic potential of NSCLC cells, indicating that TATDN1 may function as a promising prognostic biomarker and therapeutic target for the disease (10). Overexpression of SBF2 antisense RNA 1 (SBF2-AS1) is associated with lymph node metastasis, histological grade and poor overall survival, indicating that SBF2-AS1 is an independent prognostic marker and promising therapeutic target for patients with NSCLC (11). Decreased expression of tubulin alpha 4b has been reported in NSCLC cell lines and tissues, where it mediates cell proliferation and is an adverse prognostic predictor for NSCLC (12). The Pvt1 oncogene may be a candidate marker for the diagnosis and prognosis of patients with NSCLC, and may accelerate NSCLC cell proliferation by decreasing the expression of p15 and p21 (13,14). B-cell chronic lymphocytic leukemia (CLL)/lymphoma 2 is implicated in the favorable prognosis of patients with localized NSCLC via its interaction with the lncRNA metastasis-associated lung adenocarcinoma transcript 1 (15). However, elucidation of the roles of lncRNAs associated with the prognosis of patients with LUSC requires further investigation.

In order to determine the mechanisms of multiple underlying diseases, gene expression profiles of human tissues are increasingly being analyzed using bioinformatics tools (16). In the present study, RNA-sequencing data obtained from LUSC samples was downloaded from public databases. This was followed by differential expression analysis, survival analysis, co-expression network analysis and enrichment analysis to determine key lncRNAs effecting the prognosis of patients with LUSC. This may provide a basis for prognostic prediction and targeted therapy in LUSC.

Materials and methods

Data extraction and pre-processing

Using the R package of The Cancer Genome Atlas (TCGA) Biolinks package (http://www.bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html) (17), RNA-sequencing and clinical data for LUSC were extracted from TCGA http://cancergenome.nih.gov/) database. The raw data was pre-processed, and genes with no expression in >10% samples were defined as those with low expression, and subsequently filtered out. Based on the human reference genome GENCODE (hg 19) (18), ensemble gene IDs were mapped to gene symbols and the corresponding gene types were characterized. The GENCODE database classified the following types of genes as lncRNA genes: 3prime_overlapping_ncRNA, antisense_RNA, bidirectional_promoter_lncRNA, lincRNA, macro_lncRNA, non_coding, proessed_transcript, sense_intronic, sense_overlaping and TEC (18). According to this typing standard, genes from the RNA-sequencing data were grouped as lncRNA genes, coding genes and other genes.

Differential expression analysis

Data were normalized using the R package limma (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (19), and differential expression analysis for LUSC and the control samples was performed. From the genes obtained upon analysis, those with thresholds of false discovery rate (FDR) <0.01, and |log fold change (FC)| ≥3, differentially expressed genes (DEGs; including differentially expressed lncRNA genes (DELs), coding genes (DECs), and other genes (DEOs)) were selected.

Survival analysis for DELs

To analyze the prognostic features of the DELs, samples were classified into high and low expression groups with the median expression of DELs as the threshold. Kaplan-Meier (KM) survival analysis (20), was performed for clinical data to compare the survival differences between high and low expression groups. P<0.05 was considered to indicate a statistically significant difference.

Co-expression network and enrichment analysis

Co-expression network analysis is a valuable method of studying the biological functions of lncRNAs, and its theoretical basis is that the expression profiles of genes involved in the same biological function may be correlated (21). Pearson's correlation coefficients were determined between the DELs and DEOs, with the correlation coefficient |r|≥0.6 as the threshold for screening co-expression pairs. The DEL-DEG co-expression network was constructed using Cytoscape software v3.6.0 (http://www.cytoscape.org/) (22). Based on the results of the co-expression network, potential functions of prognosis-associated DELs were analyzed. Using the R package cluster Profiler (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (23), Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were predicted for those DECs co-expressed with DELs. P<0.05 was considered to indicate a statistically significant difference.

Patient sample collection

LUSC tissues and adjacent non-tumorous tissues were collected from 37 patients (age range, 40–80 years; mean age, 56 years; with 23 males and 14 females) with NSCLC who underwent surgical treatment at The Second Affiliated Hospital of Nanchang University (Nanchang, China) between May 2013 and September 2016. Following surgical resection, tissue samples were immediately frozen in liquid nitrogen and stored at −80°C. The study protocols conformed to the ethical guidelines of the 1975 Declaration of Helsinki, and were approved by the Institutional Ethical Review Committee of The Second Affiliated Hospital of Nanchang University. All patients enrolled in the study gave written informed consent.

RNA extraction and reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

Total RNA was extracted from patient tissue samples using TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and treated with DNase I (Roche Diagnostics, Indianapolis, IN, USA) according to the manufacturer's protocol. RT-qPCR was performed using the TaqMan High-Capacity cDNA Reverse Transcription Kit and the TaqMan Fast PCR Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.) with the corresponding primers, according to the manufacturer's protocol. GAPDH was used as an internal control. The RT-qPCR protocol included an initial denaturation step (95°C for 5 min) and 40 cycles of denaturation (95°C for 10 sec), annealing (60°C for 20 sec) and extension (72°C for 10 sec). The relative expression levels were calculated using the 2−ΔΔCq method as previously described. The primer sequences for surfactant associated 1, pseudogene (SFTA1P), long intergenic non-protein coding RNA 968 (LINC00968), GATA6 antisense RNA 1, (GATA6-AS1) TBX5 antisense RNA 1 (TBX5-AS1) and FEZF1 antisense RNA 1 (FEZF1-AS1) were as follows: SFTA1P forward, 5′-CAGGAGGTCACAGGGAAG-3′, and reverse, 5′-AAGCCGAATACAGTTGCC-3′; LINC00968 forward, 5′-GGCAGTTTTATTGTGGTGATT-3′, and reverse, 5′-AATGGAAGTTGACGGGATAG-3′; GATA6-AS1 forward, 5′-ACGGTTTCTGGACCTTCTC-3′, and reverse, 5′-TTGTGAACTTGTGGCTCCT-3′; TBX5-AS1 forward, 5′-GCAAAAGAAAGGGTTGGTC-3′, and reverse, 5′-AAAAGTGGGGAAACCAGAA-3′; FEZF1-AS1 forward, 5′-CTGGAGCCTTACCTGCCTT-3′, and reverse, 5′-CTGGAGGGACACACCTCAC-3′; MKP-1 forward, 5′-GGACATTTGGGCTGTGTG-3′, and reverse, 5′-CCGCTTTTGGACTGAGAGA-3′; TGF-β receptor type-2 (TGFBR2) forward, 5′-CCCCAGGTAAGGATAGCAG-3′, and reverse, 5′-CCAGGTAGGCAGTGGAAA-3′; and GAPDH forward, 5′-TCCTCTGACTTCAACAGCGACAC-3′, and reverse, 5′-CACCCTGTTGCTGTAGCCAAATTC-3′.

Statistical analysis

Statistical analysis was performed using GraphPad 6.0 statistical software (GraphPad Software, Inc., La Jolla, CA, USA). Data is shown as mean ± standard deviation. Data from two groups were analyzed using unpaired t-test and >2 groups were analyzed using one-way analysis of variance with Tukey's post hoc test. Spearman's correlation analysis was used to determine the correlations between the expression levels of LINC00968 and TGFBR2/MKP-1 in LUSC tissues. P<0.05 was considered to indicate a statistically significant difference.

Results

A total of 22,271 genes (including 3,685 lncRNAs, 16,045 coding genes and 2,541 other genes) from 551 samples (including 50 LUSC and 49 control samples) were obtained. Under FDR <0.01 and |logFC| ≥3, a total of 1,305 DEGs were identified (including 631 upregulated, and 674 downregulated genes). Among these DEGs, there were 153 DELs (99 upregulated and 54 downregulated), 1,109 DECs (504 upregulated and 605 downregulated) and 43 DEOs (28 upregulated and 15 downregulated). The volcano plot for the DEGs is displayed in Fig. 1. The clustering heat map suggested that the identified DELs may be used to distinguish between LUSC samples and controls (Fig. 2).

Figure 1.

Volcano plot of DEGs associated with lung squamous cell carcinoma. Red, blue and grey represent upregulated, downregulated and non-DEGs, respectively. DEG, differentially expressed gene.

Figure 2.

Heat map of differentially expressed long non-coding RNAs associated with lung squamous cell carcinoma. In the sample strip, red and blue represent tumor samples and non-tumor samples, respectively.

Survival analysis for identified DELs

According to the median expression level of the identified DELs, samples were divided into high and low expression groups. KM survival analysis for the two groups indicated that the expression of 22 DELs (including SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 and FEZF1-AS1) correlated with the prognosis of patients with LUSC (Table I). In addition, the expression of SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 and FEZF1-AS1 were also identified in LUSC tissues and normal tissues using RT-qPCR. As illustrated in Fig. 3, the expression levels of SFTA1P, LINC00968, GATA6-AS1 and TBX5-AS1 were significantly increased in LUSC tissues, and FEZF1-AS1 expression was significantly decreased in LUSC tissues compared with that in normal tissues (P<0.01). The KM survival curves for SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 and FEZF1-AS1 are displayed in Fig. 4. High expressions of SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 predicted poor prognosis of LUSC, while low expression of FEZF1-AS1 predicted a poor prognosis of LUSC.

Table I.

Differentially expressed long non-coding RNAs correlating with the prognosis of lung squamous cell carcinoma.

Gene	P-value
SFTA1P	0.034
LINC00968	0.047
AC0796304	0.024
TBX5-AS1	0.023
RP11-434D91	0.036
AC1337851	0.047
GATA6-AS1	0.002
RP11-613D138	0.013
RP11-385J12	0.032
CTD-2626G112	0.013
RP11-174G61	0.037
RP5-823G155	0.037
RP11-108M123	0.011
LINC00519	0.017
CTB-113P194	0.024
RP11-573D152	0.045
CDIPT-AS1	0.016
RP11-356K231	0.041
FEZF1-AS1	0.014
RP11-255A1121	0.015
RP11-357H1417	0.017
RP11-48O204	0.014

Figure 3.

Relative expression levels of long non-coding RNAs in tumor and adjacent, normal tissues. (A) SFTA1P, (B) LINC00968, (C) GATA6-AS1, (D) TBX5-AS1 and (E) FEZF1-AS1 in 37 paired LUSC tumor tissues and adjacent normal tissues were identified using reverse transcription-quantitative polymerase chain reaction. **P<0.01, ***P<0.001. SFTA1P, surfactant associated 1, pseudogene; LINC00968, long intergenic non-protein coding RNA 968; GATA6-AS1, GATA6 antisense RNA 1; TBX5-AS1, TBX5 antisense RNA 1; and FEZF1-AS1, FEZF1 antisense RNA 1.

Figure 4.

Kaplan-Meier survival curves for long non-coding RNA associations with patient survival. (A) SFTA1P, (B) LINC00968, (C) GATA6-AS1, (D), TBX5-AS1 and (E) FEZF1-AS1. Blue and red lines indicate high and low expression groups, respectively. SFTA1P, surfactant associated 1, pseudogene; LINC00968, long intergenic non-protein coding RNA 968; GATA6-AS1, GATA6 antisense RNA 1; TBX5-AS1, TBX5 antisense RNA 1; and FEZF1-AS1, FEZF1 antisense RNA 1.

The construction of a DEL-DEG co-expression network revealed 1,096 co-expressed pairs, including 93 DELs (Fig. 5). Notably, lncRNA-SFTA1P was co-expressed with the majority of the other genes identified. A total of 8 prognosis-associated DELs (including LINC00968, SFTA1P and TBX5-AS1) were enriched. The GO terms and KEGG pathways outlining the enrichment of LINC00968 are presented in Fig. 6, and details of the MAPK signaling pathway are highlighted in Fig. 7. The expression levels of TGFBR2 and MKP-1, which are involved in MAPK signaling in 37 LUSC tumor samples, were identified using RT-qPCR. Fig. 8 demonstrates that the expression of TGFBR2 negatively correlated with that of LINC00968 (P<0.001), whilst the expression of MKP-1 positively correlated with LINC0096 expression (P<0.001).

Figure 5.

Co-expression network of genes implicated in lung squamous cell carcinoma. Red and blue represent high and low gene expression, respectively. Squares indicate long non-coding RNA genes, and circles represent coding genes or other genes.

Figure 6.

(A) Functional terms and (B) pathway enrichment of long non-coding RNA long intergenic non-protein coding RNA 968. Horizontal and vertical coordinates represent the count of implicated genes and term/pathway names, respectively. Red and blue indicate low and high P-values, respectively.

Figure 7.

The mitogen-activated protein kinase signaling pathway in lung squamous cell carcinoma. Red and green represent upregulated and downregulated genes, respectively.

Figure 8.

Correlation between LINC00968 and TGFBR2/MKP-1 expression in LUSC. (A) Expression of TGFBR2 in LUSC tissues was identified by RT-qPCR analysis, and revealed a negative correlation with LINC00968 expression. (B) Expression of MPK-1 in LUSC tissues revealed a positive correlation with LINC00968 expression. LINC00968, long intergenic non-protein coding RNA 968; TGFBR2, TGF-β receptor type-2; MPK-1, mitogen-activated protein kinase 1; LUSC, lung squamous cell carcinoma; RT-qPCR, reverse transcription-quantitative polymerase chain reaction.

Discussion

In the present study, a total of 1,305 DEGs between LUSC samples and control samples were selected (including 153 DELs, 1,109 DECs and 43 DEOs). Survival analysis revealed that 22 DELs (including SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 and FEZF1-AS1) were prognosis-associated lncRNAs for LUSC. Additionally, a DEL-DEG co-expression network involving 93 DELs was constructed. Enrichment analysis for the DECs co-expressed with the DELs suggested that only 8 prognosis-associated DELs (including LINC00968, SFTA1P and TBX5-AS1) were enriched.

SFTA1P and cancer susceptibility candidate 2 mediate the functions of tumor suppressor genes and oncogenes during the development of LUSC, which may be used in diagnosis, prognostic prediction and targeted therapy (24). Increased expression of SFTA1P results in apoptosis and promotes cisplatin sensitivity in LUSC cells, indicating that SFTA1P may be considered a promising diagnostic marker, and an indicator for cisplatin sensitivity in patients with LUSC (25). Downregulation of GATA6-AS1 results in adverse outcomes for patients with gastric cancer (26), and also serves a tumor suppressing role in CLL (27). GATA6 has a regulatory effect in alveolar epithelial cells (28,29), and during lung development, GATA6-AS1 expression has been associated with that of GATA6 (30). These reports suggest that SFTA1P and GATA6-AS1 may be implicated in the prognosis of LUSC.

TBX5 expression is significantly associated with the disease stage, lymph node status and histopathological type of NSCLC, and serves a tumor-suppressing role in the disease (31,32). FEZF1-AS1 overexpression is associated with disease stage, differentiation degree and lymph node metastasis of NSCLC, and thus FEZF1-AS1 may promote tumor formation and be considered as a therapeutic target for the disease (33). Furthermore, FEZF1-AS1 expression levels were significantly upregulated in lung adenocarcinoma (LAD) samples, which resulted in unfavorable prognosis. This suggests that FEZF1-AS1 dysregulation promotes the development and progression of LAD, and may therefore be used in LAD therapy (34). The overexpression of FEZF1-AS1 is apparent in LAD tissues and cell lines, which may affect the cell cycle and proliferation of LAD cells by inhibiting the expression of p57 (35). Therefore, TBX5-AS1 and FEZF1-AS1 may also be associated with the clinical outcome of patients with LUSC.

LINC00968 inhibition suppresses the growth, invasion and migration of NSCLC cells, and is therefore an notable molecular marker for the diagnosis, prognosis and treatment success of patients with NSCLC (36). By activating the p38 MAPK and c-Jun N-terminal kinase signaling pathways, and suppressing the phosphatidylinositol 3-kinase/Akt/mammalian target of rapamycin signaling pathway, treatment with platycodin-D results in autophagy in A549 and NCI-H460 cells, and serves as an alternative compound for the treatment of NSCLC (37). Through the extracellular signal-regulated protein kinase/MAPK signaling pathway, claudin-7 suppresses the migration and invasion of growth factor-stimulated lung cancer cells (38). Pathway enrichment analysis revealed that LINC00968 was enriched for the MAPK pathway, thus, it was speculated that LINC00968 may affect the survival of patients with LUSC by regulating MAPK signaling.

In conclusion, 1,305 DEGs were identified from LUSC samples relative to corresponding control samples. The lncRNAs-SFTA1P, LINC00968, GATA6-AS1, TBX5-AS1 and FEZF1-AS1 may be associated with the prognosis of patients with LUSC. However, the effects of these key lncRNAs on patient survival rate requires further experimental verification.