Differentially Expressed miRNAs in Hepatocellular Carcinoma Target Genes in the Genetic Information Processing and Metabolism Pathways

To date, studies of the roles of microRNAs (miRNAs) in hepatocellular carcinoma (HCC) have either focused on specific individual miRNAs and a small number of suspected targets or simply reported a list of differentially expressed miRNAs based on expression profiling. Here, we seek a more in-depth understanding of the roles of miRNAs and their targets in HCC by integrating the miRNA and messenger RNA (mRNA) expression profiles of tumorous and adjacent non-tumorous liver tissues of 100 HCC patients. We assessed the levels of 829 mature miRNAs, of which 32 were significantly differentially expressed. Statistical analysis indicates that six of these miRNAs regulate a significant proportion of their in silico predicted target mRNAs. Three of these miRNAs (miR-26a, miR-122, and miR-130a) were down-regulated in HCC, and their up-regulated gene targets are primarily associated with aberrant cell proliferation that involves DNA replication, transcription and nucleotide metabolism. The other three miRNAs (miR-21, miR-93, and miR-221) were up-regulated in HCC, and their down-regulated gene targets are primarily involved in metabolism and immune system processes. We further found evidence for a coordinated miRNA-induced regulation of important cellular processes, a finding to be considered when designing therapeutic applications based on miRNAs.


