Integrative Analysis of ceRNA Network Reveals Functional lncRNAs in Intrahepatic Cholangiocarcinoma

Intrahepatic cholangiocarcinoma (ICC) is the second most common lethal liver cancer worldwide. Currently, despite the latest developments in genomics and transcriptomics for ICC in recent years, the molecular pathogenesis promoting ICC remains elusive, especially in regulatory mechanisms of long noncoding RNAs (lncRNAs), which acts as competing endogenous RNA (ceRNA). In order to elucidate the molecular mechanism of functional lncRNA, expression profiles of lncRNAs, microRNAs (miRNAs), and messenger RNAs (mRNAs) were obtained from The Cancer Genome Atlas (TCGA) database and an integrative analysis of the ICC-associated ceRNA network was performed. Moreover, gene oncology enrichment analyses for the genes in the ceRNA network were implemented and novel prognostic biomarker lncRNA molecules were identified. In total, 6,738 differentially expressed mRNAs (DEmRNAs), 2,768 lncRNAs (DElncRNAs), and 173 miRNAs (DEmiRNAs) were identified in tumor tissues and adjacent nontumor ICC tissues with the thresholds of adjusted P < 0.01 and |logFC| > 2. An ICC-specific ceRNA network was successfully constructed with 30 miRNAs, 16 lncRNAs, and 80 mRNAs. Gene oncology enrichment analyses revealed that they were associated with the adaptive immune response, T cell selection and positive regulation of GTPase activity categories. Among the ceRNA networks, DElncRNAs ARHGEF26-AS1 and MIAT were found to be hub genes in underexpressed and overexpressed networks, respectively. Notably, univariate Cox regression analysis indicated that DElncRNAs HULC significantly correlated with overall survival (OS) in ICC patients (P value < 0.05), and an additional survival analysis for HULC was reconfirmed in an independent ICC cohort from the Gene Expression Omnibus (GEO) database. These findings contribute to a more comprehensive understanding of the ICC-specific ceRNA network and provide novel strategies for subsequent functional studies of lncRNAs in ICC.


Introduction
Intrahepatic cholangiocarcinoma (ICC) is the second most common lethal liver cancer worldwide [1]. Various risk factors have been confirmed to contribute to the progression of this disease, including sclerosing cholangitis, chronic hepatitis B or C viral (HBV, HCV) infection, and fibropolycystic liver disease, and so on [2]. Surgical resection, chemotherapy, and emerging immunotherapy are alternative options to cure patients with ICC. However, due to the high rate of recurrence in this tumor [3], none of these approaches can significantly prolong long-term survival. Currently, although latest advances in genomics and transcriptomics [4][5][6] have broadened our knowledge of ICC tremendously, the molecular pathogenesis promoting ICC remains elusive [7]. erefore, there is a great need for understanding the specific molecular mechanism of tumors and for the identification of potential molecular biomarkers for ICC diagnosis and treatment.
In recent years, accumulating studies have focused on long noncoding RNAs (lncRNA) that are defined as transcripts over 200 nucleotides in length and have indicated that lncRNAs substantially affect gene expression that is dysregulated in numerous cancers. For instance, homeobox transcript antisense intergenic RNA (HOTAIR) was found to facilitate tumorigenesis by promoting phosphatase and tensin homolog (PTEN) methylation [8], and PVT1 binds EZH2 directly to silence ANGPTL4 expression by promoting cell growth and migration in cholangiocarcinoma [9]. One of the most known mechanisms for lncRNA is that it works as a competing endogenous RNA (ceRNA). is hypothesis was first proposed by Salmena et al. [10], which holds that lncRNA is involved in posttranscriptional regulation by functioning as sponges to modulate both miRNAs and mRNA expressional levels. In addition, emerging studies have confirmed that lncRNA acts as a hub gene in the regulation network between miRNA and target genes that are involved in a variety of cellular biological processes, especially in tumor proliferation and metastasis [11]. However, little is known about the comprehensive landscape of the ICC-associated ceRNA regulatory network, and very few studies were performed to investigate the lncRNAs mechanism for ICC.
In the present study, in order to elucidate the molecular mechanism of lncRNA and the lncRNA-mediated regulatory network, we performed an integrative analysis of identifying differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs), and miRNA (DEmiRNAs) of ICC using data from e Cancer Genome Atlas (TCGA) database. Next, an ICC-associated ceRNA network was successfully constructed and the overall survival (OS) analyses were carried out to identify molecules that were novel prognostic biomarkers and potential targets for ICC [12].
is study contributes to the exploration of the ceRNA regulatory mechanism of ICC and provides valuable insight for further functional research. e RNA and miRNA expression profiles and clinical follow-up data of the ICC cohort were downloaded from the TCGA database.

