Genome-wide identification of a competing endogenous RNA network in cholangiocarcinoma.

Cholangiocarcinoma (CCA) is the second widespread liver tumor with relatively poor survival. Increasing evidence in recent studies showed long noncoding RNAs (lncRNAs) exert a crucial impact on the development and progression of CCA based on the mechanism of competing endogenous RNAs (ceRNAs). However, functional roles and regulatory mechanisms of lncRNA-regulated ceRNA in CCA, are only partially understood. The expression profile of messenger RNAs (mRNAs), lncRNAs, and microRNAs (miRNAs) downloaded from The Cancer Genome Atlas were comprehensively investigated. Differential expression of these three types of RNA between CCA and corresponding precancerous tissues were screened out for further analysis. On the basis of interactive information generated from miRDB, miRTarBase, TargetScan, and miRcode public databases, we then constructed an mRNA-miRNA-lncRNA regulatory network. Kyoto Encyclopedia of Genes and Genomes and Gene Ontology analyses were conducted to identify the biological function of the ceRNA network involved in CCA. As a result, 2883 mRNAs, 136 miRNAs, and 993 lncRNAs were screened out as differentially expressed RNAs in CCA. In addition, a ceRNA network in CCA was constructed, composing of 50 up and 27 downregulated lncRNAs, 14 up and 7 downregulated miRNAs, 29 up and 25 downregulated mRNAs. Finally, gene set enrichment and pathway analysis indicated our CCA-specific ceRNA network was related with cancer-related pathway and molecular function. In conclusion, our research identified a novel lncRNA-related ceRNA network in CCA, which might act as a potential therapeutic target for patients with CCA.

INTRODUCTION

Although regarded as a rare type of cancer, cholangiocarcinoma (CCA), originating from biliary epithelium, globally accounts for 10% to 20% of all primary liver malignancies. There has been increasing interest in this tumor because of marked growth of the incidence and mortality rates for CCA in recent years.1 Since most patients with CCA are clinically asymptomatic, it is difficult for early diagnosis.2 Although several diagnosis methods including serological tests, such as serum CA19‐9 concentration, as well as radiological techniques, including magnetic resonance imaging and computed tomography, were applied for CCA, these modalities seem not always reliable and the sensitivity and specificity of these biological marker for the diagnosis of CCA is always low.3 Because majority of diagnosis appear at an advanced stage, the 5‐year survival is in general as low as 10%.4 In spite of chemotherapy acting as the standard treatment for patients with locally advanced or metastatic diseases, the effect of these chemotherapies including cisplatin and gemcitabine is limited.5 Consequently, it is urgent to detect novel biomarkers contributing to promising therapeutic strategies against CCA.

Long noncoding RNAs (lncRNAs) belonging to noncoding functional RNAs have drawn an increasing number of academic attention in recent years.6 However, lncRNAs were considered as “dark matter” across the whole genomic region in the beginning. Due to the appearance of high‐throughput sequencing technology, thousands of lncRNAs were identified recently. Actually, lncRNAs play a crucial role in diverse biological processes, especially tumorigenesis.7, 8 In recent studies, the roles of competing endogenous RNAs (ceRNAs) were inferred, based on which lncRNAs regulated various biological processes on its protein levels.9, 10 The core conception is ceRNAs regulate the large‐scale gene transcription through interaction with microRNA (miRNA) response elements in various cancers.11 Recently, several studies showed dysregulation of lncRNAs might be associated with the worse prognosis of CCA on the basis of the mechanism of ceRNAs.12 It was demonstrated that functionally crucial lncRNAs, for example, KCNQ1OT1, TUG1, H19, and HULC could play crucial impacts in the carcinogenesis of CCA based on various ceRNA mechanisms.13 Nevertheless, these studies identified potential mechanisms only based on a specific lncRNA‐miRNA‐mRNA axis. Consequently, it is urgent to identify the comprehensively competition effect of a ceRNA network in CCA, which could act as a potential therapeutic target for patients with CCA.

To address these challenges, we attempted to elucidate the specific lncRNAs of CCA involving in the ceRNA network in our study. By bioinformatic tools, a CCA‐related ceRNA network was constructed and gene set enrichment and pathway analyses were conducted to elucidate the potential regulatory roles of lncRNAs involved in this network, which might act as a promising therapeutic target for patients with CCA.

MATERIALS AND METHODS

Data processing and differential expression analysis

