Circulating microRNAs as biomarkers in bile-derived exosomes of cholangiocarcinoma.

PURPOSE: In this pilot study, using next-generation sequencing and integrated messenger RNA (mRNA) sequencing, we investigated circulating microRNA (miRNA) expression profiling from bile-derived exosomes to identify dysregulated miRNA signatures and oncogenic pathways and determine their effects on targeted mRNAs in cholangiocarcinoma (CCA). Moreover, we explored the possibility that genetic analysis using bile-derived exosomes may replace gene analysis using tissue.
METHODS: Bile was collected from a patient with perihilar CCA before curative resection. As a control, bile was collected from a patient with a common bile duct stone. Exosomes were isolated from the bile, and we performed next-generation miRNA sequencing using isolated exosomes. To evaluate miRNA-mRNA interactions, mRNA sequencing was performed using bile fluid in both patients.
RESULTS: We identified 22 differentially expressed miRNAs. More than 65% of the predicted mRNA targets of those miRNAs were actually differentially expressed between control and CCA bile samples. In functional pathway analysis, targets of 22 miRNAs were primarily enriched in mitogen-activated protein kinase, platelet derived growth factor, vascular endothelial growth factor, epidermal growth factor receptor, and p53 signaling. In particular, in the functional assessment of miRNA-mRNA interactions, RAS pathways, including downstream pathways (PI3K-AKT-mTOR and RAS-RAF-MEK-ERK), were determined to be enriched.
CONCLUSION: Circulating miRNAs in bile-derived exosomes provide new information for the development of miRNA analysis in CCA. These miRNAs may represent the oncogenic characteristics of CCA tissue, enabling them to be used instead of tissue samples for the diagnosis of CCA. Further research investigating circulating miRNAs in bile exosomes may lead to more rational, targeted approaches to treatment.

INTRODUCTION

Cholangiocarcinoma (CCA) is a malignant disease originating from cholangiocytes. A tissue biopsy is necessary for the diagnosis of CCA, and this procedure is difficult to perform routinely due to its invasiveness. The diagnostic value of serum biomarkers, such as CEA and CA 19-9, is very limited. Moreover, many patients have no specific symptoms, early diagnosis remains difficult, and overall survival for patients with CCA is dismal [1].

Exosomes are small vesicles (30–100 nm in diameter) secreted by most cell types into body fluids, including bile fluid [2,3]. Exosomes contain mainly fusion proteins and transporters, phospholipases, and other lipid-related proteins [4,5]. Tumor-derived exosomes carry nucleic acids and oncogenic proteins that contribute to the microenvironment necessary for tumor progression [6,7]. Since exosomes carry proteins, messenger RNA (mRNA) and microRNA (miRNA) can transfer signals from one cell to another [8], they may be ideal biofluids for cancer research. Since exosomes function as potent signaling molecules between the tumor and the surrounding tumor microenvironment, exosomes may represent not only tumor characteristics but also the tumor microenvironment. Therefore, studies highlighting the role of exosomes as facilitators of cancer progression have recently been increasing [9]. Numerous studies have also clearly demonstrated that exosome-containing miRNAs contribute to many events related to cancer progression, including cancer initiation, proliferation, invasion, metastasis, and induction of angiogenesis [10]. In CCA, due to the higher stability of noncoding RNA in exosomes than in body fluids, genetic analysis via exosomes would be more desirable than analysis of bile fluid itself.

However, genetic analysis using bile has not been commonly performed at present. Bile fluid is in direct contact with the tumor cells in CCA. Therefore, tumor-derived materials may be abundant in the bile fluid. Although a few studies have investigated miRNA profiling from bile-derived exosomes in CCA [11], integrative miRNA-mRNA analysis using next-generation sequencing (NGS) is rare. In this pilot study, we investigated circulating miRNA expression profiling from bilederived exosomes using NGS and integrated mRNA sequencing to show dysregulated miRNA signatures and oncogenic pathways and their effects on targeted mRNAs. Moreover, we attempted to explore the possibility that genetic analysis using bile-derived exosomes may replace genetic analysis using tissue.

