Integration analysis of miRNA–mRNA expression exploring their potential roles in intrahepatic cholangiocarcinoma

Intrahepatic cholangiocarcinoma (ICC) is the second common primary hepatic malignancy tumor. In this study, an integrative analysis of differentially expressed genes (DEGs) and miRNAs from the ICC onset and adjacent normal tissues were performed to explore the regulatory roles of miRNA–mRNA interaction. A total of 1018 DEGs and 39 miRNAs were likely involved in ICC pathogenesis, suggesting the changes in cell metabolism in ICC development. The built network indicated that 30 DEGs were regulated by 16 differentially expressed miRNA. The screened DEGs and miRNA together were probably considered the biomarkers of ICC, and their important roles in ICC pathogenesis remain to be elucidated. This study could provide a good basis to uncover the regulatory mechanism of miRNA and mRNAs in ICC pathogenesis.


Materials and methods
Patient tissues preparation and RNA extraction. The resected ICC tissues and adjacent normal tissues were collected from curative surgey for ICC patients hospitalized at the First Affiliated Hospital of Soochow University during the period of February 2018 to October 2020. The Samples of ICC tissues (ICC1, ICC2, and ICC3) and adjacent normal tissues (Nor1, Nor2, and Nor3) were from three different patients at the TNM II stage according to the condition of tumor development (Table S1). This study was approved by the Ethics Committee of the First Affiliated Hospital of Soochow University (No.2021-335) and the experiments with human samples confirmed to the guidelines set by the Declaration of Helsinki. Written informed consent was obtained from participants before their enrollment in the study. For RNA extraction, all samples of ICC and adjacent normal tissues were quickly frozen in liquid nitrogen after resecting. The total RNAs for whole transcriptome sequencing were extracted using Trizol reagent (Invitrogen). The total RNAs for small RNA sequencing were extracted using the mirVana miRNA Isolation Kit (Ambion). An Agilent 2100 Bioanalyzer (Agilent Technologies) was applied to the evaluation of RNA integrity, and the samples with an RNA integrity number (RIN) ≥ 7 were used for the subsequent analysis.
Transcriptome sequencing and analysis. The rRNA-depleted and strand-specific RNA library con-  22 . GO annotation and KEGG enrichment were performed by using R based on a hypergeometric distribution test 23-25 . miRNA-mRNA interaction network construction and GO/KEGG analysis. The normalized read counts of DEGs and differentially expressed miRNAs were extracted and used to construct miRNA-mRNA regulatory networks according to the negatively correlated miRNAs and mRNAs. GO annotation and KEGG enrichment analysis of differentially expressed miRNA with their target genes were performed using R based on the hypergeometric distribution test [23][24][25] .

Results
Analysis of the mRNA expression pattern in ICC tissues. To uncover mRNAs expression pattern in the ICC development and progression, three samples from the ICC tissues (ICC1, ICC2, and ICC3) was analyzed with RNA sequencing against the adjacent normal tissues (Nor1, Nor2, and Nor3). All the raw data were deposited into the public database with accession numbers PRJNA763019 and PRJNA763017. An average of 94,982,664 and 91,935,189 clean reads were obtained from the samples of ICCs and Nors, respectively. The mapped reads to the reference genome (GRCh38.p12) account for more than 98% of the total reads from both groups (Table S2). The read distribution indicated good sequencing depth and coverage (Fig. 1). The principal component analysis (PCA) suggested a distinct and diversified gene expression pattern in the ICC tissues in comparison with the adjacent normal tissues ( Fig. 2A). To analyze gene expression patterns, the DEGs from ICCs against Nors were visualized with the volcano plot. The criterion of fold change ≥ 2.0 and adjusted P value ≤ 0.05 filtered out 947 up-regulated genes and 1,298 down-regulated genes, whereas the up-regulated and down-regulated gene numbers were lessened to 374 and 644, respectively, by using a strict criterion of fold change ≥ 4.0 and adjusted P value < 0.01 ( Fig. 2B and Table S3). These DEGs were hierarchically clustered in an unsupervised model, and two samples were separated in the heatmap (Fig. 2C).

