Screening for key lncRNAs in the progression of gallbladder cancer using bioinformatics analyses.

The present study aimed to investigate key long non-coding RNAs (lncRNAs) and genes, and to obtain insights into their roles in the progression of gallbladder cancer (GBC). The gene expression profile and non‑coding RNA profile of GSE62335, which included five separate GBC tissue samples and five matched adjacent gallbladder normal tissue samples, was downloaded from the Gene Expression Omnibus database. The differentially expressed lncRNAs and mRNAs in the GBC tissues were identified, following which RNA binding protein analysis was performed using starBase v2.0 and the co‑expressed lncRNA‑mRNA pairs were predicted. Gene Ontology enrichment analysis for mRNAs was performed using the Database for Annotation Visualization and Integrated Discovery online tool. In addition, upstream microRNAs (miRNAs) were predicted for the co‑expressed lncRNAs and mRNAs. The results revealed that a total of 89 upregulated (13 lncRNAs and 76 mRNAs) and 261 downregulated transcripts (27 lncRNAs and 234 mRNAs) were identified in the GBC tissues. Only 9 lncRNAs had co‑expressed mRNAs, and lncRNA forkhead box P2 (FOXP2) was co‑expressed with the highest number of mRNAs, which were significant associated with the function of cell adhesion. In addition, the analysis of upstream miRNAs showed that FOXF1 adjacent non‑coding developmental regulatory RNA (FENDRR) had common upstream miRNAs, including miR‑18b‑5p, with another 119 differentially expressed genes, and that FENDRR was co‑expressed with adenomatosis polyposis coli downregulated 1 (APCDD1) and v‑kit Hardy‑Zuckerman 4 feline sarcoma viral oncogene homolog (KIT). Taken together, the results suggested that the lncRNAs FOXP2 and FENDRR may be crucial in promoting the progression of GBC via cell adhesion and regulating miR‑18b‑5p, or through interactions with KIT and APCDD1, respectively.

Introduction

Gallbladder cancer (GBC) is a rare, aggressive biliary tract malignancy, of which only ~20% of cases diagnosed in the early stages (1). The majority of patients with GBC present at an advanced stage, which is associated with a poor prognosis and poor survival outcomes for patients, with a 5-year survival rate ranging between 20 and 40% (2). Although considerable progress has been made in elucidating the genetic mechanisms involved in GBC (3,4), the molecular pathogenesis of GBC remains to be fully elucidated. Therefore, investigating novel molecules involved in the progression of GBC is like to facilitate understanding of the molecular mechanism underlying this disease and facilitate the development of novel therapeutic strategies.

Long non-coding RNAs (lncRNAs) are >200 nucleotides in length and the extent or levels of expression of lncRNAs generally appear to be specific to the tissues (5). Increasing evidence has demonstrated that lncRNAs are important regulators of diverse cellular processes and are emerging as being important in cancer (6). They have the potential to serve as therapeutic targets and prognostic indicators. For example, Wang et al (7) showed that the expression level of lncRNA-regulator of reprogramming (ROR) was upregulated in tissues from patients with GBC, and that the overexpression of lncRNA-ROR promoted the proliferation, migration and invasion of tumor cells, which was significantly associated with a poor outcome. Another study analyzed the expression of lncRNA actin filament associated protein 1 antisense RNA1 (AFAP1-AS1) using quantitative polymerase chain reaction analysis in GBC tissues and produced survival plots, which demonstrated that lncRNA AFAP1-AS1 was correlated with a poor prognosis in patients with GBC (8). Additionally, a competitive endogenous RNA (ceRNA) hypothesis has been suggested, in which lncRNAs are involved in cancer progression via interacting with microRNAs (miRNAs) (9,10). For example, a previous study reported that lncRNA Gall bladder cancer associated suppressor of pyruvate carboxylase, a target of miRNA (miR)-17-3p, negatively regulates pyruvate carboxylase-dependent cell proliferation in GBC (11). In addition, the findings of Wang et al (12) suggested that the lncRNA H19 may regulate the expression of forkhead box M1 (FOXM1) by competitively binding endogenous miR-342-3p in GBC. Although progress has been made in understanding the roles of lncRNAs in the progression of GBC, the molecular mechanisms underlying of lncRNAs require detailed investigations.

