Potential prognostic value of hsa-miR-18a in hepatocellular carcinoma and the underlying key down-stream gene CHRM2: a comprehensive bioinformatic analysis*

Bioinformatic analyses were conducted to predict target genes and carry out enrichment analysis. Validated downstream genes of hsa-miR-18a were obtained from PubMed database. Differential expression analysis was conducted within the “edgeR” R package based on the TCGA datasets. Survival analysis was performed by Kaplan-Meier survival analysis. All the visualizations were implemented by R.


Results
Bioinformatic analysis obtained a total of 90 target genes of hsa-miR-18a and revealed that target genes were involved in pathways essential for cancer onset and development such as cell cycle and PI3K/AKT signaling pathway. A review of literatures found target genes of miR-18a indeed participating in the biological processes of HCC. CHRM2 was identi ed as a special gene after the intersection analysis of TCGA differentially expressed genes (DEGs) and predicted target genes. Survival analysis validated that hsa-miR-18a and CHRM2 signi cantly affected the prognosis of HCC patients.

Conclusion
There is a strong association between hsa-miR-18a with tumors including HCC via participating essential tumor-promoting pathways including cell cycle, PI3K/AKT signaling pathway, etc. Furthermore, high miR-18a expression and low CHRM2 expression could lead to a poor prognosis in HCC. In conclusion, miR-18a could serve as an expectational prognostic biomarker in HCC.

Background
As one of the most prevalent cancer in the alimentary system, hepatocellular carcinoma (HCC) ranks sixth in terms of morbidity and fourth in terms of mortality worldwide. Due to the remained di culties in early diagnosis and the high possibility of tumor relapse and metastasis after resection, HCC has been posing a huge threat upon human health and survival over the years [1]. Hence, profound study into the mechanism of tumorigenesis and progression, discovering biomarkers with high sensitivity and speci city aiming HCC for the early diagnosis and prognosis assessment is essential.
MicroRNAs (miRNAs) are short, single-stranded non-coding RNAs of around 22 nucleotides. Expression levels of target genes are down-regulated post-transcriptionally by miRNAs through binding to the 3'untranslated region (3'UTR) of their mRNAs, by which mean miRNAs participate in the transcriptional regulation of various cellular events including cell cycle, proliferation and protein phosphorylation and play key roles in tumorigenesis and malignancy development [2,3].
It has been revealed that aberrant expression levels of miRNAs in HCC cases are bound up with their biological behaviors and clinicopathological characteristics, making miRNAs undergoing extensive investigations [4][5][6]. Furthermore, precisely because of the identity of miRNAs mentioned above, their potentials of serving as in-situ or serum prognostic biomarkers has been acquiring increasing concerns from the oncology community [7][8][9]. However, relatively bulk bioinformatic minings and experimental validations are necessary before the prognostic value of select miRNA or miRNA panel is nally authenticated in clinical practice.
Meanwhile, hsa-miR-18a has been implicated in important biological processes of several cancer types including HCC, indicating the essentiality of comprehensively identifying downstream targets together with their functions as well as their prognostic in uences upon HCC [10].
In this article, the target genes of miR-18a were predicted and analyzed for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) to explore the functions and pathways possibly involved in HCC, and the results were demonstrated through literature researches. Next, differentially expressed genes were extracted from TCGA datasets and a test was made to identify whether an intersection with predicted genes was authentically existing. Finally, survival analysis was carried out to illustrate the diagnostic and prognostic monitoring value of miR-18a, and the underlying molecular mechanism was preliminarily explored.

