Integrated analysis of ceRNA network in lung adenocarcinoma based on bioinformatics analysis

Increasing studies have revealed that long noncoding RNA (lncRNAs) can bind to microRNA (miRNAs) as competitive endogenous RNA (ceRNA), which can affect the expression of mRNAs. These lncRNA related ceRNAs theoretically play an important role in the occurrence and progress of cancer. However, the roles and functions of ceRNA network in lung adenocarcinoma (LUAD) are not completely clear.


Introduction
Lung cancer has highly incidence rate and mortality rate in the world [1]. Non-small cell lung carcinoma (NSCLC) account for 85% of lung cancer cases [2]. Lung adenocarcinoma (LUAD) is the main histological subtype of NSCLC [2]. Despite the development in diagnosis and treatment in recent years, the prognosis of LUAD patients is still poor. Statistics has shown that the 5-year survival rate is less than 20% in LUAD [3]. Insu cient understanding of the biological mechanism of LUAD has limited the improvement of therapeutic e cacy.
LncRNAs are a class of single-stranded RNA molecules that have no protein-coding capacity. In the last few days, investigators discovered that LncRNAs can participate in many important biological processes, such as epigenetic regulation, transcriptional regulation, and post-transcriptional regulation, etc. [4].
Further studies have found that lncRNAs promote or inhibit the development of LUAD by participating in a variety of cytological processes including chromatin modi cation, gene expression regulation, transcriptional activation, transcriptional interference, intranuclear tra cking, and cell differentiation, affecting cancer metastasis and invasion, proliferation and apoptosis of cancer cells, cancer angiogenesis, and regulating cancer drug resistance, providing a new strategy for LUAD treatment [7][8][9].
Numerous studies have shown that lncRNAs can adsorb multiple miRNAs, and each miRNA has multiple downstream target genes combined, forming a quite complex regulatory network of lncRNAs that become ceRNAs, which affect and regulate the expression of mRNAs and target genes. These ceRNAs associated with lncRNAs are theoretically closely related to the development and progression of cancer [10].
Therefore, lncRNAs acting as ceRNAs have diverse biological functions that deserve further exploration.
Herein, this study predicted differentially expressed genes of lncRNAs and mRNAs in LUAD patients through the GEO database. Then we constructed a ceRNA network. The most differentially expressed genes were veri ed by experiments. These lncRNAs involved in the ceRNA network may become potential therapeutic targets or diagnostic biomarkers for LUAD.

Datasets Selection
Gene expression microarray datasets were downloaded from the GEO database. The inclusion criterias were including: (1) the studies were about lung adenocarcinoma; (2) tissue samples including cancer and corresponding adjacent or normal tissues; (3) the total number of samples was not less than 40; (4) the dataset must include lncRNAs and mRNAs.

Analysis of Differential Gene Expression in LUAD
Differentially expressed lncRNAs (DElncRNAs) and DEmRNAs were obtained by gene expression pro ling of the GEO2R analysis datasets with screening conditions of | LogFC | > 1.5 and adjusted P-value < 0.05.

Constructing PPI Networks
Protein-protein interaction (PPI) analysis is an important tool for discovering LUAD-related differentially expressed genes of different modules based on protein networks. PPI networks were constructed to identify key coding genes and important gene modules in LUAD compared to paracancerous tissues. We used the String database and visualized using Cytoscape v3.8.2 to construct PPI networks.

Functional Enrichment Analysis
Functional enrichment analysis, including KEGG pathway and GO(including molecular functions, MF; cellular components, CC and biological processes, BP), were performed for the lncRNAs and mRNAs from the regulatory network with DAVID (Database for Annotation, Visualization and Integration Discovery, https://david.ncifcrf.gov/). P-value < 0.05 was set as the threshold.
2.6 Cell culture and transfection A549 cells and MRC-5 cells were provided by Procell. They were cultured in RPMI-1640 (Sigma-Aldrich, Whitehouse Station, New Jersey) medium containing 10% fetal bovine serum (FBS, Gibco, Carlsbad, California). Cells in the logarithmic growth phase were subjected to transfection at a con uency of 70-80%. Transfection vectors and Lipofectamine 6000 (Beyotime Biotechnology, China) were respectively diluted in Opti-MEM. After mixture and incubation for 20 min, the cells were plated in each well.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR).
Total RNA was extracted from NSCLC A549 cell and lung cells of human MRC-5 using trizol reagent, and was quanti ed using NanoDrop 2000 spectropho tometer. Total RNA was reverse transcribed into complementary deoxyribonucleic acid (cDNA) using PrimeScript RT reagent. The obtained cDNA was further ampli ed by quantitative PCR using SYBR Premix Ex Taq. The relative levels were quantitatively analyzed using the 2-∆∆Ct method. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal reference. Cells were seeded in a 96-well plate with 2.0×10 3 cells/well and incubated with 0, 0.5, 1, 1.5, 2 and 2.5 mg/ml, respectively. Each experiment was performed in triplicate. Absorbance at 450nm was recorded using CCK-8 kit , the viability curve was plotted. Cell viability was measured at 1 day, 2 days, 3 days, 4 days and 5 days.