Materials and Methods
is study was conducted in accordance with the publication guidelines of TCGA (http://cancergenome.nih.gov/publications/public ationguidelines). In addition, the lncRNA expression profiles from GSE89749 (including 81 ICC samples) were used as an independent validation set.

Differentially Expressed Analysis.
e raw read counts (lncRNA, miRNA, and mRNA) were first normalized by using the trimmed mean of M-values (TMM) method in the "edgeR" package of R software (Version: 3.4.3). e expression levels of lncRNA, miRNA, and mRNAs in ICC were calculated as the log 2(x + 1) of the TMM normalized level. To identify DElncRNAs, DEmiRNAs, and DEmRNAs, we compared the nontumor tissues with the ICC tumors by using the "edgeR" package [13] with a significance threshold of an absolute log2 fold change ≥2 and a false discovery rate-(FDR-) adjusted P value less than 0.05. In addition, a volcano plot was made by using the ggplot2 package to identify the differential expression genes (DEGs) with statistical significance between the tumor and nontumor groups. Hierarchical clustering and heatmap plot analyses were performed for the DEGs of DElncRNAs, DEmiRNAs, and DEmRNAs by using ComplexHeatmap [14] package in R.

Constructing the ceRNA Network Analysis.
To investigate the potential interactions between lncRNAs and mRNAs, we constructed the ceRNA coregulated network using DElncRNAs, DEmiRNAs, and DEmRNAs. e coexpression analysis was performed by using Pearson's correlation coefficient (PCC) to identify correlation pairs according to the expression levels between significantly dysregulated DElncRNAs and DEmRNAs.
e Pearson correlation cutoff value >0.7 or <− 0.7 with P < 0.001 was defined as the screening threshold for retaining the RNAs in the analysis. Finally, the coregulated network was constructed using Cytoscape (version 3.7.0) [15] to illustrate the ceRNA network.

GO and KEGG Functional Enrichment Analysis of DEGs.
e GO [16] database stores extensive information of gene sets including GO terms and the annotations of genes and provides informative pathways for substantial genes. e GO function and KEGG pathway involved in DEGs were analyzed by the online tool DAVID [17] (version: 6.7, https://david-d.ncifcrf.gov/summary.jsp). e terms with the number of enriched genes counts ≥2 and the hypergeometric significance P values <0.05 were considered significant.

Protein-Protein Interaction (PPI) Network Submodule
Analysis. In the PPI network, proteins with similar functions tend to cluster. erefore, the analysis of the functional clustering module in the PPI network may help us understand the unknown functions of proteins. In this study, the Molecular Complex Detection (MCODE) plugin in Cytoscape [18], a recursive vertex weighting scoring scheme based algorithm, was used to analyze the subnetwork module in the PPI network with a module score >0.5. By using the MCODE plugin of Cytoscape (Version: 1.4.2), the submodule of PPI network was analyzed with the default thresholds: Degree Cutoff: 2, Node Score Cutoff: 0.2, K-Core: 2, and Max Depth: 100. GO and KEGG pathway enrichment analyses for the DEGs involved in the module were then performed.

Prognostic Analysis.
In order to determine the association between DEmRNAs, DElncRNAs, and DEmiRNAs in the ceRNA network and prognostic OS in ICC, univariate Cox proportional hazards regression method was performed to analyze the relationship between the DElncRNAs and OS at a significant level of 0.05.

Validation of Prognostic lncRNA in an Independent Data
Set. To assess the prognostic lncRNA in another independent data set, lncRNA expression profiling of 81 ICC patients, based on Illumina HumanHT-12 V4.0 Expression Beadchip, was downloaded from the NCBI database of Gene Expression Omnibus (GEO) accession GSE89749. Optimal cutoff value was achieved by the function of surv_cutpoint in "survminer" R package. Kaplan-Meier plot was used to depict the survival curves and P < 0.05 was considered as significant.

Statistical Analysis.
Wilcoxon test was used for comparison between the two groups with R (Version: 3.4.3) software. Statistical significance was set at P < 0.05.

Prediction of Target miRNAs of DElncRNAs and Target mRNAs of miRNAs.
Based on the hypothesis that miRNAs interact with the lncRNAs through MREs, data of miRNAs targeting lncRNAs were from the miRcode and starBase v3.0 database. A total of 273 miRNAs that putatively target 10,349 lncRNAs, including 476,957 DElncRNA-DEmiRNA interactions, were obtained. Subsequently, we obtained the intersection of lncRNA-miRNA interactions by screening DElncRNAs and DEmiRNAs. Finally, 42 miRNAs and 196 lncRNAs were identified including 2132 DElncRNA-DEmi-RNA interactions (Table S1). Moreover, 173 identified DEmiRNAs were mapped into the mircode, miRDB, miRanda, miRTarBase, and TargetScan databases to predict their target mRNAs.
en, the intersection pairs of the DEmRNAs to miRNAs predicted by at least 3 three databases were considered as target DEmRNAs of DEmiRNAs. We obtained 2068 DEmRNAs targeted by 43 DEmiRNAs with 4623 DEmRNAs-DEmiRNAs pairs (Table S2).

Establishment of the ceRNA Network in ICCs and
Functional Enrichment Analysis. To better understand the pivotal roles of the DElncRNA, DEmiRNA, and DEmRNA in ICC, the ceRNA network was constructed on the hypothesis that lncRNAs are involved in posttranscriptional regulation by working as sponges to directly modulate both miRNAs and mRNA expressional levels. Based on this hypothesis, a total of 122 lncRNA-miRNA-mRNA interaction pairs were identified in the proposed ceRNA network, including 30 miRNAs (23 upregulated and 7 downregulated,  Table S5).

PPI Network Construction and Pathway Enrichment
Analysis of Subnetworks. To better understand the hub genes in ICC, we established a PPI network that harbored 78 genes according to the information from the STRING database v11 with scores of >0.4 (Figure 4). From Figure 4, we determined that the key hub nodes were RAC2, SYP, ERBB4, GNAO1, TNFRSF1B, CCR7, and CD4.

Discussion
Intrahepatic cholangiocarcinoma (ICC) is a serious disease affecting millions of individuals worldwide; however, the molecular pathogenesis promoting ICC remains unexplained. To date, numerous studies have been performed to reveal the underlying mechanisms of lncRNAs in regulating the development and progression of ICC, but these studies only focus on a single genetic event and may not satisfy the current requirement on discovering potential molecular targets for treating ICC. To better understand the role of lncRNA on the ceRNA network of ICC, we systematically integrated transcriptome expression profiles from TCGA to construct a lncRNA-associated ceRNA network to determine the regulatory mechanisms of lncRNAs. Moreover, we screened hub genes to identify novel prognostic biomarkers and potential treatment targets in the ceRNA network for ICC. In total, 6,738 DEmRNAs, 2,768 DElncRNAs, and 173 DEmiRNAs were identified in ICC tumors compared to nontumor tissues. We further constructed an ICC-specific network, and in order to make the results more reliable and stable, we used miRcode and starBase databases to predict target lncRNAs of DEmiRNAs and used mircode, miRDB, miRanda, miRTarBase, and TargetScan databases to predict target mRNAs of DEmiRNAs. We constructed an ICCspecific ceRNA network based on a total of 122 lncRNA-miRNA-mRNA interaction pairs, including 30 miRNAs, 16 lncRNAs, and 80 mRNAs. Pathway enrichment analysis was performed on the DEmRNAs in the ceRNA network and revealed that most of the enriched GO pathways were associated with the adaptive immune response, T cell selection, and positive regulation of GTPase activity categories.
Subsequently, 2 underexpressed and overexpressed networks were reconstructed to investigate the effects of lncRNA-mediated regulation on the ceRNA network. Notably, 2 hub regulators were highlighted in the ICC-specific ceRNA network, one being ARHGEF26-AS1, which interacted with 7 DEmiRNAs and was coregulated with 30 DEmRNAs. ARHGEF26-AS1 (ARHGEF26 antisense RNA 1) is an antisense RNA of ARHGEF26, and it has been reported that ARHGEF26-AS1 was significantly correlated with OS (P value < 0.05) in colon cancer [19]. However, no experimental evidence was found to support the biological function of ARHGEF26-AS1 in recent studies. MIAT was another important hub regulator, which interacted with 14 DEmiRNAs and was coregulated with 43 DEmRNAs. Many studies revealed that MIAT has been demonstrated as a key switcher in regulation of various biological and pathological processes in human cancers including non-small-cell lung cancer, neuroendocrine prostate cancer, and colorectal cancer [20][21][22]. e aberrant upregulation expression of MIAT in human malignancies can promote tumorigenesis through lncRNA-miRNA-mRNA networks. For example, it  was reported that MIAT regulates cell proliferation, migration, and invasion in gastric cancer through a mechanism involving the miR-29a-3p/HDAC4 axis [23]. Moreover, a PPI network harboring 78 hub genes was constructed, and according to the MCODE analysis, we finally obtained 4 modules consisting of 42 nodes and 120 edges ( Figure S1). Enrichment analysis for the subnetworks discovered that subnetwork 1 was strongly associated with the chemokine signaling pathway, cell receptor signaling pathway, and natural killer cell-mediated cytotoxicity. Subnetwork 2 was involved in renin secretion and GABAergic synapse. Subnetwork 3 was involved in the biosynthesis of amino acids, metabolic pathways, and carbon metabolism. Except for Subnetwork 4, no significantly associated KEGG pathways were observed. ese results suggest that lncRNA-associated ceRNA subnetworks may disrupt cancer immune and metabolic pathways in ICC.
To further investigate OS for these hub genes of DElncRNAs, DEmiRNAs, and DEmRNAs in the ceRNA network, Cox's proportional hazards analysis revealed that DElncRNA HULC was associated with the prognosis of patients with ICC (P value < 0.05). lncRNA HULC was found to have a lower expression level in tumors compared to nontumor tissues but the high expression of HULC in tumors was associated with poor prognosis of ICC. HULC functions as an oncogene in human cancer which promotes tumorigenesis by regulating multiple signaling pathways [24]. Especially in hepatocellular carcinoma (HCC), accumulating studies have been reported that lncRNA HULC has a potential as a promising therapeutic target in HCC and it can accelerate liver cancer by inhibiting PTEN via autophagy cooperation to miR15a [25][26][27]. However, to date, no studies have been reported about the association of HULC with ICC. With this study, we have shown that the aberrant upregulated expression of HULC affects OS in ICC and indicates a potential prognostic biomarker in ICC.
Our findings indicate that the lncRNA-associated ceRNA network plays pivotal roles in tumorigenesis and progression of ICC. Although our study provides comprehensive analysis of the competitive ceRNAs network and functional lncRNAs of ICC, several limitations of our research need to be noted. Firstly, the relationship between the ceRNA network and clinical features of ICC was not well addressed in our study due to incomplete clinical information. Secondly, the cohorts of ICC in the public databases were not large enough because ICC itself is a rare tumor and only accounts for about 3% of all gastrointestinal cancers [28], which might bias the results. irdly, additional experiments are needed to verify our results on whether ceRNA networks and functional lncRNAs are functional or specific for ICC rather than other cancer types.

Conclusions
Our study has identified several novel prognostic factors and potential treatment targets for ICC by analyzing RNA expression profiling data from the TCGA database. Our data provide a comprehensive understanding of the lncRNA-related ceRNA network in ICC but further functional experiments are needed to elucidate the underlying mechanisms.

Data Availability
e data used to support the findings of this study are included within the article. Moreover, RNA and miRNA

Conflicts of Interest
ere are no conflicts of interest to disclose. Figure S1: the ceRNA subnetworks of protein-protein interaction (PPI) network. Figure S2: Kaplan-Meier survival curves for lncRNAs HULC associated with overall survival of 81 ICC patients from GSE89749. Table S1: details of DElncRNAs-DEmiRNAs pairs. Table S2: details of DEmRNAs-DEmiRNAs pairs. Table S3: details of upregulated and downregulated miRNAs in ceRNA network. Table  S4: details of upregulated and downregulated lncRNAs in ceRNA network. Table S5: details of upregulated and downregulated mRNAs in ceRNA network. Table S6: GO pathway enrichment of mRNAs in ceRNA network.