KEGG pathway and GO enrichment analyses of DEGs in ICC tissues.
To clarify the potential function of the DEGs in ICC tissue, the up-regulated and down-regulated DEGs were subject to GO analysis and KEGG enrichment analysis, respectively. The top 10 GO terms from categories of biological process, cellular component, and molecular function were shown in Fig. 3A. The Up-regulated DEGs were mainly associated with cell growth and proliferation in the biological process. many structural cellular components such as the nucleus, cytoplasm, and chromatin were significantly enriched. Among the top 10 molecular functions, protein binding, ATP binding, and microtubule binding were significantly (P ≤ 0.05) enriched. On the contrary, the down-regulated DEGs were mainly associated with cell responses and their communications in biological processes. The components of extracellular exosomes, endoplasmic reticulum, and peroxisomes were significantly enriched, and the molecular functions included activity modulation by small compound binding. www.nature.com/scientificreports/ analysis revealed 30 and 56 metabolic pathways for the enrichment of the up-regulated and down-regulated DEGs, respectively. The top 20 pathways were categorized into the cellular process, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. To note, many pathways associated with metabolisms of amino acids and xenobiotics were repressed, and the activated pathways were related to several human disease factors, including microRNAs (miRNA) in cancer (Fig. 3B).
Of the enriched gene set (BRCA1, CCNE2, CDCA5, CDKN2A, E2F1, E2F2, KIF23, SOX4, TGFB2, TRIM71; 10/202 of hsa06206), breast cancer type 1 susceptibility protein BRCA1 specifically mediates the formation of Lys-6-linked polyubiquitin chains and plays a central role in DNA repair by facilitating cellular responses to DNA damage, E3 ubiquitin-protein ligase TRIM71 can behave the post-transcriptional mRNA repression in a miRNA independent mechanism. The E3 ubiquitin-protein ligase activity is required for its tumor suppressor function, and all these enriched genes are engaged in miRNA interactions for the development of some cancers. Altogether, these results revealed that ICC development did change cell metabolism, and the DEG pattern was probably linked to miRNA regulation.

Analysis of the miRNA expression pattern in ICC tissues.
To identify the differentially expressed miRNAs in ICC tissues, small RNA libraries were constructed from the samples of ICC tissues. The clean reads were in a range of 13,741,533-15,498,480 in both ICCs and Nors, and more than 92.1% could be mapped to the reference genome of humans. To note, around 30-40% of reads were aligned to miRbase version 22.0 in both samples (Fig. 4A), and their average abundance was dozen times higher than that of the total reads for both samples (Fig. 4B). The PCA plots showed a distinct miRNA expression pattern of the ICCs to that of the Nors, which was in line with the observation of their mRNA expression patterns (Fig. 2B). A total of 521 aligned (or predicted) miRNAs with baseMean of ≥ 10 were filtrated at the criterion of fold change ≥ 2.0 and P value < 0.05 to screen the deferentially expressed miRNAs. As a result, there were 12 significantly up-regulated and 27 significantly down-regulated miRNAs (Fig. 4C, and Table S4). Within an unsupervised model, these differentially expressed miRNAs could be grouped according to the sample source (Fig. 4D), and they were also clearly separated in the heatmap (Fig. 4E).

Target prediction for the differentially expressed miRNAs and function analysis.
To understand the potential roles of the differentially expressed miRNAs in ICCs, their target mRNAs were predicted with miRanda algorithms (Table S5). Twenty-one down-regulated miRNAs were predicted to regulate 365 targets via 380 interactions, and eleven up-regulated miRNAs probably interplay with 309 predicted targets via 322 interactions (Fig. 5A,B). These targeted candidates were further subjected to GO terms and KEGG enrichment analysis. The GO enrichments indicated that the critical biological process, cellular component, and molecular function focused on the Wnt signaling pathway, calcium modulating pathway, and those related to cell proliferation, axons, synapses, and protein stabilization (Fig. 5C). The KEGG pathway analysis revealed that these targeted genes were categorized into environmental information processing, human disease, and organismal system (Fig. 5D). Many pathways are correlated to cancer development, hinting at the crucial roles in the ICC pathogenesis.