Cell migration assay
Cells were seeded in a 6-well plate with 5.0×10 4 cells/well and incubated. Mark drew a straight line in the 6-well plate. The dropped cells were washed with PBS. The scratch width of cells was observed at the time of 0h and 48h (n=5).

Cell invasion assay
For the invasion assay, 4.0×10 4 cells were added to the top of the Transwell membrane in the upper chamber coated with Matrigel. A complete medium was added to the bottom wells as stimulation of invasion. The cells were then incubated at 37°C for 48h. The cells on the upper surface were washed away, while the cells on the bottom surface were xed with 20% methanol (5 minutes) and stained with 0.1% crystal violet. Then stained cells were observed and counted under an inverted microscope, and ve different microscopic views were randomly selected for analysis, including different multiples, and experiments were performed in the triplicate assay (n=3).

Cell apoptosis assay
Apoptosis was analyzed by ow cytometry with Annexin V-PI apoptosis detection kit. Cells were treated with siRNA for 48h. Cells were then collected and stained with AnnexinV-PI according to the instructions. The early and late apoptotic cells were detected by ow cytometry after double staining.

Statistical analysis
Student's t-test and one-way ANOVA analysis were performed to analyze the data using SPSS 22.0 software. P < 0.05 was considered statistically signi cant.  The Venn diagram showed the intersection of DEmRNAs and DElncRNAs, with a total of 494 DEmRNAs and 6 DElncRNAs. (Figure 2).  3.3 GO and KEGG pathway analysis GO and KEGG pathway enrichment analyses revealed that the DEmRNAs involved in ceRNA network were remarkably associated with 37 BPs, including BMP signaling pathway, positive regulation of cell migration, and negative regulation of axon extension involved in axon guidance. The enrichment of MFs is mainly related to heparin-binding, signaling protein receptor binding, chemorepellent activity, SMAD binding, growth factor activity, receptor binding. We found that the most enriched CCs were heparinbinding, semaphorin receptor binding, and chemorepellent activity. KEGG enrichment analysis showed that DEmRNAs were signi cantly enriched in regulating stem cell pluripotency, TGF-β signaling pathway and protein digestion and absorption. In addition, we constructed a PPI network to identify the interaction between mRNA proteins in the ceRNA network. We found some genes with higher composite scores, including FGF2, COL1A2, BMP2, KLF4, CXCL5 and MMP11, were mainly enriched in "TGF-β signaling pathway", as shown in Figure 6.  The expression of lncRNA AFAP1-AS1 in A549 cells was 1.37 ± 0.15, which was signi cantly higher than that in MRC-5 by 0.57 ± 0.13 ( P < 0.05), as shown in Figure 6. The expression of lncRNA AFAP1-AS1 in siRNA-NC, siRNA-1, siRNA-2 and siRNA-3 group were 1.23 ± 0.11, 0.55 ± 0.09, 0.31 ± 0.04 and 0.44 ± 0.02 (F = 90.05, P < 0.001). Compared with siRNA-NC group, the expression of lncRNA AFAP1-AS1 decreased signi cantly in siRNA-2 group, as shown in Figure 6. The results of tumor growth in LUAD cells with altered expression of AFAP1-AS1 were shown in Figure 8. On the fourth and fth days, the growth rate of silencing AFAP1-AS1 was signi cantly lower than that of the control group and the difference was statistically signi cant. The effect of AFAP1-AS1 on the migration ability of LUAD cells was detected by cell migration assay. As shown in Figure 9, in A549 cell line, the number of cell migration in AFAP1-AS1 siRNA-2 group was signi cantly lower than that in the control group. The down-regulation of AFAP1-AS1 on the invasion of LUAD cells was detected by Transwell experiment. The number of cell invasion in down-regulation of AFAP1-AS1 was lower than that in the control group, as shown in Figure 10 (P < 0.05).  The results of apoptosis in LUAD cells with down-regulation of AFAP1-AS1 were shown in Figure 11. The suppression of AFAP1-AS1 signi cantly increased the apoptosis of A549 and the difference was statistically signi cant. Discussion Statistic showed that LUAD caused about 500,000 deaths each year all over the world [7]. Up to now, the exact pathogenesis of LUAD has remained incompletely understood. Therefore elucidating the molecular mechanisms underlying LUAD is of great importance for identifying new therapeutic targets to improve the clinical outcomes of patients. Accumulating studies indicated that lncRNAs could in uence the occurrence and progression of cancer through multiple pathways [11][12]. In 2011, Salmena et al. [13] rst proposed lncRNA-miRNA-mRNA interactions which provided a basis for the in-depth exploration of the mechanisms of various diseases, particularly cancer. Increasing studies were carried out based on this theory. Sun et al. [14] found that lncRNA H19 induced miR-675-5p up-regulation, while the level of P53 was signi cantly down-regulated, and p53 was identi ed as a target gene of miR-675. In addition, lncRNA LINC00355 and LINC00466 can target CCNE1 and HOXA10 by sponging miR-195 and miR-144, respectivelyto promote the progression of LUAD [15][16].
To the best of our knowledge, there are few ceRNA regarding LUAD. Thus, we obtained microarray datasets from the GEO database using bioinformatics analysis to mine LUAD-related lncRNAs and mRNAs with differential expression. By constructing the lncRNA-miRNA-mRNA ceRNA network, a total of 6 lncRNAs, 22 miRNAs and 55 mRNAs were identi ed, which included 83 nodes and 155 edges. The results of GO enrichment analysis indicated that RNAs associated with ceRNAs are mainly involved in cancer-related biological processes, for example, positive regulation of DNA transcription, signaling protein receptor binding and extracellular space. In addition, DEmRNAs involved in the ceRNA network were signi cantly abundant in three KEGG pathways (signaling pathway regulating stem cell pluripotency, TGF-β signaling pathway, protein digestion and absorption). Enrichment analysis suggests that by regulating these biological processes and pathways, LUAD-speci c ceRNA networks may be involved in cancer processes. Which suggested that the lncRNA-miRNA-mRNA ceRNA network may regulate the biological processes and pathways of LUAD.
LINC01207 was found to be a novel LUAD-speci c lncRNA that was signi cantly up-regulated in LUAD, and the expression level of LINC01207 correlated with TNM stage. Both in vivo and in vitro experiments have shown that LINC01207 can promote LUAD cell proliferation [17]. Acha-Sagredo et al. [18] found that LINC00968 was down-regulated in LUAD and could inhibit the progression of LUAD through the miR-21-5p/SMAD7 signaling axis. Unfortunately, there were few studies about the roles of PCAT19 and LINC00312 in LUAD. Therefore, further research is needed to clarify the role of these genes in LUAD.
The actin ber-associated protein 1-antisense RNA1 (AFAP1-AS1), which is located on the antisense strand of chromosome 4 of humans, plays a vital role in mammalian development [19][20]. LncRNA AFAP1-AS1 has been described as an oncogene in various cancers, including gastric, colorectal, liver, hepatocellular and lung cancers [21]. High expression of AFAP1-AS1 was found to be an independent risk factor for poor patient prognosis in LUAD and could promote LUAD cell apoptosis by regulating the miR-545-3p/hepatocarcinoma-derived growth factor axis [22]. To test the biological function of AFAP1-AS1, we designed siRNA that e ciently silenced AFAP1-AS1 expression in A549 cells. The expression level of AFAP1-AS1 in LUAD cells was signi cantly higher than that in normal lung cells. As predicted, the silence of AFAP1-AS1 inhibited proliferation, suppressed metastasis ability, reduced invasive ability and promoted apoptosis of LUAD cell line. These results indicated that AFAP1-AS1 might function as an oncogene and exhibit important role development and progression of lung cancer.