RNA‐seq profiles of patients with CCA were obtained from The Cancer Genome Atlas (TCGA) data set (https://cancergenome.nih.gov/), including 36 CCA tissues and 9 corresponding precancerous tissues. On the basis of genetic annotation downloaded from the Ensembl database (http://asia.ensembl.org/index.html), 19 667 protein‐coding genes and 14 367 lncRNAs were identified. Mature miRNAs expression profiles instead of pre‐miRNAs were obtained from UCSC Xena (http://xena.ucsc.edu/) and 1870 miRNAs were identified. Because all data profiles were downloaded from TCGA, no additional approval from the Ethics Committee were needed. Moreover, we followed the TCGA publication guidelines and data access policies (http://cancergenome.nih.gov/pulications/publicationguidelines) to perform all analysis in our study. Further, the differentially expressed mRNAs, lncRNAs, and miRNAs were identified separately based on the DESeq2 R package based on the subsequent criteria: |log2 fold change| ≥2 and P < .01.14

ceRNA network construction

The differentially expressed miRNAs, lncRNAs, and mRNAs were utilized to construct the ceRNA network. To maximize data reliability, interactions between mRNAs and miRNAs, which were validated by all of the three databases, miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/), miRDB (http://www.mirdb.org/), and Targetscan (http://www.targetscan.org/), were accepted for further analysis. In addition, the miRNA‐lncRNA interactive prediction was downloaded from the miRcode (http://www.mircode.org/) database. Taking not only lncRNA‐miRNA but also miRNA‐target gene interactive prediction into account, a CCA‐related ceRNA network on the basis of mRNA‐miRNA‐lncRNA axis was constructed. The ceRNA network was visualized by the Cytoscape software.15

Gene set enrichment and pathway analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted by the clusterProfiler R package to identify the biological function of the ceRNA network.16 The outcomes of the GO and KEGG analyses were visualized by the package of the GOplot R software.17

Survival analysis

The clinical information of patients with CCA was also downloaded from TCGA. The median expression level of each lncRNA acted as the cutoff to separate patients into the high or low expression group. Univariate Cox proportional hazards regression analysis was conducted to identify candidate genes involving in the ceRNA network significantly associated with overall survival. Kaplan‐Meier analysis was utilized to paint the survival curves and it was compared by the log‐rank test based on the threshold of P < .05.18

Validation of prognostic lncRNA in an independent data set

To evaluate the prognostic lncRNA originated from TCGA cohort, gene expression profiling (GSE89748) of CCA based on Illumina HumanHT‐12 V4.0 Expression Beadchip, including 72 patients, was downloaded from the NCBI Gene Expression Omnibus database.19 Optimal cut‐off was ascertain by the surv_cutpoint function of the survminer R package. Kaplan‐Meier analysis was utilized to paint the survival curves and it was compared by the log‐rank test based on the threshold of P < .05.18

Statistical processing

Statistical analyses were conducted by the R software (3.3.5 version). Continuous variables were examined by Student's t test and categorical variables by χ 2.

RESULTS

Differentially expressed mRNAs, miRNAs, and lncRNAs in CCA

Differential expression analysis was performed between 36 CCA and 9 corresponding precancerous tissues. A total of 2883 mRNAs (Figure 1A,B) 993 lncRNAs (Figure 2A,B; Table S1) were screened out as differentially expressed RNAs in CCA. Among them, 1736 upregulated mRNAs and 1147 downregulated mRNAs were screened out; 639 lncRNAs were overexpressed and 354 lncRNAs were low‐expressed. In addition, in terms of the miRNA, 136 miRNAs reached the inclusion criteria between CCA and corresponding precancerous samples, containing 54 low‐ and 82 overexpressed miRNAs (Figure 3).

Figure 1

Identification of differentially expressed mRNAs. A, The volcano plot showed that a total of 1736 overexpressed mRNAs and 1147 low‐expressed mRNAs (tumor vs normal tissues) were screened out. B, Circos plot of genome‐wide differentially expressed mRNAs (green represents low‐expressed mRNAs and red represents overexpressed mRNAs). mRNAs, messenger RNAs

Figure 2

Identification of differentially expressed lncRNAs. A, The volcano plot showed that a total of 639 over‐ lncRNAs and 354 low‐expressed lncRNAs (tumor vs normal tissues) were screened out. B, Circos plot of genome‐wide differentially expressed lncRNAs (green represents low‐expressed lncRNAs and red represents overexpressed lncRNAs). lncRNAs, long noncoding RNAs

Figure 3

Identification of differentially expressed miRNAs. The volcano plot showed that a total of 82 over‐ and 54 low‐expressed lncRNAs (tumor vs normal tissues) were screened out. lncRNAs, long noncoding RNAs; miRNAs, microRNAs

Construction of the ceRNA network in CCA

Taking not only miRNA‐lncRNA but also miRNA‐target gene interactivity into account, an integrated mRNA‐miRNA‐lncRNA network was constructed (Figure 4), including 50 overexpressed and 27 low‐expressed lncRNAs (Figure 5A), 14 overexpressed and 7 low‐expressed miRNAs (Figure 5B), 29 upregulated and 25 downregulated mRNAs (Figure 5C). This ceRNA regulatory network contained some well‐described biomarkers, for example, CBX2, HULC, and miRNA‐424. However, most molecules involving in this ceRNA network were required to be further investigated to ascertain the potential biological mechanism in CCA.

Figure 4

Competing endogenous RNA (ceRNA) network in CCA. The red indicates the risky RNAs, and the blue indicates the protective RNAs in CCA. The diamond represents lncRNAs, triangle represents miRNAs, and circle represents mRNAs. CCA, cholangiocarcinoma; lncRNAs, long noncoding RNAs; miRNAs, microRNAs

Figure 5

Heatmap of differentially expressed RNAs involving in ceRNA network. A, lncRNAs; (B) miRNAs; (C) mRNAs. lncRNAs, long noncoding RNAs; miRNAs, microRNAs; mRNAs, messenger RNAs

Functional enrichment analysis

On the basis of functional enrichment analysis, the biomarkers involving in this identified ceRNA network were remarkably correlated with five KEGG pathways, including “Central carbon metabolism in cancer”, “PI3K−Akt signaling pathway,” and “Cellular senescence” (Figure 6B). In terms of the GO categories, ceRNA network–related genes were significantly enriched in a number of molecular functional processes, such as “transcription coactivator binding,” “core promoter binding,” and “beta−catenin binding” (Figure 6A).

Figure 6

A, Gene Ontology (GO) terms of genes included in the ceRNA network. B, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of genes included in the ceRNA network

Survival analysis of lncRNAs, miRNAs, and mRNAs involving in the ceRNA network

To identify candidate ceRNA‐based biomarkers correlated with overall survival in CCA, lncRNAs, miRNAs, and mRNAs involving in the ceRNA network were explored in 36 patients with CCA. Therefore, univariate survival analysis was conducted to detect prognostic RNAs. As a result, due to the small sample size, only lncRNA HULC was identified to be prognosis associated in the TCGA cohort (Figure 7). In addition, lncRNA HULC in CCA was also validated to be significantly associated with overall survival in an independent cohort (GSE89748; Figure S1).

Figure 7

Kaplan‐Meier survival curves for lncRNA HULC associated with overall survival. lncRNAs, long noncoding RNAs

DISCUSSION

CCA, originating from biliary epithelium, contributes to a large number of deaths every year all around the world due to a high risk of mortality and morbidity. Because of high rates of recurrence and metastasis during the treatment of CCA,3 it is urgent to identify potential biomarkers of CCA to improve prognosis in patients. In recent studies, increasing evidence demonstrated that lncRNA plays a crucial biological role in the initiation, progression, metastasis, and recurrence of tumor.20, 21, 22, 23 Moreover, lncRNAs could also serve as the promising biomarkers for prognosis and diagnosis in various tumors.24, 25, 26, 27 Increasing research has showed that dysregulation of lncRNAs were correlated with the carcinogenesis of CCA, as well as decreasing the survival rate of patients.28, 29 Therefore, it is necessary to identify potential lncRNAs acting as promising prognostic biomarkers, which could function as a potential therapeutic target for patients with CCA.30, 31 Recently, it was reported that lncRNAs, also called miRNA sponges, could depress the activity of miRNAs by binding and sequestering them.32, 33 At the same time, lncRNAs could also regulate the expression of specific mRNAs. On the basis of the ceRNA hypothesis that lncRNA could act as ceRNA to interact with mRNA through competing with the corresponding miRNA.34 Recent experimental evidence validated that lncRNA, which possesses sequences similar to its targeted miRNA, could sequester these miRNA away from corresponding mRNA. Previous studies indicated that lncRNA KCNQ1OT1 was overexpressed in CCA tissues compared with paracancerous tissues and its positive role in cell proliferation, invasion, and the epithelial‐mesenchymal transition was also detected. KCNQ1OT1 could function as a ceRNA to promote carcinogenesis of CCA through mediating miR‐140‐5p/SOX4 signaling pathway.35 Zeng et al13 demonstrated lncRNA TUG1, acting as a ceRNA, could sponge miR‐145 to improve cancer progression, which could provide a new therapeutic strategy for treating CCA. In addition, lncRNA HULC and lncRNA H19, promoting cellular migration and invasion of CCA, could target CXCR4 and IL‐6 via ceRNA patterns of sponging miR‐372/miR‐373 and let‐7a/let‐7b, respectively.36

However, a whole genomic wide analysis of CCA‐related miRNA and lncRNA, especially based on high‐through sequencing, was still lacking. Therefore, it was urgent to identify the underlying regulatory mechanisms of the mRNA‐miRNA‐lncRNA ceRNA network in the progression of CCA. In the present study, a total of 2883 mRNAs, 136 miRNAs, and 993 lncRNAs were screened out as differentially expressed RNAs based on the TCGA database. Then, a CCA‐related ceRNA network was constructed according to specific criteria. This lncRNA‐miRNA‐mRNA ceRNA network was composed of 50 over‐ and 27 low‐expressed lncRNAs, 14 over‐ and 7 low‐expressed miRNAs, and 29 over‐ and 25 low‐expressed mRNAs. Further, GO and KEGG pathway analysis were utilized to identify predominant biological themes of differentially expressed genes involving the ceRNA network. Our results of GO analysis indicated that nine GO items was significant with P < .05. These significant GO terms involved “steroid metabolic process,” “transcription coactivator binding,” “neurotrophin receptor binding,” “core promoter binding,” “RNA polymerase II proximal promoter sequence−specific DNA binding”, “phosphatidylinositol 3−kinase activity,” “beta−catenin binding,” “core promoter sequence−specific DNA binding,” and “steroid binding.” The pathway analysis further indicated that five pathways were enriched, including “PI3K−Akt signaling pathway,” “Proteoglycans in cancer,” “Central carbon metabolism in cancer,” “Prolactin signaling pathway,” and “Cellular senescence.” It suggested that these biological processes and pathways were crucial factors for the ceRNA network to promote progression of CCA. Finally, through survival analysis, lncRNA HULC was identified as a prognostic biomarker among genes involving the ceRNA network, which could guide individualized therapy strategies.

Among these lncRNAs involving in the ceRNA network, Li et al37 have reported that Linc00483 acting as ceRNA regulated proliferation and apoptosis via activating mitogen‐activated protein kinases in gastric cancer. In addition, recent studies indicated that lncRNA UCA1 played important roles in the progression of various cancers, including papillary thyroid carcinoma, melanoma, esophageal squamous cell carcinoma, and breast cancer.38, 39, 40, 41 In the previous study, HULC could act as an oncogene through negatively regulating miR‐125a‐3p in ovarian carcinoma.42 Importantly, lncRNA HULC, which was activated by oxidative stress, could enhance CCA cellular migration and invasion by a ceRNA mechanism.36 However, the roles of all of these lncRNAs, except HULC, were never reported in CCA. Therefore, it is worthwhile to further explore the underlying mechanisms of these ceRNA‐network‐related lncRNAs in CCA.

Although our research constructed the ceRNA network in CCA for the first time, a few potential limitations also existed. First, although we have exert our utmost effort to mine the data resource, merely 36 CCA samples were obtained from TCGA, which led to selection bias. Second, our research was merely based on the samples of CCA, other vital characteristics, such as the status of IDH1 mutation and history of intrahepatic duct calculus were not taken into consideration. Third, the biologic potential mechanisms of these lncRNAs involving in ceRNA network in CCA remains to be further validated in the future. Finally, because this research was merely based on TCGA from the USA, prospective research in different populations and center containing larger patient sizes are needed to further validate our result.

In conclusion, the lncRNA‐miRNA‐mRNA ceRNA network were constructed for the first time in CCA and the underlying mechanisms involving in the progression of CCA were analyzed. Our results indicated that lncRNAs played a crucial role in the progression of CCA. LncRNA HULC was identified as a prognostic biomarker, which could contribute to individualized therapy improvement. Experimental studies in the future are needed to validate the potential biological mechanisms of these ceRNA‐network‐related lncRNAs in CCA.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

Supplementary information

Click here for additional data file.

Supplementary information

Click here for additional data file.

Supplementary information

Click here for additional data file.