METHODS

Patients

A total of 2 patients were enrolled at Keimyung University Dongsan Hospital between August 2018 and March 2020. Bile was collected from a 61-year-old male with perihilar CCA with a type IIIa Bismuth classification before curative resection. As a control, bile was collected from a patient with a common bile duct stone, this patient was an 87-year-old female. This prospective study was approved by the Institutional Review Board of Keimyung University Dongsan Medical Center (No. DSMC 2018-06-054) and included the secondary use of humanderived materials. Written informed consent was obtained from all patients.

Bile collection and preparation

Bile samples from both groups were collected at the time of percutaneous transhepatic biliary drainage (PTBD) for relieving cholangitis or obstructive jaundice (for CCA) and for relieving cholangitis by removing common bile duct stones (for the control). In both patients, bile fluid was collected after resolving cholangitis at least 3 days after the PTBD procedure. Collected bile fluid was centrifuged immediately at 16,000 ×g for 10 minutes at 4℃ to obtain supernatant (bile); each pellet was resuspended in chilled phosphate-buffered saline (2-fold) and then centrifuged at 16,000 ×g for 5 minutes at 4℃ to obtain a bile pellet, as described in the protocol by Abi Zabron Imperial College, London, England. All experiments were performed with Avanti J-25I (Beckman, Pasadena, CA, USA), and bile and pellets were stored at −80℃. After separation, we used these samples in all experiments.

Chemicals and reagents

An Exo2d-EV isolation kit was purchased from EXOSOMEplus (Seoul, Korea). QIAzol Lysis Reagent and a miRNeasy Serum/Plasma Kit were purchased from QIAGEN (Hilden, Germany). The TaqMan miRNA Reverse Transcription Kit and TaqMan 2X Universal PCR Master MiX were purchased from Life Technologies (Carlsbad, CA, USA). SMARTer smRNA-Seq kit purchased from Illumina Inc. (San Diego, CA, USA). Ethanol and all other chemicals were obtained from Sigma (St. Louis, MO, USA).

Exosome isolation

Exosomes were isolated using an Exo2d-EV isolation kit according to the manufacturer's instructions. Briefly, 1 mL of bile was centrifuged at 3,000 ×g for 15 minutes to remove cell debris and larger vesicles. The supernatant was mixed with Exo2d and incubated for 15 minutes at 37℃. After centrifugation at 3,000 ×g for 30 minutes, the exosome pellet was resuspended in 100 μL of phosphate-buffered saline for further analysis.

miRNA isolation and quantitative real-time polymerase chain reaction

miRNA purification was performed using the miRNeasy Serum/Plasma Kit according to the manufacturer's instructions. Bile exomes were lysed in QIAzol Lysis Reagent. After the addition of chloroform, the lysate was separated into aqueous and organic phases by centrifugation at 12,000 ×g for 15 minutes. RNA was partitioned to the upper aqueous phase, while DNA was partitioned to the interphase and proteins were partitioned to the lower organic phase or the interphase. The upper aqueous phase was extracted, and ethanol was added to provide appropriate binding conditions for all RNA molecules from approximately 18 nucleotides upwards. The sample was subsequently applied to the RNeasy MinElute spin column, where the total RNA bound to the membrane, and phenol and other contaminants were efficiently washed away. High-quality RNA was subsequently eluted in a small volume of RNase-free water.

SMARTer smRNA library