Results
The workflow in Fig. 1 was employed to elucidate the pathways and clinically relevant gene targets that are modulated by dys-regulated miRNAs in 100 HCCs. The steps of the workflow were as follows: (1) Identify differentially expressed miRNAs between tumorous and adjacent non-tumorous liver from the miRNA expression profiles.
(2) Integrate miRNA and mRNA profiles to select potentially relevant miRNAs in HCC. (3) The pathways and clinical significance of these miRNAs and their gene targets were examined. miRNAs and gene targets differentially expressed in HCC. We assessed the levels of 829 human mature miRNAs. We found thirty-two to be significantly differentially expressed (absolute fold change (FC) > 1.5 and FDR < 0.05) between the tumorous and adjacent non-tumorous liver of HCC patients (Fig. 2). Six of these were up-regulated while 26 were down-regulated in the tumors.
To identify the relevant gene targets of these differentially expressed miRNAs in HCC, we employed an integrated approach, which combines miRNA with gene expression profiles as well as in silico prediction of likely gene targets of miRNAs. As miRNAs modulate gene expression through mRNA degradation or translational repression, we first examined a negative correlation in expression between miRNA and their gene targets.
We determined the negative correlation between the expression of each of the differentially expressed miRNA and all genes in the human genome and ranked genes according to the strength of the relationship. An example of a strong negative correlation between miR-93 and PTPRN2 is shown in Supplementary Fig. S1. We then compared this ranked list of negative gene correlation to the set of in silico targets predicted by the MiRanda program 17 using gene set enrichment analysis (GSEA) 18 . GSEA computes an enrichment score by stepping through the ranked list of genes. If a gene is an in silico predicted target, the enrichment score is increased. Otherwise, the enrichment score is decreased. An example of an enrichment plot for miR-93 is shown in Supplementary Fig. S2. We found eight of the 32 differentially expressed miRNAs are significantly enriched with negatively correlated predicted gene targets (FDR < 0.01; Table 1, MiRanda column). Interestingly, most of these miRNAs were previously reported to be associated with HCC 3 . We also compared the same ranked list of negative gene correlation with the set of in silico targets predicted by MirTarget2 19 , which identified seven miRNAs with statistically significant enrichment scores (Table 1, MirTarget2 column). Six of these miRNAs (miR-21, miR-26a, miR-93, miR-122,  (Table 1) and were subjected to further analysis.
Pathways affected by relevant miRNA targets. We first examined if any of the six abovementioned miRNAs are coordinately dys-regulated in the tumors of HCC patients using Spearman's correlation. As evident in Supplementary Table S2, only a pair of miRNAs (miR-93 and miR-221) that is up-regulated and a pair of miR-NAs (miR-26a and miR-130a) that is down-regulated in the tumors of HCC patients are correlated (Spearman's correlation coefficient: 0.6 and 0.56, respectively), suggesting that most of these miRNAs are likely regulated independently in HCC patients.
To elucidate the pathways that are likely to be affected by the relevant gene targets of the three up-regulated and three down-regulated miRNAs, we first examined the differentially expressed gene targets of these differentially expressed miRNAs. Fewer differentially expressed gene targets were identified for down-regulated miRNAs than for up-regulated miRNAs (Fig. 3). The list of differentially expressed gene targets of down-regulated and up-regulated miRNAs, their correlation coefficient, and their MiRanda and MirTarget2 predictions are given in Supplementary Tables S3 and S4. For down-regulated miRNAs, out of 148 gene targets, few genes (6, 15, and 4 gene targets respectively) are common between two miRNAs (Fig. 3). On the other hand, for up-regulated miR-NAs, out of 517 gene targets, a larger number of genes (46,27, and 44 gene targets respectively) were common between two and 14 between all three miRNAs. These data suggest that while there is a subset of gene targets that   1 Only miRNAs that showed statistically significant enrichment against both target prediction programs (MiRanda and MirTarget2, GSEA FDR < 10 -3 ) were considered for further analysis.
are commonly modulated by two or more significantly differentially expressed miRNAs, the majority of the gene targets of these miRNAs are unique. Pathways modulated by relevant gene targets of differentially expressed miRNAs were determined using KEGG pathways as defined in MSigDB 20 . Interestingly, although the proportion of differentially expressed gene targets that are common among either the three up-regulated or the three down-regulated miRNAs is relatively low (Fig. 3), all three up-regulated miRNAs (miR-21, miR-93, and miR-221) coordinately target a set of pathways that is distinct from all three down-regulated miRNAs (miR-26a, miR-122, and miR-130), which also coordinately target a set of pathways ( Table 2). We found only the Cell Growth and Death pathway is targeted by all six miRNAs ( Table 2). A closer examination of the Cell Growth and Death pathway revealed that the gene targets of all six differentially expressed miRNAs modulate TP53 signaling. However, the gene targets of all three up-regulated miRNAs are involved in apoptosis, whereas the gene targets of all three down-regulated miRNAs are involved in cell cycle (Table 3). This is consistent with tumorigenesis down-regulating apoptosis and up-regulating cell proliferation.
Gene targets of all three up-regulated miRNAs (miR-21, miR-93, and miR-221) were found to coordinately target primarily metabolism pathways, particularly the amino acid, carbohydrate, and lipid metabolism pathways, as well as the immune system ( Table 2). On the other hand, gene targets of all three down-regulated miRNAs (miR-26a, miR-122, and miR-130) primarily target Genetic Information Processing, in particular, DNA replication and repair, transcription, cell growth, and nucleotide metabolism ( Table 2). This is consistent with the process of tumorigenesis in which pathways not immediately pertinent to the proliferation of the tumor cell, such as amino acid, carbohydrate, or lipid metabolism, are down-regulated so resources can be redirected to ensure the proliferation of the tumor cells.
Next, we derived gene networks based on differentially expressed targets using the Ingenuity Pathway Analysis (IPA) program (http://www.ingenuity.com). We found that the network of genes in the metabolism pathways, dys-regulated by the three up-regulated miRNAs, is primarily centered on nodal molecules. IL6 and CXCL8 are nodal molecules of metabolism pathways, both of which are targets of miR-93 (Fig. 4a). We also found that CXCL8 is an important nodal molecule, together with TP53 in the gene network of the inflammation and immune response pathways (Fig. 4b). Moreover, TP53 and TP73 play central roles in the network of genes in the Cell-Cycle Progression and Apoptosis pathways dys-regulated by all six miRNAs (Fig. 4d). The gene network dys-regulated by the three down-regulated miRNAs in the Genetic Information Processing pathways including Cell Proliferation, DNA Replication, RNA transcription, and DNA metabolism centers on PTEN/AKT, Histone H3, and Cyclin A/A2 (Fig. 4c).

Relevant differentially expressed miRNA targets associated with various clinical characteristics.
We next evaluated if differentially expressed miRNAs or their relevant differentially expressed gene targets are associated with various clinical characteristics, such as gender, age, hepatitis infection status, liver cirrhosis, stage, grade, size, relapse status, encapsulation, and Alpha Fetoprotein (AFP) levels (Tables 4 and 5). As evident in Table 4, miR-122 showed significant association with AFP levels (FDR < 0.2). Similarly, miR-221 was significantly correlated with age (FDR: 0.13) and hepatitis infection (FDR: 0.15). More specifically, HBV-positive patients showed a significantly higher miR-221 expression compared to HBV-negative patients. Next, we associated the differentially expressed targets of the six miRNAs with various clinical characteristics. While only a single differentially expressed gene target of down-regulated miR-122 was found to be significantly associated with stage of HCC, 13, 2, and 5 down-regulated gene targets of up-regulated miRNAs were found to be significantly associated with stage of cancer, tumor grade, and encapsulation. No miRNA was significantly associated with relapse status using the Cox proportional hazards model (Supplementary Table S9), although miR-130a showed the lowest p-value. When another cohort of 373 HCC patients from The Cancer Genome Atlas (TCGA) 21 was assessed, miRNA-130a was found to be significantly associated with patient outcome (time-to-death, 89 events) (HR: 0.52, FDR: 0.02, Supplementary Table S9 and Figure S4).
We employed IPA to obtain a glimpse of the major cancer-related pathways dys-regulated by these genes associated with various clinical characteristics. Interestingly, greater than 75% of the genes dys-regulated by up-regulated miRNAs that are associated with a lower stage of cancer, tumor grade, or encapsulation were found to be in cancer or related pathways (primarily cell cycle, cell death and survival, as well as inflammatory/immune disease/response).

Discussion
Through profiling, miRNAs have been shown to be differentially expressed in the tumors of HCC patients (see review 5 ). Gene targets of many of these dys-regulated miRNAs have been identified. However, thus far, only very few studies have attempted to identify which one/combination of the differentially expressed miRNAs is more likely to play roles in HCC through miRNA-gene interaction networks.
In this study, we integrated miRNA with gene expression profiles, as well as miRNA target prediction programs to identify a subset of miRNAs that are highly correlated with appropriately dys-regulated predicted gene targets. We hypothesized that this subset of miRNAs is most likely to play roles in HCC through their gene targets. We also examined the correlation of these miRNAs or their gene targets with clinical characteristics, as well as pathways that are perturbed by the appropriately dys-regulated predicted gene targets of these miRNAs.
Of the 32 miRNAs that were differentially expressed between the HCC tumorous and non-tumorous livers (Fig. 2), we found that six are highly correlated with appropriately dys-regulated predicted gene targets (Fig. 3). Four (miR-21, miR-221, miR-26, and miR-122) of these six miRNAs were previously reported to be differentially expressed in HCC in several studies and to dys-regulate important cancer pathways 12,[22][23][24] , highlighting the importance of these miRNAs in HCC. We further found these six miRNAs are mainly independently regulated in HCC patients except for two pairs of miRNAs (miR-93 with miR-221 and miR-26a with miR-130a), which are moderately correlated in expression (Supplementary Table S2).
Three (miR-26a, miR-122, and miR-130a) of these six miRNAs were down-regulated while the other three miRNAs (miR-21, miR-93, and miR-221) were up-regulated in the tumors of HCC (Fig. 2). Interestingly, the total number of down-regulated predicted gene targets of the three up-regulated miRNAs (517) was ~3.5 times more Organismal Systems   Table 3. Detail of Pathway group "Cell growth and death" enriched in miRNA Targets.
than the up-regulated gene targets of the three down-regulated miRNAs (148). In addition, the greater percentage of common gene targets was found among the up-regulated miRNAs (34-52%) ( Fig. 3 and Supplementary  Table S4) than among the down-regulated miRNAs (28-31%) ( Fig. 3 and Supplementary Table S3). This perhaps suggests that the up-regulated miRNAs play a more diverse but coordinated role in HCC, down-regulating the expression of a greater number of different gene targets. Some miRNAs showed significant association with clinical characteristics. Expression of miR-122 was significantly correlated with AFP levels. Indeed, miR-122 was previously implicated in the regulation of AFP and high AFP levels are associated with a more aggressive cancer 25 . Similarly, expression of miR-221 was significantly correlated with age and hepatitis infection. Furthermore, several of the predicted dys-regulated gene targets of these miRNAs were associated with clinical characteristics. The down-regulated gene targets of up-regulated miRNAs were found to be associated with a better prognosis, including a lower stage of cancer and lower grade of tumor, as well as complete encapsulation, and are in cancer or cancer-related pathways, including cell cycle, cell death, and survival as well as inflammatory/immune disease/response. On the other hand, the sole up-regulated gene targets of the down-regulated miR-122 that showed clinical relevance was found to be associated with a higher stage of cancer.
Notably, although the majority of gene targets of the six miRNAs were unique, pathway enrichment analysis revealed that relevant gene targets of all six up-regulated miRNAs (miR-21, miR-93, and miR-221) were enriched in metabolism and immune system processes ( Table 2) while relevant gene targets of all three down-regulated miRNAs (miR-26a, miR-122, and miR-130a) were enriched in genetic information processing pathways including cell proliferation, DNA replication, and metabolism, as well as RNA transcription. These observations suggest built-in redundancy of different gene targets regulated by different miRNAs targeting the same set of pathways. This redundancy may pose challenges for any single miRNA to be used as a target for cancer therapy.
The pathways targeted by these six miRNAs represent the various hallmarks of cancer. The predicted gene targets of the three up-regulated miRNAs reside in the metabolism pathway, which has recently generated great interest [26][27][28] . In fact, these up-regulated miRNAs (miR-21, miR-93, and miR-221) have previously been implicated in the regulation of lipid metabolism [29][30][31][32] . Closer examination of the pathways leading to the dys-regulation of lipid, amino acid, and carbohydrate metabolism identifies proinflammatory molecules IL6 and CXCL8, both of which are down-regulated by miR-93, as the two nodal molecules that play key roles in metabolism (Fig. 4a). Interestingly CXCL8 was also found to be one of the two key nodal molecules, together with TP53 in the inflammation and immune response pathways (Fig. 4b). As a result, the down-regulation of these pro-inflammatory molecules, IL6 and CXCL8 in cancer cells could perhaps perturb their intrinsic metabolic pathways, as well as help these cells evade the immune response, which potentially leads to a survival advantage of these malignant cells 33,34 .
In addition to being a nodal molecule for inflammation and immune response, TP53, together with TP73, were identified as key molecules involved in the cell-cycle progression and apoptosis pathways. Targets of all six   reported to be involved in TP53-mediated apoptosis, suggesting that miR-93 plays a role in the evasion of apoptosis, a hallmark of cancer 39 . As mentioned earlier, the three down-regulated miRNAs (miR-26a, miR-122, and miR-130a) target cell proliferation and genetic information processing including DNA replication, transcription, and nucleotide metabolism (Table 2 and Fig. 4). Consistent with our results, miR-26a 40 and miR-122 41 were previously reported to inhibit cell growth in liver cancer cells, probably through dys-regulation of DNA replication 42 and metabolic processes 43 respectively. Closer examination reveals that pathways modulated by these down-regulated miRNAs center around key molecules PTEN/AKT and Histone H3. Both miR-26a 44 and miR-122 45 have been reported to target the AKT pathway in HCC. Furthermore, miR-130a was reported to target MET 46 and to stimulate apoptosis in non-small cell lung cancer 46 . As MET also plays a role in AKT signaling 47 , down-regulation of these three miRNAs (miR-26a, miR-122, and miR-130a) may perturb the AKT axis, leading to the disruption of cell death.
Indeed, our data confirms that MET expression was found to be ~1.26 fold higher in the tumorous compared to the adjacent non-tumorous tissues and was moderately correlated with miR-130a (rho = 0.21, Supplementary  Table S5). In addition, the involvement of Histone H3 as a nodal molecule suggests that these pathways may be epigenetically regulated. The gene targets of these miRNAs have been reported to either be direct binding partners of Histone H3 (NCAPD2 48 , CHAF1A 49 , and PARP1 50 ) or activated by post-translational modification of Histone H3 (CCNE1 51 , Cyclin A 52 CDKN2B 53 , TACC2 54 ), highlighting the complex gene networks regulated by these miRNAs.  In summary, these data reveal six miRNAs that may play important roles in tumorigenesis through the dys-regulation of gene targets. These six miRNAs generally modulate unique gene targets but coordinately target a few common pathways, with partially overlapping key nodal molecules, highlighting the redundancy of miRNAs and gene targets that focuses on reprogramming critical pathways to facilitate tumor development. Although the redundancy poses a challenge to the identification of a single target for therapy, deeper analyses of the pathways modulated by these miRNAs suggest that key nodal molecules such as IL6, CXCL8, TP53, TP73, PTEN, and AKT may serve as potential targets for intervention. Many of these molecules have already been explored as targets of therapy in other cancers 33,[55][56][57][58] and future focus could be directed towards these molecules in HCC.
In addition to potential targets for intervention, this analysis also highlights several dys-regulated predicted gene targets as potential biomarkers for cancer prognosis, since some of the down-regulated predicted gene targets of the up-regulated miRNAs were found to be associated with better disease prognosis, while the sole up-regulated gene target of the down-regulated miRNA was associated with poorer prognosis.  Table 5.

Methods
Gene and miRNA expression profiling. We sent samples to Miltenyi Biotec, CA, for gene and miRNA expression microarray profiling. Briefly, total RNA isolation for tumorous and adjacent non-tumorous tissue samples was performed using the standard RNA extraction protocols (Trizol). We assessed total RNA integrity using the Agilent 2100 Bioanalyzer platform (Agilent Technologies), and only samples with RNA integrity number (RIN) > 6 were profiled for gene and miRNA expression 59 .
For gene expression profiling, we used 1 μ g of each total RNA sample as starting material. We amplified and labeled total RNA samples using the Agilent Low RNA Input Linear Amp Kit (Agilent Technologies) following the manufacturer's protocol. We measured cRNA yields and dye-incorporation rates with the ND-1000 Spectrophotometer (NanoDrop Technologies). We labeled non-tumorous samples with Cy3 and tumorous samples with Cy5. We performed hybridization according to Agilent's 60-mer oligo microarray processing protocol using the Agilent Gene Expression Hybridization Kit (Agilent Technologies). Briefly, we combined and hybridized 825 ng of the corresponding Cy3-and Cy5-labeled fragmented cRNA overnight (17 hours, 65 °C) to Agilent Whole Human Genome Oligo Microarrays 4 × 44K (G4112F). We washed microarrays with 6x SSPE buffer containing 0.005% N-lauroylsarcosine for 1 min at room temperature followed by a second wash with pre-heated 0.06x SSPE buffer (37 °C) containing 0.005% N-lauroylsarcosine for 1 min. The second washing step was performed with acetonitrile for 30 seconds. Finally, we detected fluorescence signals using Agilent's DNA microarray scanner (Agilent Technologies).
For miRNA profiling, we labeled samples according to the undisclosed miRXplore ™ user manual. We hybridized fluorescent-labeled samples to miRXplore ™ microarrays using the a-Hyb ™ hybridization station. We labeled non-tumorous samples with Hy3 and tumorous samples with Hy5. Finally, we detected fluorescence intensity signals using an Agilent laser scanner (Agilent Technologies).
Pre-processing gene expression profiles. We loaded raw data of gene expression profiles into R using the R/Bioconductor Limma package 60 . We corrected data for background noise (normexp) 61 , followed by within (loess) and between array (quantile) normalization 62 . We detected no batch effects between expression profiles using principal component analysis (PCA).
Pre-processing miRNA expression profiles. We loaded the raw miRNA expression profiles into R using Limma. Similar to mRNA gene expression profiles, we background corrected and normalized miRNA raw data within and between arrays. Next, we corrected arrays for batch effects using ComBat 63 . We determined differentially expressed miRNAs between tumorous and adjacent non-tumorous tissues of 100 HCC patients using Limma. To account for multiple testing, we employed the Benjamini-Hochberg method to calculate FDR values 64 . MiRNAs with an absolute fold change greater than 1.5 and an FDR smaller than 0.05 were considered significantly differentially expressed.
Targets of differentially expressed miRNAs. For every differentially expressed miRNA, we calculated the negative Spearman's correlation coefficient between the miRNA and all genes in expression across tumorous and adjacent non-tumorous tissues. We ranked genes by their negative correlation with the miRNA. Next, we downloaded 65 and used miRNA target predictions based on the MiRanda target prediction program as a gene set for enrichment analysis. For each miRNA, we detected enrichment of its predicted targets in the ranked list of genes through GSEA's pre-ranked function (GSEA version 2.0.10) 18 . Similarly, we calculated enrichment scores for target predictions based on the MirTarget2 program 66 . Only miRNAs with a FDR < 10 -3 in both programs, MiRanda and MirTarget2, were considered for further analysis. Next, we calculated target prediction scores for all genes using the MiRanda target prediction command line tool (version 3.3a) 17 . For that purpose, we downloaded mature miRNA sequences from miRBase 67 (release 19) and 3'UTR exon sequences from the UCSC Scientific RepoRts | 6:20065 | DOI: 10.1038/srep20065 genome browser database 68 (hg19). We normalized prediction scores to zero by subtracting the minimum score achieved by any gene. Next, we calculated the product between the negative correlation with the miRNA and the normalized target prediction score and used it as a ranking metric (Supplementary Tables S5 and S6).
Finally, based on this ranking metric, we identified biological themes for each miRNA using GSEA′ s pre-ranked function. The KEGG (excluding human disease) gene sets from MSigDB were employed in our analyses. Then we summarized enrichment results by consolidating pathways to pathway groups and pathway groups to categories based on the KEGG orthology conceptual hierarchy, implemented in the Bioconductor keggorthology package (version 2.1.0). We ranked pathway groups by their total number of underlying pathways, and we considered up-and down-regulated miRNAs separately. For ties, we ranked groups with at least one pathway enriched in all three up-regulated miRNAs, down-regulated respectively, and then more significant FDRs higher. We rearranged categories according to the highest-ranking pathway group in up-regulated miRNAs.
Gene networks with functional annotation. We collected genes for each pathway group across three up-regulated miRNAs and three down-regulated respectively. We only considered differentially expressed genes of the GSEA core enrichment set. IPA was employed to infer regulatory networks and functional annotation for differentially expressed genes associated with selected pathways in each category in Table 2.
Association of differentially expressed miRNAs and their targets with clinico-pathological features. MiRNAs and their relevant MiRanda predicted targets with significant differential expression (absolute FC > 1.5, FDR < 0.05) were associated with clinical features, namely gender, age, hepatitis infection status, liver cirrhosis, stage, grade, size, relapse status, encapsulation, and Alpha Fetoprotein (AFP) levels. We evaluated association through the log-rank test and Cox's proportional hazards model (relapse status) or the Kruskal-Wallis rank sum test (all other clinical parameters). We obtained the genes' functional annotation through IPA.