Discussion
In this study, the integrative analysis of mRNA and miRNA expression in ICC was performed with transcriptome sequencing. The potential functions of DEGs and differential miRNAs were comprehensively explored to understand the miRNA-mRNA regulatory mechanism on the ICC pathogenesis and progression. There are 1018 DEGs identified from the ICC tissues, suggesting the aberrant transcriptional expression and their crucial roles in the ICC etiology and pathogenesis. Functional and KEGG analysis showed that these DEGs s were related to functions of cell proliferation, protein metabolisms, and extracellular matrix structural constituent, as well as many human diseases-relevant signaling pathways. Some genes could be further screened as prognostic markers of the ICC. Moreover, the upregulated DEGs were enriched in pathways in cancer, human papillomavirus (HPV)    1,[26][27][28] . We guessed that factors like alcohol use, chemical contamination, and metabolic syndrome could be the shared risks for ICC and other liver diseases. Molecular pathway level analysis found that the dysregulated expression of DEGs in ICC was associated with the metabolism of aberrant tRNA, amino acid, and lipoprotein 29 . The increasing pieces of evidence also revealed that aberrantly expressed miRNAs were common in various types of human cancer and play important roles in tumorigenesis [30][31][32][33] . Thank to the high-throughput sequencing of small RNAs, 39 differentially expressed miRNAs were filtrated from the ICC tissues. The clustering analysis showed the distinct expression profiles of miRNAs between ICC tissues and adjacent normal tissues, suggesting a potential role of these miRNAs as biomarkers for ICC diagnosis and therapeutic targets. It was to note that cholangiocarcinoma was originated from the epithelium of the bile duct, and most cancerous tissue was consisted adenocarcinoma cells from bile duct [34][35][36] . Such an aberrant expression profile of miRNA-mRNA might be related to bile duct lesions. Given the analyzed ICC specimens were at TNM II stage, the clinical stage of perhaps also contributes to this observation. The miRNA expression of miR-433, miR-22, miR-21, miR-125b-5p, miR-551b-3p, and miR-182 have been reported to regulate the progress or pathogenesis of ICC via regulating target genes expression [37][38][39][40] . For example, miR-182 was maily expressed in ductular reaction cells and was associated with cholangiocyte damage and ductular reaction. Knockdown of miR-182 resulted in reduced bile acid accumulation and cholestasis, even reduced live injury and inflammation 41 . In this regard, It was hypothesized this proposed miRNA/mRNA network might be infection would reinforce in ICC progression 32 , whereas our results did not suggest HCV infection is involved in ICC progression. As none analyzed case suffered from HCV infection in this study, it remained to elucidate its contribution by sampling the HCV infected cases. The present study also predicted the miRNA-mRNA interaction in ICC pathogenesis via GO term and KEGG analysis. Of 39 miRNAs, 32 candidates were predicted to regulate 615 targets, which were involved in the Wnt signaling pathway, calcium modulating pathway, and those related to cell proliferation, axons, synapses, and protein stabilization. These results demonstrated that differentially expressed miRNAs may function as a crucial regulator in ICC pathogenesis by modulating the target genes. The regulation was also founded on the important amino acid and fatty acid metabolic pathways, which was in line with the previous observation 29,42 . But be aware that the analyzed samples were collected from surgical resection rather than laser capture microdissection, so the moderately precise separation of ICC specimens would cause condemination of epithelial cells from www.nature.com/scientificreports/ stromal cells, and increase the background noise to discriminate the clonal origins of ICC [43][44][45][46] . Moreover, ICC cells interact with a complex milieu of supporting stromal cells that form the tumor microenvironment. single cell atlases provide important insight into understanding the contribution of different cell clusters to the malignant process in a spatial view 46 . Such an analysis distinguish the clonal origins of the subtupes of the combined HCC-ICCs 47 . Therefore, a more comprehensive biological and clinical insight into ICC could be excepected by sampling operation with aser capture microdissection and futher integration of signle cell atlases. In conclusion, the present study suggested a set of mRNAs and miRNAs involved in ICC pathogenesis by integrative analysis of their expression profiles in the ICC tissues against adjacent normal tissues. It would be a good case to uncover the regulatory mechanism of miRNA and mRNAs in ICC pathogenesis, which could assist us with development of new therapies for ICC.

Data availability
Raw sequence reads for mRNA (Accession: PRJNA763017) and miRNA (Accession: PRJNA763019) were deposited in NCBI database. Other data are available from the corresponding author upon reasonable request. www.nature.com/scientificreports/ Reprints and permissions information is available at www.nature.com/reprints.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.