This article is mentioned in:

Introduction

Lung cancer is responsible for ~26% of all cancer mortalities worldwide (1). Despite the continuous increase in survival rates for most types of cancer, the 5-year survival rate of lung cancer has remained low, at only 18% (1). Lung cancer can be divided pathologically into two types. Most lung cancer cases (>80%) correspond to non-small cell lung cancer (NSCLC), with the remainder being small-cell lung cancer (2). The pathogenesis of NSCLC is not fully understood, although studies about the molecular mechanisms of lung cancer continue to provide new insights into the disease.

Long noncoding RNAs (lncRNAs) are attracting increasing attention for their potential roles in numerous cancer types (3,4). For example, in NSCLC, two lncRNAs-plasmacytoma variant translocation 1 (PVT1) and urothelial carcinoma-associated 1 (UCA1)-were found to promote oncogenesis and tumor growth, respectively (5,6). LncRNAs may participate in tumorigenesis through the competitive endogenous RNA (ceRNA) network (7,8) by interacting with microRNAs (miRNAs) (9). Differential miRNA expression has been detected in precancerous and cancerous tissues (10,11). Therefore, it is possible that the ceRNA network might be important in carcinogenesis and cancer development. Furthermore, knowledge about ceRNA interactions might provide insights into the potential mechanisms and biological functions of lncRNAs in cancer. One study demonstrated that the carcinogen H19 acted as a ceRNA to alter expression levels of certain miRNAs (miR-138 and miR-200a) and genes [zinc finger E-box binding homeobox 1 (ZEB)1/ZEB2] in colorectal cancer (8). In addition, ceRNA networks have been identified in other cancer types, including hepatocellular, breast, pancreatic and gastric cancer (12–15).

Using bioinformatic analysis, it was previously demonstrated by the authors that the lncRNA LINC00968 was differentially expressed between normal lung and tumor tissues (16). In the present study, microarray analysis was employed to identify which miRNAs were differentially expressed with LINC00968 overexpression. Among the five identified differentially expressed miRNAs, miR-9-3p was selected for further analysis. In total, nine prediction algorithms were applied to identify target genes of miR-9-3p, which were then subjected to functional enrichment analysis (17–20). Finally, the lncRNA-miRNA-mRNA regulatory axis of LINC00968 was validated by bioinformatic and correlation analyses.

Materials and methods

miRNA microarray profiling of NSCLC cells

The NSCLC A549 cell line was obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). A549 cells are adenocarcinoma human alveolar basal epithelial cells, maintained in F-12 culture medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA), trypsin-EDTA solution (Gibco; Thermo Fisher Scientific, Inc.) and D-Hank's solution (Shanghai GenePharma Co., Ltd., Shanghai, China). A549 cells were transduced with lentiviruses encoding LINC00968 for 24–72 h, then transduction efficiency was examined by an inverted fluorescent microscope and subsequent experimentation was performed. Considering the effects of green fluorescent protein (GFP)-containing lentiviral transfection on subsequent experiments such as apoptosis detection, a GFP-free lentivirus was selected. Optimal transfection conditions (achieving 95% efficiency) were when the lentivirus titer reached 1×108 TU/ml, diluted to 1:5. miRNA profiles were determined by using miRCURY LNA Array (version 18.0; Exiqon, Inc., Woburn, MA, USA) on an Axon GenePix 4000B microarray scanner (Axon Instruments; Molecular Devices, LLC, Sunnyvale, CA, USA). Results of the microarray analysis indicated significant downregulation of miR-9-3p (log2|fold change (FC)| >0.585, P<0.05).

Comparison with Gene Expression Omnibus (GEO) dataset