In the present study, a set of bioinformatics approaches were used to comprehensively analyze the publicly available microarray data from the Gene Expression Omnibus (GEO) database, which included five separate GBC tissue samples and fived matched adjacent gallbladder normal tissue samples. The differentially expressed lncRNAs and mRNAs were first identified in separate GBC tissues, and compared with those of matched adjacent normal gallbladder tissues. Subsequently, co-expressed lncRNA-mRNA pairs were identified, followed by functional enrichment analysis for these mRNAs. The shared upstream miRNAs regulating the co-expressed lncRNA and mRNAs were also predicted. The present study is aimed to examine key lncRNAs potentially involved in the progression of GBC and to elucidate their molecular mechanisms in the development of GBC.

Materials and methods

Affymetrix microarray data

The microarray data of GSE62335 was downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), which was deposited by Ma et al (13) on 15th October, 2014. In total, 10 samples were used to develop the array data, which included five separate GBC tissues and five matched adjacent normal gallbladder tissues. The raw data and annotation files were downloaded based on the platform of the GPL16686 Affymetrix Human Gene 2.0 ST Array (Affymetrix Inc., Santa Clara, CA, USA) for further analysis.

Data preprocessing and screening of differentially expressed lncRNAs and mRNAs

All the raw data were preprocessed using the robust multichip average method in the oligo package (available through Bioconductor version 3.0; http://www.bioconductor.org) (14). The comparison of differentially expressed lncRNAs and mRNAs in separate GBC tissues with normal tissues were screened using the Limma package version 3.22.7 of Bioconductor version 3.0 (http://www.bioconductor.org/packages/release/bioc/html/limma.html) (15). P<0.05 and |log2 fold-change (FC)|>0.58 were defined as the cut-off values for screening.

RNA binding protein (RBP) analysis

The starBase v2.0 database (http://starbase.sysu.edu.cn/) (16,17) is an experimentally supported database, which provides the most comprehensive protein-RNA, miRNA-mRNA and miRNA-lncRNA interactions supported by large-scale cross-linking immunoprecipitation (CLIP)-Seq (HITS-CLIP, PAR-CLIP, iCLIP and CLASH) data. The newly developed starBase v2.0 database has been used in several studies to identify between miRNAs and lncRNAs, which have also been validated by experiments (18,19). In the present study, RBPs directly interacting with differentially expressed lncRNAs and mRNAs were identified based on the information in starBase v2.0. The lncRNAs and mRNAs, which were able to bind to the same RBP, were then selected for subsequent analysis.

Co-expression network construction and functional enrichment analysis

The Pearson correlation coefficient (20,21) is widely used as a co-expression similarity measure. In the present study, the lncRNAs and mRNAs found to interact with the same RBP were selected for co-expression analysis, and the co-expressed lncRNA-mRNA pairs were screened out with the absolute value of a Pearson correlation coefficient >0.85. The co-expression network was then constructed with the significant lncRNA-mRNAs using Cytoscape software version 3.0.2 (http://www.cytoscape.org/) (22).

The Gene Ontology (GO) database (http://www. geneontology.org/) is widely used for the unification of biological functions and provides gene annotation terms for large-scale genomic or transcriptomic data (23,24). The Database for Annotation Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) (25) is used for systematically associating the functional terms with gene or protein lists. The DAVID Bioinformatics Resources consist of the DAVID Knowledgebase and a set of integrated, web-based functional annotation tool suites, which can be used to screen large gene lists in detail by comprehensively integrating several different biological angles to extract associated biological meanings (26). In the present study, the DAVID online tool was used for GO enrichment analysis of mRNAs in the co-expressed lncRNA-mRNA pairs identified. P<0.01 was considered the cut-off value.

Upstream miRNA analysis

Based on the information obtained on the miRNA-lncRNA and miRNA-mRNA pairs in starBase v2.0, the upstream miRNAs of lncRNAs and mRNAs were respectively predicted. The shared upstream miRNAs of the co-expressed lncRNAs and mRNAs were then identified according to the their expression levels.

Results

Screening of differentially expressed lncRNAs and mRNAs

Using the Limma package with P<0.05 and |log2 FC|>0.58 as cut-off thresholds, a total of 89 upregulated transcripts (13 lncRNAs and 76 mRNAs) and 261 downregulated transcripts (27 lncRNAs and 234 mRNAs) were identified in the separate GBC tissues, compared with the normal tissues. The heatmap of the differentially expressed transcripts is shown in Fig. 1.

Figure 1.

Heatmap of differentially expressed transcripts. The left vertical strips represent the distribution of transcripts, yellow strips represent mRNAs and green strips represent long non-coding RNAs. Red indicates upregulated transcripts and green indicates downregulated transcripts.

RBP analysis

A total of 2,759 RBP-RNA pairs were initially identified according to the information of starBase v2.0, among which 1,747 were screened out as they were able to interact with not only lncRNAs, but also mRNAs.

Co-expression and functional enrichment analysis

In the present study, a total of 287 co-expressed lncRNA-mRNA pairs were screened with the absolute Pearson correlation coefficient value of >0.85. Of these co-expressed lncRNA-mRNA pairs, 65 were involved in the 1,747 identified RBP-RNA pairs, from which a co-expression network was constructed (Fig. 2). Notably, all mRNAs within the 65 lncRNA-mRNA pairs were linked with only nine lncRNAs, including FOXP2, clusterin (CLU) and coxsackie virus and adenovirus receptor pseudogene 3 (CXADRP3), with co-expressed mRNAs (Table I). FOXP2 was also able to co-express with the highest number of mRNAs, followed by CLU.

Figure 2.

Co-expression network of lncRNAs and mRNAs. Red nodes represent upregulation and blue node represent downregulation. Diamond-shaped nodes represent mRNAs and circular nodes represent lncRNAs. LncRNAs, long non-coding RNAs.

Table I.

Number of co-expressed mRNAs corresponding to lncRNAs.

LncRNA	Co-expressed mRNAs (n)
AC011526.1	1
CLU	11
CXADRP3	4
FENDRR	2
FOXP2	21
LINC00578	10
MIR22HG	1
PTGER3	6
SCARNA22	9

CLU, clusterin; CXADRP3, coxsackie virus and adenovirus receptor pseudogene 3; LncRNAs, long non-coding RNAs; FOXP2, forkhead box P2; FENDRR, FOXF1 adjacent non-coding developmental regulatory RNA; LINC00578, long intergenic non-protein coding RNA 578; MIR22HG, MIR22 host gene; PTGER3, prostaglandin EP3 receptor; SCARNA22, small Cajal body-specific RNA 22.

As shown in Fig. 3, the co-expressed mRNAs connected with FOXP2 and CXADRP3 were significantly enriched in several GO terms, including cell-matrix adhesion, cell-substrate adhesion and cell adhesion. The co-expressed mRNAs, which had interactions with CLU were significantly associated with the regulation of cell proliferation and cell surface receptor signaling pathway terms.

Figure 3.

GO enrichment analysis of mRNAs corresponding to lncRNAs. Red strips represent enriched GO terms. Black strips represent non-enriched GO terms. The abscissa axis above represents lncRNAs. The longitudinal axis represents enriched GO terms. GO, Gene Ontology; LncRNAs, long non-coding RNAs.

According to the information from starBase v2.0, only lncRNA FOXF1 adjacent non-coding developmental regulatory RNA (FENDRR) was able to interact with upstream miRNAs. FENDRR shared common upstream miRNAs with 119 differentially expressed mRNAs. However, based on the correlation information in the starBase version 2.0 database, not all 119 differentially expressed mRNAs exhibited a correlation with FENDRR; thus, it was not possible to construct the ceRNA interaction network with these 119 differentially expressed mRNAs, FENDRR and the shared upstream miRNAs.

However, FENDRR was revealed to be co-expressed with adenomatosis polyposis coli downregulated 1 (APCDD1; correlation coefficient, 0.905) and v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog (KIT; correlation coefficient, 0.846). Notably, upstream miRNA analysis showed that miR-18b-5p, miR-28-5p and miR-4735-3p were common upstream miRNAs of FENDRR and KIT, whereas miR-18b-5p was the common upstream miRNA of FENDRR and APCDD1.

Discussion

In the present study, a comprehensive biological informatics approach was used to examine key lncRNAs involved in GBC and to elucidate their molecular mechanisms in the development of GBC. In total, 13 upregulated and 27 downregulated lncRNAs were identified in the separate GBC tissues, compared with the normal tissues. Among these lncRNAs, FOXP2 was found to be downregulated in the GBC tissues and co-expressed with the highest number of differentially expressed mRNAs, which were significantly enriched in GO terms associated with cell adhesion. FENDRR was also a downregulated lncRNA, which was found to share common upstream miRNAs with 119 differentially expressed mRNAs, including miR-18b-5p, and was co-expressed with KIT and APCDD1. These results suggested that lncRNA FOXP2 and FENDRR were likely to be important in the pathogenesis of GBC via cell adhesion and regulating miR-18b-5p, or the interaction with KIT and APCDD1, respectively.

Cell adhesion is an important event in biological processes, which determines the polarity of cells and is the morphological hallmark of malignant tumors (27). Cell-cell adhesiveness is usually reduced in several types of cancer, which is considered to be essential for cancer invasion and metastasis (27). Increasing evidence has highlighted the roles of cell adhesion molecules in perineural invasion in GBC (28). The overexpression of epithelial cell adhesion molecule antigen is an independent prognostic marker in GBC (29). In addition, E-cadherin is a transmembrane cell-adhesion protein mediating intercellular adhesion in epithelial tissues (30). It is reported that loss of the expression of E-cadherin is important in the progression from tumor formation to invasion and metastasis in several types of cancer (31). Kohya et al (32) confirmed that E-cadherin has a critical role in the proliferation, motility and invasion of GBC cells. Therefore, cell adhesion may be an important biological process involved in GBC. By contrast, lncRNAs have been shown to regulate gene expression through several mechanisms (33). LncRNAs have also been shown to allosterically modify RBPs, which regulate transcription (33). In the present study, lncRNA FOXP2 was found to be downregulated in separate GBC tissues, compared with that in matched adjacent normal gallbladder tissues, and the majority of mRNAs co-expressed with lncRNA FOXP2 were significantly associated with the function of cell adhesion. Therefore, it was hypothesized that lncRNA FOXP2 may be crucial in the progression of GBC via co-expression with important molecules in cell adhesion.

KIT is a member of the class III receptor tyrosine kinase family, which can be bound to and activated by stem cell factor (34). KIT receptor activation can cause the phosphorylation and activation of the Shc adaptor protein (35) and Ras (36). The Shc adaptor protein can link c-KIT to the Ras/mitogen-protein kinase (MAPK) pathway (37). K-ras gene mutation is frequently found in GBC cells (38). Activated Ras induces gefitinib-resistance by activating the epidermal growth factor receptor signaling pathways in human GBC cells (39). In addition, KIT be crucial in the majority of gastrointestinal stromal tumors via activation of the phosphatidyl-inositol-3-kinase (PI3-K)/AKT survival pathway. The activated MAPK and AKT pathways can lead to the uncontrolled growth of GBC epithelium, and promote the invasion and metastasis of GBC cells (40). By contrast, APCDD1 is a membrane-bound glycoprotein, which can interact in vitro with WNT3A and low-density lipoprotein receptor-related protein 5, which are two essential components of the Wnt signaling pathway (40). A previous study demonstrated that APCDD1 is the activated downstream target of the WNT/β-catenin pathway, which may enhance the proliferation of pediatric Wilms' tumor cells (41). Mutations of components of the Wnt/β-catenin pathway, including β-catenin exon 3 mutations, have been found to be associated with the progression of GBC (42). In addition, Wnt signaling can cross-talk with the E-cadherin cell adhesion system, which is involved in the morphogenesis of several types of cancer, including GBC (27,32). In the present study, FENDRR, as a downregulated lncRNA in GBC, was co-expressed with KIT and APCDD1. In a previous study, the expression profiles of lncRNAs in infantile hemangioma were assessed, and significant associations of FENDRR and its associated mRNAs were identified (43). Based on the results of the present study, and previous evidence confirming the importance of KIT and APCDD1 in the progression of GBC, FENDRR may be crucial in the progression of GBC through interaction with these two molecules by targeting several signaling pathways, including the Wnt signaling pathway.

In the present study, the analysis of upstream miRNAs also showed that miR-18b-5p was a common upstream miRNA of FENDRR, KIT and APCDD1. miR-18b-5p acts as an oncogene, and has been reported to be involved in the development and pathogenesis of human gastric cancer (44). The expression of miR-18b-5p is associated with the tumor protein 53 status, and may predict clinical outcome in patients with head and neck squamous cell carcinoma (45). Although the roles of miR-18b-5p in the development of GBC remain to be fully elucidated, the present study hypothesized that miR-18b-5p may be involved in the progression of GBC via targeting KIT and APCDD1, and lncRNA FENDRR may be involved in the development of GBC via regulating miR-18b-5p.

However, there were several limitations to the present study. Firstly, the size of samples was relatively small, therefore, further investigations with a higher number of samples are required. Secondly, statistical validation using other datasets, including meta-analysis (46,47) may be used to cross-check the results of the present study. Thirdly, due to insufficient material, experimental validation was not performed in the present study. The inclusion of these as a central component of future investigations can assist in confirming results using different approaches, including quantitative real-time polymerase chain reaction analysis.

In conclusion, lncRNA FOXP2 and FENDRR may be crucial in promoting the progression of GBC via cell adhesion and regulating miR-18b-5p or interaction with KIT and APCDD1, respectively. Therefore, these two factors may be recommended as potential targets in the management of GBC. The findings of the present study shed new light on the molecular mechanism of GBC and may provide a foundation for the development of credible therapeutic approaches. However, in the absence of experimental validation in the present study, further analysis of the microarray data and genetic experiments are required.