Materials And Methods
Mature sequence of miR-18a in several species Sequences of miR-18a in different species were obtained from miRBase (http://www.mirbase.org) [11].
Using custom R scripts for mature sequence identi cation.
Expression level of target genes in normal tissue TiGER database (http://bioinfo.wilmer.jhu.edu/tiger/)[16] was used to obtain the expression level of target genes in various normal tissues. Only a gross picture of the expression was observed because several target genes could not be found in TiGER database.
GO and KEGG pathway enrichment analysis DAVID online tool (https://david.ncifcrf.gov/tools.jsp) [17] was used to conduct GO enrichment to identify the biological processes, cellular components and molecular functions of target genes as well as the KEGG pathway information [18]. The KEGG database (http://www.genome.jp/kegg/) was used for generating a certain pathway map. p < 0.05 was considered signi cant.

Protein-protein interaction (PPI) network construction
Interactions between target genes of hsa-miR-18a were obtained from the online database resource Search Tool for the Retrieval of Interacting Genes (STRING) (https://www.string-db.org) [19].

Data selection and statistical analysis
A total of 424 cases were extracted from TCGA database (https://www.cancer.gov/) with 374 HCC tissues and 50 normal tissues. The genes expressed differentially between tumor and normal tissues were analyzed and identi ed by "edgeR" R package. Survival curves were depicted with the Kaplan-Meier Plotter online tool (http://kmplot.com/analysis/) [20,21]. p < 0.05 was considered signi cant.

Data visualization
Visualization was conducted in R. The gene expression matrix, enrichment plots and differential expression boxplot were depicted by "ggplot2" R package [22]. The venn diagram was drawn with "VennDiagram" R package [23].

Literature research on hsa-miR-18a
A research based on PubMed database was performed to obtain preliminary details of the effect of miR-18a in various tumors ( Table 1). The results suggested a possible role for miR-18a as an onco-miRNA in cancers such as cervical cancer and gastric cancer because of the positive correlation between its increased expression and the tumorigenesis and progress. A study focused on nasopharyngeal cancer evidenced an increased expression of miR-18a in nasopharyngeal tissues, which can promote tumor proliferation and metastasis through down-regulating SMG1 expression [24,25]. However, another research demonstrated that decreased miR-18a level led to an over-expression of SREBP1, promoting the invasion of breast cancer through enhancing the epithelial-to-mesenchymal transition (EMT) [26]. Taken together, the results in Table 1 show that miR-18a has different characters in various tumors [27]. Mature sequence of miR-18a in various species Sequences of miR-18a in various species were obtained from miRbase database. As presented in Table 2, miR-18a is highly conserved during evolution and the mature sequence is UGUUCUAAGGUGCAUCUAGUGCAGAUAGUGAAGUAGA. The result suggests an important role of miR-18a in organisms. Target genes prediction and their expression map in normal human tissues Online databases TargetScan, miRDIP, miRWalk and miRDB were utilized and the numbers of target genes predicted were 5543, 982, 3078 and 944, respectively ( Fig. 1a.). The original prediction data are provided in Table S1. Only target genes which were presented by all four databases simultaneously were received as the gene set for further research. To elevate the comprehensiveness and feasibility for further work, ve genes which were proved to be the targets of miR-18a in miRecords database were added into the gene set. We obtained 90 target genes for further study ultimately as shown in Table 3.
Next, to further re ne our exploration into the function of miR-18a, we depicted a rough expression map of miR-18a in normal tissues base on TiGER database (Fig. 1b). After clustering and statistic analysis, no tissue speci city was observed within target genes of miR-18a. Genes included in TiGER database were provided in Table S2. The result reinforces the importance of unveiling the association between miR-18a and neoplastic diseases including HCC.
* Target genes obtained from miRecords database.

GO, KEGG enrichment analysis and PPI of miR-18a target genes
In order to elucidate the mechanisms of miR-18a in cancer, GO and KEGG enrichment analyses and PPI network construction were carried out further. 85 target genes obtained GO annotations, 37 target genes obtained KEGG annotations and 44 target genes participated in PPI network.
In terms of GO enrichment, nucleus (GO:0005634) and nucleoplasm (GO:0005654) were the most represented GO terms in the cellular component subontology (P < 0.05). In biological molecular function, DNA binding (GO:0003667), protein complex binding (GO:0032403) and protein serine/threonine kinase activity (GO:0004647) were relatively abundant terms (P < 0.05). The top enriched biological process terms were transcription, DNA-templated (GO:0006351), cell cycle (GO:0007409) and protein phosphorylation (GO:0006468) (P < 0.05). Table 4 presented full details and Fig. 2a visualized the number and p-value of GO enrichment.
KEGG pathway analysis revealed 12 signi cantly enriched pathways as shown in Table 5, especially in tumorigenesis-related pathways and metabolic pathways including cell cycle (hsa04110), PI3K pathway (hsa04151), Wnt signaling pathway(hsa04310), arginine biosynthesis (hsa00220) and angiogenesis (hsa04370) (P < 0.05). The number and p-value of all KEGG terms were shown in Fig. 2b. Notably, we observed in the pathway map provided by KEGG database that CCND1 and FGF played an essential role in cancer (Fig. 2c).
A total of 90 nodes and 51 edges were revealed in the PPI network with a PPI enrichment p-value < 0.05 (Fig. 2d). ATM, CCND1, CCND2, NR3C1 were identi ed as hub genes which might play key roles in the network and the function of miR-18a. Crucial genes unveiled in KEGG and PPI analysis made it reasonable to believe that these genes accompanied with miR-18a could be fruitful in the association with cancer.  information on the target genes (n = 14) and pathways (n = 3) modulated by miR-18a (Table 6). Entries retrieved for our literature research were connected by "OR" and were provided in Table S3. 11 out of 14 target genes validated by published studies focused on the function of tumorigeneses and progress while all three pathways identi ed participated in tumor-promoting functions including promoting EMT, metastasis and inhibiting DNA damage repair and apoptosis. Of note, a recent research conducted by Zheng et. al demonstrated that miR-18a was involved in tumorigenesis and progression of HCC via directly targeting NEDD9, which suggested miR-18a of holding latent value in its association with HCC.  [49] Potential of miR-18a and target gene CHRM2 as prognostic predictors of HCC patients To explore the association between miR-18a with HCC, differentially expressed genes (DEGs) between HCC cases (n = 371) and controls (n = 50) were deposited (| logFC | > 2 and FDR < 0.05 ) based on the TCGA transcriptomic datasets (Fig. 3a,b). 144 un-regulated DEGs and 7 down-regulated DEGs were identi ed, respectively (Table S4). Consistent with the predicted target genes of miR-18a, the KEGG enrichment of up-regulated DEGs presented a latent HCC-promoting feature by harboring signi cant annotations in cell cycle and p53 signaling pathway (P < 0.05) (Fig. 3c). Other characteristics were also provided by GO and KEGG enrichment on DEGs, including some down-regulated liver-speci c metabolism pathways (xenobiotic metabolism, P450 activities, gluconeogenesis, etc) and elevated tumor-promoting features (DNA replication, nuclear division, cell cycle phase transition, etc) ( Fig. 3d ~ f).
Next, the intersection between DEGs compared to predicted target genes was examined and the only gene obtained was Cholinergic Muscarinic Receptor 2 (CHRM2), dramatically (Fig. 4a). The expression level of CHRM2 in HCC tissues was signi cantly lower than in paired adjacent normal tissues (Fig. 4b). The Kaplan-Meier Plotter online tool was used to evaluate the prognostic role of miR-18a and CHRM2 in HCC based on RNA-seq data (Fig. 4c, d). Notably, a high level of miR-18a and low level of CHRM2 were authenticated as prognosticators for poor overall survival (OS) and the results were statistically signi cant (P < 0.01). In order to unravel the possible mechanism of CHRM2 in affecting the HCC prognosis, PPI network of CHRM2 was constructed which contained 11 proteins (enrichment P-value = 3.84e-13). Intriguingly, proteins interacting with CHRM2 indeed engaged PI3K/AKT pathway (GNB1, GNB5, GNG2, CHRM2), Ras signaling pathway (GNB1, GNB5, GNG2) and Pathways in cancer (GNB1, GNB5, GNG2, GNAI1) after KEGG enrichment analysis from STRING database (Fig. 4e). Taken together, these results highlighted a potential impact of miR-18a and CHRM2 on HCC prognosis and indicated a latent value as prognostic biomarkers in HCC tissues.

Discussion
Hepatocellular carcinoma is the most prevalent form of liver cancer and accounts for the fourth most common malignant tumor in China [50]. Bearing the third highest cause of cancer-related mortality, HCC poses a gigantic threat on the public health [51]. Due to the smoldering onset and slowly progressive course, the early manifestations of HCC are di cult to detect which leads to the dilemma that a large proportion of HCC patients cannot be diagnosed properly until the mid-late stage, missing the best period for cures [52]. Moreover, the heterogeneity between HCC patients has been greatly appreciated recently. As displayed in the study of Gao et. al [53], differences in clinical features and prognosis of HCC patients can be nely demonstrated through transcriptomic and proteomic discrepancies. Thus, the individual variability suggests the clear demand for biomarkers more e cacious for HCC prognostic prediction, which is essential for guiding the treatment strategies.
A huge number of miRNAs have been reported to be involved in cancer onset and progression and are differentially expressed among various cancer tissues [54][55][56]. Through participating in various biological processes including modulating oxidative stress, enhancing Warburg effect, blocking apoptosis, etc, miRNAs play crucial roles in cancer onset, invasion and metastasis [57][58][59]. However, it is important to note the truth that miRNAs can serve as both oncogenes and tumor-suppressing genes depending on several parameters such as miRNA types and tumor classi cation[28-31, 34-36]. Owing to the characteristics above, miRNAs and their panels have been receiving ascending insights into their potential as prognosticators[60, 61].
Through literature research, we found that miR-18a has been reported up-regulated in several cancers in previous studies. Yuan et. al validated an activated Wnt/β-catenin signaling pathway mediated by upregulated miR-18a, which elevated gastric cancer metastasis [30]. Similar tumor-promoting progresses corresponding to high miR-18a levels were unveiled in cancers including HCC. The study of Avencia et. al uncovered that up-regulated miR-18a led to a low SOCS5 level, contributing to the HCC tumorigenesis through enhancing the activity and stability of mTOR signaling pathway [34]. However, the underlying mechanism of miR-18a on the HCC prognosis, which is determinant to prognosticator screening, has not yet been revealed clearly. Hence, carrying out bioinformatic analysis on miR-18a might help to address this problem.
Our mature sequence identi cation of miR-18a depicted a well-conserved orthologous sequence UGUUCUAAGGUGCAUCUAGUGCAGAUAGUGAAGUAGA across several species, underscoring the evolutionarily conserved and important functions.
Then by conducting GO analysis, we found that target genes of miR-18a were mainly functioning in the nucleus, and were enriched in biological processes involving cell cycle, transcription modulating, protein However, given the small number of researches directly focused on the association between miR-18a and HCC, we tried to obtain a rough estimate of the possible target gene-depend functions of miR-18a through investigating published studies in PubMed database. As expected, several target genes predicted in our analysis involving SOX6, NOTCH2, RORA, NEDD9 [33,[37][38][39], etc, have been validated to mediate the effects of miR-18a in multiple types of cell lines and tissues mainly through promoting proliferation and migration. Notably, in the study of Zheng et. al, it was demonstrated that miR-18a/NEDD9 axis played a directly promoting role in HCC progression, making the association between miR-18a and HCC more convincible.
For further exploring the value of miR-18a as an in-situ biomarker of HCC, TCGA data were acquired to support our analyses. The down-regulated gene CHRM2 was present in both TCGA DEGs in HCC and predicted target genes of miR-18a. Astonishingly, a signi cant prognostic in uence of the muscarinic receptor, CHRM2, was observed through the Kaplan-Meier Plotter overall survival data in HCC that low level of CHRM2 meant a better prognosis. Actually, the correlation between CHRM2 and cancer has received little study at present. The establishment of the linkage between CHRM2 and non-small cell lung cancer was implemented by Zhao et. al through NF-κB signaling pathway [69]. However, in their study, inactivation of CHRM2 led to a suppressed NF-κB signaling pathway, nally weakened the epithelialmesenchymal transition process, which indicated a tumor-promoting role of CHRM2. To explain this contradiction, proteins interacting with CHRM2 was identi ed and their KEGG enrichment suggested that they might be involved in PI3K/AKT signaling pathway and Ras signaling pathway, giving the possibility for CHRM2 of participating in the miR-18a-mediated cancer-related biological processes because of the targeting of miR-18a and the shared pathway with miR-18a. Indispensably, survival analysis of miR-18a was conducted and high miR-18a level was consistent with a signi cantly poor prognosis, suggesting miR-18a as a potential in-situ prognostic factor, which was comparatively apprehensible based on the results above. Despite the prognostic signi cance of CHRM2 and miR-18a in HCC, we should recognize the limitation of our study that mechanisms of CHRM2 in HCC remain largely unknown, and the procancer effect of miR-18a lacks holistic uncovering. The boundedness of our study needs further in vivo and in vitro veri cations hereafter.

Conclusions
Altogether, miR-18a might be involved in pathways contributing to the onset and progression of HCC. Its target gene CHRM2 is signi cantly down-regulated in HCC. What's more, up-regulated miR-18a and downregulated CHRM2 are indicative of poor prognosis in HCC. We speculated that CHRM2 might participate in the modulation of miR-18a in HCC based on KEGG analysis of related proteins. The expression level of miR-18a could thus be introduced as a biomarker for HCC prognosis. Availability of data and materials The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Ethics approval and consent to participate All data sets were derived from The Cancer Genome Atlas (TCGA) data portal and are publicly accessible from the TCGA website.