The RNA isolated from each sample was used to construct sequencing libraries with the SMARTer smRNA-Seq Kit for Illumina following the manufacturer's protocol. In brief, input RNA was first polyadenylated to provide a priming sequence for an oligo (dT) primer. Complementary DNA (cDNA) synthesis was primed by the 3′ smRNA dT primer, which incorporated an adapter sequence at the 5′ end of each RNA template. This approach added nontemplated nucleotides that were bound by SMARTer smRNA Oligo-enhanced with locked nucleic acid technology for greater sensitivity. In the template-switching step, PrimeScript RT (TAKARA Bio Inc., Shiga, Japan) uses the SMARTer smRNA Oligo as a template for the addition of a second adapter sequence to the 3′ end of each first-strand cDNA molecule. In the next step, full-length Illumina adapters (including index sequences for sample multiplexing) were added during polymerase chain reaction (PCR) amplification. The forward PCR primer bound to the sequence added by the SMARTer smRNA Oligo, while the reverse PCR primer bound to the sequence added by the 3′ smRNA dT primer. The resulting library cDNA molecules included sequences required for clustering on an Illumina flow cell. The libraries were gel-purified by blue pippin and validated by checking their size, purity, and concentration on an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were pooled in equimolar amounts and sequenced on an Illumina HiSeq 2500 instrument to generate 51-base reads. Image decomposition and quality value calculations were performed using the modules of the Illumina pipeline.

miRNA sequencing

