Integrative Analysis of Transcriptome and MicroRNome Profiles in Cholangiocarcinoma

Aim: To identify deregulated expression of miRNAs and target genes in cholangiocarcinoma (CCA), as well as drug-gene interactions that may be useful for patient treatment. Methods: We analyzed quantitative transcriptome and miRNA deep sequencing expression data from 45 samples, 36 CCA and 9 histologically normal biliary tissues, obtained from the public repository The Cancer Genome Atlas (TCGA). Bioinformatic methods were used to identify and integrate differential expression of miRNA and transcriptome profiles in CCA vs. normal tissues. Deregulated miRNA and corresponding target genes were identified and mapped into miRNA-mRNA networks. Results: Results showed 64 differentially expressed miRNAs (48 over and 16 under-expressed) between CCA and corresponding normal biliary tissues. Additionally, 432 genes (180 over and 252 under-expressed) were identified between CCA and normal samples. We identified individual miRNAs with the largest number of gene targets. Among these, miR-125a was over-expressed and had the highest number of direct interactions with 33 mRNA targets. miRNAs miR-122 and miR-139 were the under-expressed miRNAs with the highest number of interactions (9 targets each). miR-122 was found to modulate the expression of the transcription factor FOXM1, known to be involved in tumorigenesis and the matrix metalloproteinase MMP7, an important mediator of tumor invasion. miR-148 and miR-194 were predicted to modulate NQO1, which is up-regulated in cancer and associated with treatment resistance in cholangiocarcinoma. Conclusion: The novelty of our study is the identification of complex deregulated networks of miRNAs and target genes in CCA. miRNAs with a large number of targets may have a higher functional impact on cell regulation. These findings contribute for a better understanding of CCA biology. Identified miRNAs and target genes are potential candidates for the design of validation strategies towards the characterization of clinically applicable biomarkers; such biomarkers may be useful for the development of molecularly-targeted therapeutics that can benefit patients.


Introduction
Cholangiocarcinoma (CCA) is a malignant tumor of the liver that arises from epithelial cells of the bile ducts, both within (~10%) and outside (~90%) of the liver [1,2]. Although rare, it is an aggressive tumor associated with a poor prognosis. Survival rates around 15% at 2 years after diagnosis have been reported in the USA [2]. A recent epidemiology study in the USA showed that mortality rates by CCA have increased significantly in the past decade: for individuals >25 years old, there was a 36% increase in mortality between 1999-2014: 2.2 to 3 per 100,000, and deaths were the highest among African Americans (45%), followed by Asians (22%) and whites (20%) [3]. In Brazil, 8,700 deaths by CCA combined with all liver cancers were reported in 2013 [4]. In Asia, the Japanese population has relatively higher rates of CCA, with an incidence of 3.5 and 3.05 cases per 100,000 in Osaka and Hiroshima, respectively. Other countries such as Singapore, Philippines and Vietnam have a lower incidence of 1.45, 1.2 and 0.1 respectively, per 100,000 people [5]. The incidence and mortality by intrahepatic CCA in Japan are similar to those other countries, however, incidence of extrahepatic CCA is higher [6].
Risk factors for development of CCA include chronic inflammation of the biliary tree, seen in primary sclerosing cholangitis, congenital fibropolycystic diseases of the biliary tract, HCV infection, previous exposure to Thorotrast, used in radiography of the biliary tree and infection of the biliary tracked by Clonorchis sinensis, common in southeast Asia [2]. Cholangiocarcinoma is often difficult to diagnose [7] and treat with conventional therapies due to the high intra and intertumoral genetic heterogeneity [8]. Surgical treatment, when possible, remains the best option to achieve better survival. Therefore, there is an urgent need for the development of genomic-driven treatment strategies for patients with this aggressive form of cancer.
MicroRNAs (miRNAs) are small RNAs of 18-22 nucleotides in length, non-protein coding and are potent gene expression regulators. A canonical mode of action of miRNAs is through binding to 3'-UTR of mRNA, leading to translational inhibition or mRNA degradation [9]. To date, 2,588 mature miRNA sequences have been identified in the human genome according to miRBase (http://www.mirbase. org/cgi-bin/mirna_summary.pl?org=hsa) [10]. More recently, a deep sequencing analysis of 13 different cell types identified 3,707 novel mature miRNAs, suggesting that there may be several more miRNA sequences still unknown in the human genome [11].
Expression and functional assessment of the role of miRNAs in CCA have been shown in previous studies. However, most studies, as reviewed in [24], have been mainly focused on candidate miRNAs rather than performing a large comprehensive miRNA expression analysis. Therefore, there is still a need to better understand the regulatory pathways involved in the pathogenesis of this aggressive cancer.
A recent study by Farshidfar et al. reported The Cancer Genome Atlas (TCGA) biomarker analysis of CCA through integrative analysis of copy number, gene expression, methylation and somatic mutations. This study identified distinct histological subtypes enriched for the presence of isocitrate dehydrogenase (IDH) mutations and compared global genomic changes occuring in cancers of the liver, pancreas and the biliary tree [25]. Our study used this same cholangiocarcinoma -TCGA dataset however it was focused on the role of global miRNA expression changes in transcriptome regulation. Our integrative analysis of transcriptome and miRNA expression data contributes to better understand a number of functionally relevant pathways in CCA. Such studies may lead to the identification of clinically applicable biomarkers to improve diagnosis, patient treatment and survival.