Data for miR-9-3p were searched and assembled in the GEO database (https://www.ncbi.nlm.nih.gov/gds/; date accessed: 4th June 2017), using the following search terms: (Lung OR pulmonary OR respiratory OR bronchi OR bronchioles OR alveoli OR pneumocytes OR ‘air way’) AND (cancer OR carcinoma OR tumor OR neoplas* OR malignan* OR adenocarcinoma) AND (microRNA OR miRNA OR ‘micro RNA’ OR ‘small temporal RNA’ OR ‘noncoding RNA’ OR ncRNA OR ‘small RNA’). For microarrays, the following inclusion criteria were used: i) Lung adenocarcinoma and adjacent noncancerous tissues (or normal lung tissues) were included in each dataset; ii) sample organism was Homo sapiens; and iii) expression levels of miR-9-3p (hsa-miR-9* or hsa-miR-9-3p) in the experimental and control groups were provided or could be calculated.

The mean and standard deviations of expression values were extracted to estimate miR-9-3p levels in the test and control groups using STATA (version 12.0; StataCorp LP, College Station, TX, USA). P<0.05 or I2 value >50% indicated significant heterogeneity. Standard mean differences (SMDs) with 95% confidence intervals (CIs) were used to evaluate continuous outcomes. Funnel plots were generated to evaluate publication bias. P<0.05 was considered to indicate statistical significance.

Bioinformatic prediction of miRNA targets and functional enrichment analysis

The putative target genes of miR-9-3p were predicted using nine miRNA prediction algorithms: TargetScan (http://www.targetscan.org/), miRDB (http://www.mirdb.org/), RNA22 (https://cm.Jefferson.edu/rna22/), miRMap, mirTarBase (http://mirtarbase.mbc.nctu.edu.tw/), DIANA-microT (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/iinde), GeneCards (http://www.genecards.org), TarBase (http://diana.imis.athenainnovation.gr/DianaTools/index.php?r=tarbase/index) and TargetMiner (https://www.isical.ac.in/~bioinfo_miu/final_html_targetminer/hsa-miR-9-3p.html). Only target genes found in at least three databases were considered for further analysis. To clarify potential roles of the target genes, the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) was used to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. The method used to identify significant regulatory functions and enrichment pathways was based on the P-value as previously described (21).

Bioinformatic analysis of the association between LINC00968 and miR-9-3p using The Cancer Genome Atlas (TCGA) database

Level 3 RNAseq and miRNAseq data of lung adenocarcinoma were downloaded from TCGA (http://cancergenome.nih.gov/). An R package was used to distinguish differentially expressed RNAs and miRNAs between lung cancer and adjacent non-tumor lung tissues. For all P-values, the false discovery rate (FDR) was applied to correct for multiple hypothesis testing. The threshold of log2|FC|≥1 and FDR <0.05 was considered significant. Genes with an absolute correlation coefficient >0.1 were regarded as correlated genes.

Construction of protein-protein interaction (PPI) network

Next, how closely connected the predicted target genes were with LINC00968 and miR-9-3p was determined. The predicted genes were integrated with the correlated genes to obtain prospective objective genes. To analyze interactions of the objective genes, the Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org) was used to construct the PPI network (22), using a medium CI of 0.4 as a threshold.

Statistical analysis

The SPSS (version 22.0; SPSS Inc., Chicago, IL, USA) and R (version 3.4.0) (23) software packages were used to analyze all experimental data. Begg's test and Egger's test were used to detect publication bias in meta-analysis. Receiver operating characteristic (ROC) curves were constructed to evaluate the sensitivity and specificity of the molecular biological indicator in diagnostic tests. Differences in overall survival between two groups were estimated and compared using the Kaplan-Meier method with the log-rank test. FDR was used to correct the P-value for multiple testing. Pearson's correlation analysis was used to identify the association between LINC00968, miR-9-3p and CCNA2 in TCGA database. P<0.05 was considered to indicate a statistically significant difference.

Results

Downregulation of miR-9-3p expression in LINC00968-overexpressing A549 cells

The miRNAs that were differentially expressed in LINC00968-overexpressing A549 cells were profiled using miRNA microarray analysis. A total of five differentially expressed miRNAs were identified: One miRNA (miR-3675-3p) was upregulated and four miRNAs (miR-9-3p, miR-22-5p, miR-668-3p and miR-4536-3p) were downregulated in LINC00968-overexpressing A549 cells. miR-9-3p (FC=0.569, P=0.044) was selected for further analysis.

Confirmation of miR-9-3p downregulation in the GEO database

A meta-analysis was performed to verify miR-9-3p expression in other lung adenocarcinoma tissue microarrays. A total of seven eligible microarrays were identified. As significant heterogeneity was detected among the microarrays (I2=59.3%, P<0.05), the random-effects model was used to assess the pooled SMD and 95% CIs. Pooled SMD was 0.86 (95% CI, 0.40–1.32, P=0.022; Fig. 1A), which suggested that the expression of miR-9-3p in lung adenocarcinoma tissues was higher compared with normal lung tissues. Funnel plots revealed no significant publication bias in the meta-analysis (Begg's test, P=0.681; Egger's test: P=0.681; Fig. 1B).

Figure 1. Forest plot and Begg's funnel plot of miR-9-3p expression in microarrays from GEO database. (A) Forest plot evaluating differences in miR-9-3p expression between lung adenocarcinoma tissues and normal lung tissues. Experimental group, high miR-9-3p expression; control group, low miR-9-3p expression. (B) Funnel plot for publication bias test of GEO microarrays. CI, confidence interval; GEO, Gene Expression Omnibus; miR, microRNA; s.e., standard error; SMD, standard mean difference.

Functional enrichment analysis of miR-9-3p target genes

A total of nine miRNA target prediction algorithms were used to identify target genes of miR-9-3p. Of the 7856 predicted target genes, 2047 genes overlapped in at least three datasets. To elucidate the biological function of miR-9-3p, functional enrichment analysis was implemented using these 2,047 target genes. Positive regulation of transcription from RNA polymerase II promoter (GO:0045944) was the most significant biological process term. Cytoplasm (GO:0005737) was the most significant cellular component term, and protein binding (GO:0005515) was the most significant molecular function term (Table I). A schematic of the top five most significant GO processes are presented in Fig. 2A.

Figure 2. Results of functional enrichment analysis. (A) Top five most enriched Gene Ontology terms. (B) Top five most enriched Kyoto Encyclopedia of Genes and Genomes pathways. Blue, high significance; red, low significance; AMPK, AMP-activated protein kinase; FoxO, forkhead box O1.

Table I. Significant GO terms most strongly enriched by target genes of miR-9-3p.

Table I.

Significant GO terms most strongly enriched by target genes of miR-9-3p.

GO term	Count	P-value	FDR
Biological processes
GO:0045944: Positive regulation of transcription from RNA polymerase II promoter	175	1.90×10−11	3.62×10−8
GO:0045893: Positive regulation of transcription, DNA-templated	97	6.97×10−8	1.33×10−4
GO:0000122: Negative regulation of transcription from RNA polymerase II promoter	123	3.09×10−7	5.87×10−4
GO:0006366: Transcription from RNA polymerase II promoter	93	7.81×10−7	1.48×10−3
GO:0006351: Transcription, DNA-templated	277	1.03×10−6	1.95×10−3
Cellular components
GO:0005737: Cytoplasm	690	3.07×10−13	4.72×10−10
GO:0005654: Nucleoplasm	400	6.87×10−12	1.05×10−8
GO:0005634: Nucleus	703	7.26×10−12	1.11×10−8
GO:0005923: Bicellular tight junction	29	1.26×10−5	1.93×10−2
GO:0016020: Membrane	291	1.78×10−5	2.73×10−2
Molecular function
GO:0005515: Protein binding	1133	1.65×10−21	2.72×10−18
GO:0003730: mRNA 3′-UTR binding	22	8.11×10−9	1.34×10−5
GO:0043565: Sequence-specific DNA binding	97	5.74×10−8	9.48×10−5
GO:0003700: Transcription factor activity, sequence-specific DNA binding	156	1.03×10−7	1.71×10−4
GO:0001077: Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	54	1.36×10−7	2.25×10−4

[i] FDR, false discovery rate; GO, Gene Ontology; miR, microRNA; UTR, untranslated region.

The biological relevance of KEGG pathways that were enriched in these target genes was analyzed. The targets were involved in many cancer-associated pathways (Table II). The top five KEGG pathways that were enriched in the target genes are listed in Fig. 2B. KEGG pathway analysis identified many important signaling pathways [e.g., AMP-activated protein kinase (AMPK) and forkhead box O1 signaling pathways] known to be involved in the development, invasion, and metastasis of cancer. These findings indicate the biological relevance of miR-9-3p.

Table II. Kyoto Encyclopedia of Genes and Genomes pathways most strongly enriched by target genes of miR-9-3p.

Table II.

Kyoto Encyclopedia of Genes and Genomes pathways most strongly enriched by target genes of miR-9-3p.

Terms	Count	P-value	FDR
hsa04152: AMPK signaling pathway	33	3.61×10−7	4.74×10−4
hsa05215: Prostate cancer	27	3.74×10−7	4.91×10−4
hsa04068: FoxO signaling pathway	34	1.14×10−6	1.49×10−3
hsa05200: Pathways in cancer	70	3.77×10−6	4.95×10−3
hsa05218: Melanoma	22	4.87×10−6	6.39×10−3
hsa04390: Hippo signaling pathway	35	6.78×10−6	8.91×10−3
hsa04550: Signaling pathways regulating pluripotency of stem cells	32	2.45×10−5	3.22×10−2
hsa05214: Glioma	18	2.10×10−4	2.75×10−1
hsa03015: mRNA surveillance pathway	22	2.70×10−4	3.54×10−1
hsa04151: PI3K-Akt signaling pathway	57	2.87×10−4	3.77×10−1

[i] AMPK, AMP-activated protein kinase, catalytic, α-1; FDR, false discovery rate; FoxO, forkhead box O1; miR, microRNA; PI3K, phosphoinositide-3-kinase.

PPI network

To understand the interactions between target genes and their associations with LINC00968 and miR-9-3p, the predicted genes were integrated with the correlated genes. A total of 6,515 and 2,677 genes were correlated with LINC00968 and miR-9-3p, respectively, from the TCGA database. Using a Venn diagram, the 2,047 predicted target genes were integrated with the 6515 LINC00968-associated genes and the 2677 miR-9-3p-associated genes, uncovering 120 objective genes. To depict the complex PPI relationships among these objective genes, the PPI network of the top 30 genes were visualized by STRING (Fig. 3A). The number of links between genes was calculated, and a list of the top 30 most connected genes was obtained. Among them, cyclin A2 (CCNA2) had the greatest connectivity (Fig. 3B) and was subjected to further validation.

Figure 3. Analysis of functional protein association. (A) PPI network of 120 objective genes. (B) Top 30 most connected genes of PPI that are arranged by connectivity. PPI, Protein-protein interaction.

Validation of regulatory axis of LINC00968 in TCGA database

Predictive performances of LINC00968, miR-9-3p and CCNA2 in lung adenocarcinoma were estimated by ROC curve analysis. Areas under the ROC curve (AUCs) of components in this regulatory axis among lung adenocarcinoma patients from TCGA database were as follows: LINC00968, AUC=0.988 (95% CI, 0.977–0.998, P=0.000), miR-9-1, AUC=0.983 (95% CI, 0.967–0.999, P=0.000), miR-9-2 AUC=0.982 (95% CI, 0.968–0.997, P=0.000), miR-9-3 AUC=0.981 (95% CI, 0.968–0.995, P=0.000), and CCNA2, AUC=0.960 (95% CI, 0.921–0.999, P=0.000). The data suggested that the three signatures had good predictive performance in lung adenocarcinoma patients.

A Kaplan-Meier survival analysis was performed to determine the effects of the three components on predicting clinical outcome (Fig. 4). The median gene expression (which, for LINC00968, CCNA2, miR-9-1, miR-9-2 and miR-9-3 were 20.70, 694.90, 1009.58, 986.66 and 1010.02, respectively) was used to separate the patients into high and low expression groups in TCGA datasets. The number of patients expressing LINC00968 and CCNA2 in the high-expression and low-expression groups in the present study was 267. The number of patients expressing miR-9-3p in the high-expression and low-expression groups was 260. Low expression of LINC00968 (Fig. 4A) or high expression of CCNA2 (Fig. 4E) was a significant predictor of poor prognosis for overall survival in lung adenocarcinoma patients. These findings suggested that LINC00968 and CCNA2 had good potential as diagnostic and prognostic indicators.

Figure 4. Overall survival analysis of LINC00968, miR-9-3p and CCNA2 expression in lung adenocarcinoma from The Cancer Genome Atlas database. Survival curves for (A) LINC00968, (B) miR-9-1, (C) miR-9-2, (D) miR-9-3 and (E) CCNA2. Red and green lines indicate high and low levels of expression, respectively. x-axis, overall survival time (year). y-axis, overall survival probability. The number of patients expressing LINC00968 and CCNA2 in the high-expression and low-expression groups is 267. The number of patients expressing miR-9-3p in the high-expression and low-expression groups is 260. CCNA2, cyclin A2; miR, microRNA.

Finally, the underlying mechanisms of the LINC00968 regulatory axis in lung adenocarcinoma were investigated. In our previous GEO results, miR-9-3p was detected to be significantly upregulated in tumor tissues and had oncogenic potential in lung adenocarcinoma. Therefore, the correlation between LINC00968 and miR-9-3p (including miR-9-1, miR-9-2, and miR-9-3 in TCGA database) was examined. Negative correlations were identified between LINC00968 and miR-9-3p (Fig. 5A-C) and between LINC00968 and CCNA2 (Fig. 5D), whereas a positive correlation between miR-9-3p and CCNA2 (Fig. 5E-G) was identified in patients with lung adenocarcinoma. Therefore, it is concluded that the regulatory role of LINC00968 in lung adenocarcinoma might involve the LINC00968/miR-9-3p/CCNA2 axis.

Figure 5. Correlation analysis among LINC00968, miR-9-3p, and CCNA2 expression in lung adenocarcinoma based on TCGA database. Correlations (A) between LINC00968 and miR-9-1, (B) between LINC00968 and miR-9-2, (C) between LINC00968 and miR-9-3, (D) between LINC00968 and CCNA2, (E) between miR-9-1 and CCNA2, (F) between miR-9-2 and CCNA2, (G) between miR-9-3 and CCNA2. CCNA2, cyclin A2; miR, microRNA; TCGA, The Cancer Genome Atlas.

Discussion

In recent years, microarray studies, RNA sequencing, and quantitative reverse-transcription PCR have been used to investigate the differential expression profiles of lncRNAs in various types of cancer. LncRNAs regulate diverse biological processes and act as oncogenes or tumor suppressors in many cancer types (24–28). As such, lncRNAs have great potential for use as molecular biomarkers in cancer diagnosis and prognosis. For example, the lncRNA, metastasis associated in lung adenocarcinoma transcript 1, was found to promote the proliferation and metastasis of lung cancer (29). The authors previously reported that the lncRNA LINC00968 was differentially expressed in normal vs. lung tumor tissues in the GEO database (16). However, the regulatory mechanisms of LINC00968 in lung cancer have been unclear.

LncRNAs and miRNAs may have important interactions that contribute to cancer development. Evidence for this ceRNA hypothesis has been provided by several reports. For instance, lncRNA cancer susceptibility candidate 2 (CASC2) was shown to inhibit cell proliferation and tumor growth by regulating miR-18a in colorectal cancer (30). C032469 regulated telomerase reverse transcriptase expression by binding to miR-1207-5p in gastric cancer (31). Growth Arrest Specific Transcript 5 suppressed carcinogenesis in NSCLC by repressing miR-135b expression (32). In light of this research, the present authors postulated that LIN00968 may act as a ceRNA in NSCLC. The expression profile of miR-9-3p in NSCLC A549 cells that overexpressed LINC00968 was detected. Finally, bioinformatic analysis was used to reveal the potential regulatory role of LINC00968 in lung adenocarcinoma.

miR-9-3p has been demonstrated to regulate biological functions in multiple types of cancer. For example, Yang et al (33) reported that H2O2-induced apoptosis of glioma cells was facilitated by miR-9-3p through downregulation of Herpud1 expression. In another study, miR-9-3p inhibited cell proliferation by targeting the PDZ-binding motif TAZ in hepatocellular carcinoma (34). miR-9-3p was detected to suppress tumor biological behaviors in gastric cancer (35). Other studies have suggested regulatory roles for miR-9-3p in tumor initiation, growth and progression in NSCLC (36–38).

Results of the present analysis of miR-9-3p target genes indicated that the most enriched GO term was protein binding, a key process in carcinogenesis. For example, adiponectin was detected to inhibit proliferation of A549 cells through its downregulation of binding to cyclic AMP response element (39). Another study described a fusion protein that could be used to improve antitumor activities in NSCLC by targeting epidermal growth factor receptor (EGFR) and insulin-like growth factor 1 receptor (IGFR1) (40). In the present study, KEGG pathway analysis indicated that the most enriched pathway was the AMPK signaling pathway. Phosphoinositide-3-kinase/Akt/mechanistic target of rapamycin kinase and Raf/Mek/extracellular signal-regulated kinase, which are regulated by AMPK, are crucial signaling pathways related to cell growth, survival and molecular expression (41). Furthermore, numerous studies have demonstrated that AMPK is an important mediator of NSCLC progression (42–44). Taken together, all of these results indicate that miR-9-3p participates in the tumorigenesis process and acts as an efficient regulator of NSCLC development.

As an important regulator of the cell cycle and cell proliferation, CCNA2 has been consistently demonstrated to be closely associated with the pathogenesis of lung cancer. For instance, coiled-coil domain containing 106 promoted proliferation of A549 and H1299 cells by upregulating CCNA2 expression (45). In another report, miR-137 and miR-30a inhibited tumor growth by inducing cell-cycle arrest and downregulating cell cycle-associated regulators, including CCNA2 (46,47). In addition, the gene E2F1 is an important transcription factor and is one of the most important pro-metastatic genes. Mei et al (48) observed that E2F1 had a significant effect on the prognosis of patients and was closely associated with patient survival, suggesting that E2F1 has an important role in the progression of clear cell renal cell carcinoma and hepatocellular carcinoma.

In the present study, the correlation analysis results suggested the existence of tight connections among LINC00968, miR-9-3p and CCNA2 in lung adenocarcinoma. Additionally, results of the Kaplan-Meier survival analysis indicated that LINC00968 may be a prospective biomarker for the diagnosis and prognosis of lung adenocarcinoma. Overall, these findings indicate that the LINC00968/miR-9-3p/CCNA2 regulatory axis deserves further research and investigation. The continuous development of novel research tools for identifying signal transduction pathways and discovering new therapeutic targets was inspired by the current investigation of molecular tumor profiles. Using tools to identify the prognostic and diagnostic markers in the present study presents an analysis concept based on existing tools, a number of which are available online, including R packages, Cytoscape, GO and the KEGG (49). Furthermore, the combination of the application of these tools and experimental verification will have a more important role in the identification of miRNAs. Therefore, the regulatory mechanisms of LINC00968 should be experimentally verified.

In conclusion, LINC00968 might have important functions in NSCLC via the LINC00968/miR-9-3p/CCNA2 regulatory axis. Results of the present study describe a novel approach to research the underlying mechanisms of LINC00968 in the tumorigenesis and development of NSCLC.

Acknowledgements

Not applicable.

Funding

The present study was supported by Guangxi Key Project of Science and Technology (grant no. 1598012-30), the Fund of the Guangxi Provincial Health Bureau Scientific Research Project (grant no. Z2013201), the Fund of the National Natural Science Foundation of China (grant nos. NSFC81660488, NSFC81360327 and NSFC81560469), and the Natural Science Foundation of Guangxi, China (grant no. 2015GXNSFCA139009). The funders had no role in the study design, the data collection and analysis, the decision to publish, or the preparation of the manuscript.

Availability of data and materials

The datasets generated and analyzed during the current study are available in GEO (https://www.ncbi.nlm.nih.gov/gds) and TCGA (http://cancergenome.nih.gov/).

Authors' contributions

SKL and GC conceived the idea for the study. DYL, WJC and JS analyzed and interpreted the data. DYL was a major contributor in writing the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that there are no competing interests.

References