The raw reads of small RNAs were preprocessed to eliminate adapter sequences. Adapters in the raw reads were trimmed using the Cutadapt program (https://cutadapt.readthedocs.io). If a sequence was matched to more than the first 5 bp of the 3′ adapter sequence, it was regarded as an adapter sequence and was trimmed from the read. Trimmed reads, whose length was longer than 18 bp, were selected for mapping reliability. Next, the remaining reads were classified into nonadapter reads whose adapter sequences were not sequenced. In this analysis, trimmed reads and nonadapter reads were combined and regarded as processed reads for downstream analysis. To minimize sequence redundancy for computational efficiency, processed reads were clustered by a sequence. A unique cluster consists of reads whose sequences and length are the same. To eliminate rRNA, reads that aligned to the 45S pre-ribosomal RNA (rRNA) and mitochondrial rRNA of Homo sapiens were excluded. Sequence alignment and detection of known and novel miRNAs were performed using the miRDeep2 software algorithm (https://www.mdc-berlin.de/content/mirdeep2-documentation). rRNA-filtered reads were aligned to the mature and precursor miRNAs of Homo sapiens obtained from miRBase v22.1 (https://www.mirbase.org) using the miRDeep2 quantifier module. The miRDeep2 algorithm is based on the miRNA biogenesis model; it aligns reads to potential hairpin structures to determine if their mapping context is consistent with Dicer processing and assigns scores representing the probability that hairpins are true miRNA precursors. Novel miRNAs were predicted from the mature, star and loop sequences according to the RNA fold algorithm using miRDeep2. The RNA fold function used the nearest-neighbor thermodynamic model to predict the minimum free-energy secondary structure of an RNA sequence. Next, uniquely clustered reads were sequentially aligned to the reference genome, miRBase v22.1 and noncoding RNA database RNA central release 14.0 to identify known miRNAs and other types of RNA for classification.

mRNA sequencing library preparation and sequencing

We preprocessed the raw reads from the sequencer to remove low-quality and adapter sequences before analysis and aligned the processed reads to Homo sapiens (GRCh38) using HISAT v2.1.0 (https://ccb.jhu.edu/software/hisat2/index.shtml). HISAT utilizes 2 types of indexes for alignment (a global, whole-genome index, and tens of thousands of small local indexes). These 2 types of indexes are constructed using the same Burrows-Wheeler transform graph FM index as Bowtie2 (http://bowtie-bio.sourceforge.net/index.shtml). Because of its use of these efficient data structures and algorithms, HISAT generates spliced alignments several times faster than Bowtie and BWA (http://bio-bwa.sourceforge.net), which are widely used. The reference genome sequence of Homo sapiens (GRCh38) and annotation data were downloaded from the National Center for Biotechnology Information database. Next, the transcript assembly of known transcripts was processed by StringTie v2.1.3b (http://ccb.jhu.edu/software/stringtie/). Based on the results, the expression abundance of transcripts and genes was calculated as the read count or FPKM value (fragments per kilobase of exon per million fragments mapped) per sample.

RESULTS

Identification of differentially expressed miRNAs (DEmiRNAs) and genes (DEGs)

Originally, 2,656 mature miRNAs were identified in 2 samples. After trimmed mean of M-value normalization using the calNormFactors method in edgeR in the R library (https://rdrr.io/bioc/edgeR/man/calcNormFactors.html), miRNA counts per million reads with 0 in both samples were removed. For the remaining 537 miRNAs, further analysis was performed. Statistical analysis for identifying DEmiRNAs between 2 samples was performed using fold change, exactTest using edgeR (https://bioconductor.org/packages/edgeR) per comparison pair. The significant DEmiRNAs were selected on conditions of |fold change| ≥ 2. In this process, 52 mature miRNAs with |fold change| ≥ 2 and P < 0.05 (Fig. 1) were identified, and finally, a total of 22 DEmiRNAs with false discovery rate (FDR) < 0.05 were obtained (Table 1).

Fig. 1

Heat map of the 1-way hierarchical clustering using Z-score for normalized value (log2 based) based on 52 mature microRNAs satisfying |fold change| ≥ 2 and raw. P < 0.05.

Table 1

List of differentially expressed miRNAs

miRNA, microRNA.

For analysis of miRNA-mRNA interactions, DEGs between tumor and control tissues were analyzed using the edgeR package in the R library. After preprocessing, the significant DEGs were selected on conditions of |fold change| ≥ 2 with FDR < 0.05, and a total of 1,286 significant DEGs were identified.

miRNA-mRNA regulatory network

The miRNA-mRNA regulatory network was analyzed with Mienturnet (http://userver.bio.roma1.it) based on miRTarBase. Seventeen of 22 DEmiRNAs had predicted target genes (FDR < 0.05), and a list of them is listed in Supplementary Table 1. More than 65% of the predicted mRNA targets of 22 DEmiRNAs were actually differentially expressed between normal and cancer bile samples (Supplementary Table 1, Fig. 2).

Fig. 2

MicroRNA (miRNA)-messenger RNA (mRNA) regulatory network. Among 22 differentially expressed miRNAs (DEmiRNAs), 14 netDEmiRNA had a network of miRNAs–mRNAs that has a threshold for the minimum number of miRNA-target interactions of 4 and false discovery rate < 0.01. Bold blue boxes indicate differentially expressed mRNAs between normal and cancer bile samples. CCNT2, cyclin T2; ALDH5A1, aldehyde dehydrogenase 5 family member A1; DICER1, dicer 1, ribonuclease III; CNOT6, CCR4-NOT transcription complex subunit 6; ITGB1, integrin subunit beta 1; AKT3, AKT serine/threonine kinase 3; CTNNBIP1, catenin beta interacting protein 1; DNMT1, DNA methyltransferase 1; ZEB2, zinc finger E-box binding homeobox 2; WDR37, WD repeat domain 37; MYB, MYB proto-oncogene, transcription factor; DCBLD2, discoidin, CUB and LCCL domain containing 2; RPS6KA3, ribosomal protein S6 kinase A3; EGR1, early growth response protein 1; NLK, nemo like kinase; PLAG1, PLAG1 zinc finger; BCL2L11, BCL2 like 11; CASP8, caspase 8; BMPR2, bone morphogenetic protein receptor type 2; VHL, Von Hippel-Lindau syndrome; IRS1, insulin receptor substrate 1; CUL5, cullin 5; FOS, Fos protooncogene; PIK3CG, phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit gamma; NCOA3, nuclear receptor coactivator 3; SIRT1, sirtuin 1; PTEN, Phosphatase and tensin homolog; BCL2, B-cell lymphoma 2; VEGFA, vascular endothelial growth factor A; CDK6, cyclin dependent kinase 6; MYC, MYC proto-oncogene; MMP2, matrix metallopeptidase 2; DNMT3B, DNA methyltransferase 3 beta; LOX, lysyl oxidase; ADAM12, ADAM metallopeptidase domain 12; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; SERPINH1, serpin family H member 1; COL4A2, collagen type IV alpha 2 chain; FBN1, fibrillin 1; LAMC2, laminin subunit gamma 2; COL5A2, collagen type V alpha 2 chain; TDG, thymine-DNA glycosylase; SPARC, secreted protein acidic and rich in cysteine; TET2, Tet methylcytosine dioxygenase 2; IRS1, insulin receptor substrate 1; BMPR2, bone morphogenetic protein receptor type 2; PIK3R1, phosphoinositide-3-kinase regulatory subunit 1; COL4A1, collagen type IV alpha 1 chain; AKT2, AKT serine/threonine kinase 2; LRP6, LDL receptor related protein 6.

Pathway and functional analysis

Functional annotation and miRNA pathway clustering of 22 DEmiRNAs was performed using miRNA Enrichment analysis and Annotation Tool (miEAA; https://www.ccb.uni-saarland.de/mieaa2), which is based on Gene Trail [12]. Functional assessment showed that their targets were mainly enriched in RAS, mitogen-activated protein kinase (MAPK), mechanistic target of rapamycin (mTOR), platelet-derived growth factor (PDGF), vascular endothelial growth factor (VEGF), epidermal growth factor receptor (EGFR), and p53 signaling and were dysregulated. In addition, inflammation and cell cycle-related pathways were also dysregulated. The top 30 significantly enriched pathways with FDR < 0.05 are listed in Table 2.

Table 2

Top 30 significantly enriched pathways of 22 differentially expressed miRNAs

miRNA, microRNA; Jak-STAT, Janus kinase and signal transducer and activator of transcription; TGF, transforming growth factor; MAPK, mitogen-activated protein kinase; TCR, T-cell receptor; IL, interleukin; DDR, DNA damage-regulated; VEGF, vascular endothelial growth factor; EGF, epidermal growth factor; EGFR, epidermal growth factor receptor; ECM, extracellular matrix; ErbB, erythroblastic leukemia viral oncogene homolog; Wnt, wingless and int-1.

To assess the functional pathways of the DEGs, pathway enrichment analyses were performed using the online database Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8; https://david.ncifcrf.gov) [13].

In the functional assessment of miRNA-mRNA interactions, RAS pathways, including downstream pathways (PI3K-AKT-mTOR and RAS-RAF-MEK-ERK), were enriched in both miRNA and mRNA functional pathway analysis (Fig. 3) [14].

Fig. 3

RAS pathway in functional assessment of microRNA (miRNA)-messenger RNA (mRNA) interaction. Highlighted red-colored and yellow-colored boxes indicate differentially expressed miRNA and mRNA signaling in the present study, respectively. This image was obtained by Kyoto Encyclopedia of Genes and Genomes with copyright permission (Kanehisa [14]). RAS, rat sarcoma virus; IKK, inhibitor of nuclear factor-κB (IκB) kinase; NF-κB, nuclear factor kappa-light-chain-enhancer of activated B cells; PI3K, phosphoinositide 3-kinase; PIP3, phosphoinositide 3-kinase; AKT, protein kinase B; BAD, BCL2 associated agonist of cell death; BCL-X, B-cell lymphoma-extra large; AFX, FOXO4; FasL, Fas ligand; AF6, ALL1-fused gene from chromosome 6 protein; Raf-1, rapidly accelerated fibrosarcoma-1; MAPK, mitogen-activated protein kinase; MEK, mitogen-activated protein kinase kinase; ERK, extracellular-signal-regulated kinase; KSR, kinase suppressor of Ras; PLA, poly-lactic acid; ELK, ETS transcription factor; ETS, electron transport chain; Repac, another Epac subfamily member called.

DISCUSSION

miRNAs posttranscriptionally regulate gene expression via mRNA decay and translational repression [15]. Several miRNAs are involved in carcinogenesis by acting either as oncogenes or tumor suppressor genes [16]. Dysregulation of miRNAs has been observed in CCA tissues, suggesting the potential roles of miRNAs in CCA progression [17,18]. Considering that early diagnosis by tissue biopsy is limited in CCA, liquid biopsy for the detection of cancer-derived circulating miRNA in blood or body fluid may represent an ideal biomarker method for diagnostic screening [19,20]. However, the diagnostic accuracy of circulating miRNA in blood varies in many studies [21,22]. Moreover, discrepancies among the results obtained with blood and tissue samples are often observed [21]. Therefore, investigation of circulating miRNA in blood may not represent genetic characteristics of tumor tissue, and the diagnostic value of blood circulating miRNA for CCA remains limited. In the case of CCA, bile fluid directly contacts tumor cells, and tumorderived materials may be abundant in bile, and considering that miRNAs enclosed in extracellular vesicles can increase their half-life in circulation [23], miRNA profiling in bile-derived exosomes may be an ideal biofluid for replacing tissue biopsy.

The present study revealed that several miRNAs up- or downregulated in bile exosomes, such as miR-191, miR-200b, miR-29a, miR-29b, and miR-21, were potential biomarkers of CCA tissue in previous studies [24]. miR-21 and miR-200b were overexpressed in CCA tissue and cell lines, and their high expression is associated with advanced CCA [22,25]. Similarly, miR-192 expression is associated with lymph node metastasis and a poor prognosis in CCA [22]. In the present study, more than 65% of predicted mRNA targets of DEmiRNAs were actually differentially expressed between normal and cancer bile samples. In addition, several oncogenic pathways of DEmiRNAs, including the RAS, VEGF, EGFR, p53, and transforming growth factor-beta signaling pathways, which are commonly seen in CCA tissues, were also enriched in bile exosomes. In particular, RAS pathways, including downstream pathways (PI3K-AKT-mTOR and RAS-RAF-MEK-ERK), were notably enriched. RAS proteins, and their downstream enzymes are activated in CCA tissues, and oncogenic mutations in RAS genes, including KRAS and NRAS, occur. The presence of KRAS mutations is crucial and is related to tumor progression and poor survival [26,27]. This finding suggests that bile-derived exosomes represent the oncogenic characteristics of CCA tissues and can replace tissue biopsy, which is more invasive.

We analyzed several novel potential biomarkers of bile exosomes. miR-10a-5p, which is overexpressed in pancreatic cancer and colorectal cancer, was also overexpressed in the present study. It is known that inhibition of miR-10a can suppress CCA through downregulation of the PI3K-AKT pathway [28]. Since the PI3K-AKT pathway was notably enriched in bile exosomes in the present study, miR-10a-5p in bile exosomes can be a potential diagnostic biomarker, as well as an individualized treatment target. miR-181a-5p, which is associated with cancer progression in cervical cancer and renal cell carcinoma [29,30], could be another potential bile-derived exosome biomarker in CCA. The combination of different markers, such as circulating miRNA in bile exosomes and blood CA 19-9, may improve the early diagnosis of CCA.

However, this study was limited, as only 2 cases were investigated. Nevertheless, considering that there are still no studies of miRNA-mRNA interactions by NGS in bile exosomes, our pilot study has value to further efforts to explore potential biomarkers. Further studies are warranted to validate our findings. In addition, we investigated molecular subtypes which have subtype-specific clinical features based on mRNA expression in previous study [26]. Although we could not show the relevance between clinical and genetic features in the present study because of small samples, we hope relationship between clinical information and miRNA expression in bile exosome can be evaluated by integrative clinical-molecular analysis with large samples in the future.

In conclusion, the results of this study indicate that circulating miRNAs in bile-derived exosomes provide new knowledge for the developing field of miRNA analysis in CCA. These exosomes may represent the oncogenic characteristics of CCA tissue and may be used instead of tissue biopsy for the diagnosis of CCA. Further research investigating circulating miRNAs in bile exosomes may lead to more rational, targeted approaches to treatment.