Materials and Methods
Global transcriptome and miRNA expression data were retrieved from 45 samples (36 CCA and 9 histologically normal biliary tissues) profiled by TCGA (http://cancergenome.nih.gov/). Raw data (level 3, read counts) were quantified by counting transcriptome sequencing (RNA-Seq) and miRNA sequencing (miRNA-Seq) reads. Quantitative data were used to calculate differential mRNA and miRNA expression levels in CCA relative to normal biliary tissues. R/Bioconductor EBSeq [26] and edgeR package [27] were used in their default settings, for mRNA (gene) and miRNA expression analyses and Multidimensional scaling (MDS) analysis, respectively. A list of differentially expressed (DE) genes and miRNAs was generated and fold change (FC) expression values were determined. Genes and miRNAs were identified as over-or under-expressed with FDR ≤ 0.05. Differentially expressed genes with FC ≥ 3 and miRNAs with FC ≥ 1, based on their quartiles (Q3), were used in subsequent analyzes. The lists of differentially expressed genes and miRNAs in CCA compared to normal tissues were overlapped and miRNA target genes were identified. Identification of miRNA target genes was performed using microRNA Data Integration Portal (mirDIP) (http://ophid.utoronto.ca/mirDIP) [28]. External data validation was performed to compare differentially expressed miRNAs identified through analysis of TCGA data with the current literature (Supplementary Table 1). Furthermore, we searched DisGeNET, a database containing gene-disease associations (http://www.disgenet. org/web/DisGeNET/menu/home) [29] with the goal of enriching deregulated genes in CCA and that were also found deregulated in the TCGA dataset.
Next, in order to achieve a systematic understanding of deregulated networks associated with disease biology, miRNA target genes and their respective proteins were used to construct protein-protein interaction (PPI) networks. PPI networks were identified using STRING, a search tool for the retrieval of interacting Genes/Proteins (http://string-db. org/) [30]. The DGIdb database (http://dgidb.genome.wustl.edu/) [31] and STITCH (http://stitch.embl.de/) [32] were used to identify potential interactions between genes and proteins. Visualization and analysis of network properties of graphs such as Degree [33] were performed in Cytoscape [34,35].

Results
Results presented herein are in whole based on our analysis and interpretation of data generated by the TCGA research network (http:// cancergenome.nih.gov/). We identified 64 miRNAs differentially expressed between CCA and corresponding normal tissues from the TCGA dataset (Supplementary Table 2). Of these, 48 were over-and 16 were under-expressed, as described in Table 1. In addition, we identified 432 differentially expressed genes between CCA and normal tissues with FDR<0.05. Of the 432 genes, 180 were over-and 252 were underexpressed. Supplementary Table 3 lists all genes identified including this 432 gene subset. The top 20 over-and under-expressed genes are highlighted in Table 2.
To visualize the separation between the different classes of samples (normal and tumor), a MDS graph was elaborated with the data set of miRNAs available for cholangiocarcinoma ( Figure 1). Among the 48 over-expressed miRNAs, miR-125a was identified to have the highest number of direct interactions (33 mRNA targets), followed by miR-27a (17 targets), miR-127 (15 targets) and miR-200a (14 targets) ( Figure  2). A subset of genes was found to interact with a large number of miRNAs, such as the under-expressed genes ACADSB, which interacts with 19    15 miRNAs. Predicted interactions between under-expressed miRNAs and over-expressed genes in CCA are shown in Figure 3. Among the under-expressed miRNAs, miR-122 and miR-139 were found to have the highest number of interactions (9 targets each), followed by miR-424 and miR-194 with 8 targets each. miR-122 was predicted to target the transcription factor FOXM1, the matrix metalloproteinase MMP7 and the TGFß regulator PMEPA1. Among these genes, FOXM1 and MMP7 were specifically deregulated in CCA cells. In addition, miR-148 and miR-194 were predicted to target NAD(P)H quinone dehydrogenase 1 (NQO1) gene. The over-expressed genes PMEPA1 and ITGA2, ITGA3 and IGSF3 were found to be targeted by at least five miRNAs. PPI networks including enrichment analysis of biological functions were constructed for over- (Figure 4) and under-expressed genes ( Figure 5). Furthermore, we identified potential drugs with validated drug-gene interaction data ( Figure 6).

Discussion and Conclusion
Since miRNAs are potent gene expression regulators targeting multiple genes in several pathways, precise identification of miRNA expression in cancer compared to paired normal tissues can shed light into the mechanisms that determine cell fate and function related to disease development and progression.
Herein, deregulated expression of miRNAs was found to modulate a wide range of genes in disease-associated pathways. Among these, miR-365, which plays a role in chromatin modification, was downregulated in CCA compared to normal biliary tissues (Supplementary Table 2) and has also been described in other tumor types [36]. miR-365a also regulates genes related to epithelial to mesenchymal transition, cell cycle progression, and PI3K/AKT and MAPK pathways. Notably, a chromatin modifier gene expression signature was reported in cholangiocarcinoma tumor subtypes harboring IDH mutations [25].
Previous studies in other cancer types [37-40] have also provided supporting evidence for the involvement of genes and miRNAs identified as deregulated in CCA. Interestingly, miR-24 up-regulation was also detected in Acute Myeloid Leukemia [41] and in hepatocellular carcinoma being correlated with larger tumor size, higher microvessel density and tumor dedifferentiation, which are important factors associated with cancer cell invasion [42].
Among known miRNA clusters, miR-183 (miR-183, miR-96 and miR-182) was over-expressed in CCA compared to normal biliary tissues. miR-183 cluster has been reported as highly expressed in cancer   [43]. miR-182 alone was found to be up-regulated in hepatocellular carcinoma [44,45] and suggested to act as a metastasis-inducing miRNA in the liver [46]. miR-182 may be a good candidate for functional validation studies as a therapeutic target in CCA.
miR-200 family was also over-expressed in CCA; this miRNA family consists of five members distributed in two chromosomes: cluster I on chromosome 1 (miR-200a, miR-200b, and miR-429) and cluster II on chromosome 12 (miR-200c and miR-141) [47]. Some members of this cluster, namely, miR-200a, miR-200b and miR-200c have increased expression in CCA [48]. miRNAs miR-21, miR-141 and miR-200b were reported as highly expressed in CCA and suggested to have a role in chemotherapy, since inhibition of miR-21 and miR-200b resulted in increased sensitivity of tumor cells to gemcitabine and inhibition of miR-141 led to decreased tumor cell growth [49]. miR-200 family has been also reported to be up-regulated in bladder carcinoma [50] and both miRNA clusters, miR-183 and miR-200 are positively regulated in ovarian cancer [51]. These miRNAs were suggested as potentially useful biomarkers for diagnosis and treatment of lung cancer [52].
The miR-148/152 cluster (miR-148a, miR-148b and miR-152) has been reported as commonly down-regulated in cancer [53], including CCA [54]. Down-regulated expression of miR-148a and miR-152 were reported in breast [55] and gastric cancer [56,57], and under expression of miR-148a was reported in pancreatic cancer [58,59] and colorectal carcinoma [60]. Of note, in vivo replacement of miR-148a in hepatocellular carcinoma cells was shown to suppress and prevent tumor growth and reduce liver fibrosis [61]. Interestingly, a predicted target of the under-expressed miR-148, NQO1 gene, was found to have a large number of interactions and to be up-regulated in CCA. NQO1 over-expression occurs in different tumors compared to their normal surrounding tissues [62] and was recently associated with poor prognosis in serous ovarian cancer [63]. Cholangiocarcinoma is highly resistant to chemotherapy and NQO1 protein, NAD(P)H-quinone oxidoreductase 1, has been shown to play a cooperative role in cancer chemoresistance [64]. Recent data by Zeekpudsa et al. [65] showed that NQO1 knock-down increased the response of CCA cells to the chemotherapeutic drugs 5-fluorouracil and gemcitabine. Altogether, these findings suggest a potential role for NQO1 in treatment of patients with CCA.
Under-expression of miR-101 has been identified herein as well as in previous studies in liver cancer [66,67] and CCA [68]. Down-  regulated miR-101 was associated with hepatitis B virus infection [69]. In CCA, miR-101 may modulate vascular endothelial growth factor (VEGF) and cyclooxygenase-2 (COX-2) expression levels leading to a decrease in tumor cell growth [68]. miR-101 is a putative regulator of genes linked to cell adhesion, including ITGA2, which has increased expression in CCA [70], as well as ITGA3 and CDH6.
According to our analyses, miR-125a was the over-expressed miRNA with the highest number of regulated target genes. Up-regulated miR-125a was also reported in prostate cancer [71]. However, other studies have reported under-expression of miR-125a and suggested a tumor suppressive role for this miRNA [72][73][74][75]. miR-125a can regulate the transcription of ALB, PROC and MST1, which are deregulated in CCA according to DisGeNET data. Interestingly, MST1 inactivation in mammalian liver has been demonstrated to increase hepatocyte proliferation and hepatomegaly followed by the development of hepatocellular carcinoma and CCA [76]. Among the under-expressed miRNAs, miR-122 targets the highest number of genes. miR-122 is down-regulated in hepatocellular carcinoma [77] and has a role in tumor development associated with the hepatitis C virus [78], since this miRNA is needed for viral replication [79]. miR-122 has been also shown to negatively regulate FOXM1, which has increased expression in CCA [70] and has been demonstrated to promote epithelial to mesenchymal transition and metastasis in hepatocellular carcinoma [80]. Inhibition of FOXM1 was suggested to have a role in cancer treatment [81] and miR-122 has been indicated as a potential target for replacement therapy in hepatocellular carcinoma [82].
Among the 180 genes found to be over-expressed in our analysis of this TCGA dataset, SPP1 (Osteopontin) was a hallmark. SPP1 overexpression has been reported in intrahepatic CCA [83]. SPP1 is involved in tumor cell adhesion, migration and invasion [84] and its serum levels were suggested to have prognostic value in hepatocellular carcinoma [85]. SPP1 over-expression has also been found in lung cancer [86], oral squamous cell carcinoma [87], colorectal [88], gastric and liver cancers [89]. SPP1 protein was associated with metastasis [90], demonstrated to interact with integrins [91] and matrix metalloproteinases (MMPs) in cancer [92].
MMPs are known extracellular matrix remodeling proteins with key roles in tumor cell migration, invasion, apoptosis, angiogenesis and immune cell migration [93]. MMPs have been suggested as prognostic biomarkers in CCA [94]. MMP11 is potentially regulated by the miRNAs let-7c and miR-148a, both with reduced expression levels as identified in our data analysis. MMP11 promotes cellular proliferation, angiogenesis and decrease cancer cell sensitivity to NK cells [95]. MMP11 is highly expressed in gastric tumors and siRNA inactivation of MMP11 led to decreased cancer growth [96]. Motrescu and Rio [97] showed that MMP11 is expressed in the tumor microenvironment in the adipose tissue around the tumor, acting in adipocyte dedifferentiation and leading to accumulation of non-malignant peritumoral fibroblastlike cells, which support tumor cell survival. In normal liver, MMP11 demonstrates constitutive expression [98] having IGF-BP1 (insulin-like growth factor binding protein 1) as a substrate and thus cooperating to increase the availability of IGF-1 [99]. Interestingly, MMP11 was the most significantly up-regulated gene in a CCA microarray study [100].
Our transcriptome analyses also identified MMP7 (Matrilysin) overexpression. One MMP7 substrate is E-cadherin, an adherens junction protein involved in cell to cell interaction [101]. MMP7 has roles in tumorigenesis by cooperating with cellular growth and proliferation, having anti-apoptotic effects and inducing epithelial to mesenchymal transition and cell migration (reviewed in Gialeli) [95]. Abnormal MMP7 protein levels were detected in CCA [99]. MMP7 is potentially regulated by miR-122, which was broadly under-expressed and targeting a large number of genes as identified herein. MMP7 has been suggested as an useful biomarker in lung cancer [101] and CCA [98]. MMPs may be targeted by synthetic inhibitors such as Marimastat, a broad spectrum MMP inhibitor [95] that was experimentally tested and shown useful for treatment of polycystic liver disease, which presents cholangiocyte overgrowth by MMP hyperactivity [102]. In addition, Batimastat and Prinomastat are synthetic MMP inhibitors and Neovastat is a natural MMP7 inhibitor [95]. Despite this and other published studies, there is a need to understand, at the protein level, the effects of miRNA-regulated and disease-associated pathways. Nevertheless, our study contributes for a better understanding of miRNA-modulated signaling pathways involved in CCA pathobiology. Future in vitro and in vivo functional validation studies targeting miRNAs and pathways are required for the development of novel treatment strategies for patients.