Dynamic recruitment of microRNAs to their mRNA targets in the regenerating liver

microRNAs (miRNAs), 22–23 nucleotide noncoding RNAs, contribute to the control of diverse developmental, growth, and disease processes [1, 2]. Abnormal expression of miRNAs has been established in cancers, with miRNAs functioning as either tumor suppressors or oncogenes. MicroRNAs decrease expression of mRNA targets by either destabilization of mRNA or inhibition of protein translation [3]. MicroRNAs are thought to target mRNAs through binding of nucleotides at position 2–8 of the miRNAs (the so-called ‘seed region’) to a complementary sequence in the mRNA [4, 5].

While differential expression of miRNAs has been determined in multiple contexts, the validation of physiologic miRNA targets has proven to be difficult. Modulation of miRNA levels using gain- and loss-of-function approaches have identified some mRNA targets; however, concerns exist regarding the potential for indirect and off-target effects, particularly in the case where a miRNA is overexpressed. Currently available miRNA target prediction algorithms produce large numbers of potential targets; however, very few of these targets have been experimentally validated. Many current algorithms mine sequences limited to the 3^′UTR of mRNAs, and therefore do not identify miRNAs that target the coding region and 5^′UTR, even though miRNA targeting to all mRNA regions has now been experimentally validated [5–11]. Furthermore, the ability of miRNAs to target mRNAs with only partial complementarity [12] indicates that in many cases identification of miRNA targets may not be possible based on seed match complementarity alone [10]. Consequently, better methods are needed to identify bona fide miRNA targets.

Recently, Pandolfi and colleagues have proposed an ‘mRNA code’, in which competing endogenous mRNAs (ceRNAs) communicate with each other through crosstalk between shared microRNA regulatory elements (MRE) [10, 13]. An important prediction of this model is that competition for a common pool of miRNAs could result in de-repression of one mRNA following robust activation of a second mRNA containing the same MREs. Therefore, identification of MREs and mRNA:miRNA targeting relationships is an essential step towards understanding the role of miRNAs in the regulation of complex biological responses.

The partial hepatectomy model in rodents has been used extensively to investigate the mechanisms responsible for hepatic growth and proliferation, and is the best in vivo model of synchronous cell cycle progression in mammals [14, 15]. Following surgical removal of two thirds of the liver, hepatocytes and nonparenchymal cells rapidly reenter the cell cycle, replicate, and restore the original mass of the liver within 10–14 days [14]. This process involves a complex regulatory cascade of cytokine signals and transcriptional regulators that coordinate cell cycle progression while maintaining homeostasis [15, 16]. An essential contribution of miRNAs in this regenerative response has been supported by a recent study in which mice with genetic deletion of the DROSHA cofactor DGCR8, a factor required for microRNA biogenesis, exhibited markedly impaired hepatocyte proliferation after partial hepatectomy [17]. Although changes in expression of miRNAs after partial hepatectomy and in liver-graft models have been reported using array-based assays [17–24], only a small number of targets have been validated [22, 23, 25, 26]. Importantly, the identification of the subset of mRNAs that are regulated by miRNAs in the regenerating liver is far from complete, in part due to the large number of possible mRNA:miRNA targeting relationships predicted by computational approaches.

We hypothesized that recruitment of miRNAs to the RISC, or RNA-induced silencing complex, which contains the partially base-paired miRNA and its mRNA target, would be a more informative method to assess miRNA activity during liver regeneration, since changes in the overall expression level of miRNAs reported thus far do not take into consideration changes in target mRNA expression or MRE accessibility. To identify which miRNAs participate in the regenerative process and to identify their targets, we used a technique pioneered by the Darnell laboratory in which UV cross-linked miRNA:mRNA complexes are immunoprecipitated with an antibody to Argonaute, an essential component of the RISC, and then subjected to deep sequencing analysis [27–29]. We found that a subset of miRNAs is dynamically recruited to the RISC during liver regeneration, and that for the majority of these miRNAs, RISC recruitment did not correlate with changes in overall expression. Furthermore, we exploited the comparison of RNA-seq and RISC recruitment data to identify those mRNAs whose recruitment to the RISC was highly enriched relative to their overall abundance, and which are therefore likely to be regulated by miRNAs.

Changes in expression of miRNAs, mRNAs and/or accessibility of microRNA regulatory elements (MREs) to their complementary miRNAs are likely to affect miRNA recruitment to the RISC [4]. We hypothesized that miRNA activity would be more accurately determined by quantification of miRNA abundance in the RISC than by changes in miRNA levels alone. To this end, we determined miRNA recruitment at selected time points posthepatectomy previously shown to correspond to the G₁ (1h), S (36H) and M (48H) phases of the hepatocyte cell cycle, and in the quiescent liver, corresponding to G₀[30, 31]. Cell cycle stage was confirmed by activation of genes associated with G₁ (jun, fos, myc), S (cyclin D1) and M (Foxm1) cell cycle phases determined by RNA-seq analysis of samples from these time points (Additional file 1: Figure S1). Following UV-crosslinking of the RISC to microRNA and mRNA, we quantified miRNAs immunoprecipitated with an antibody that recognizes Argonaute 1 through 4 proteins and therefore immunoprecipitates all RISCs [32] by ultra-high throughput sequencing. We calculated the loading microRNA relative to all others for each replicate as reads per million (RPM). We then employed a stringent cut-off of 100 RPM through which we detected 226 miRNAs in the RISC complex at one or more time points examined (Additional file 2: Table S1 contains loading data for the 632 miRNAs that had any aligning reads). An important issue for any comparison across time is global normalization of the data. Because the average levels of RISC-associated microRNAs were relatively unchanged over time, we employed additional global normalization, as well as a stricter 1,000 RPM cut-off, to identify differentially-loaded microRNAs.

Using these criteria, 16 miRNAs exhibited a significant increase (9) or decrease (7) in RISC recruitment at one or more time points posthepatectomy relative to the quiescent liver, with a dynamic range exceeding five-fold in an interval as short as one hour (Figure 1A). We did not detect differences in RISC recruitment of mir-21-5p, a microRNA that had been shown previously to be significantly induced 18 hours posthepatectomy [17]. This may reflect a lack of correspondence between overall expression of mir-21-5p and RISC recruitment. It is also conceivable that mir-21-5p RISC enrichment is increased selectively at 18 hours posthepatectomy, a time point that was not examined in the current study. Our HITS-CLIP analysis did, however, predict several previously confirmed mir-21-5p targets in the regenerating liver including the anti-apoptotic factor, Btg2 [17] and the NF-κB inhibitor Pellino 1 [18].

The change in recruitment was extremely rapid for the majority of these miRNAs. For example, RISC recruitment of mir-215 increased approximately 5.6-fold by 1 h posthepatectomy, and was also similarly enriched at the time point corresponding to hepatocyte S phase (36 h) relative to quiescent liver. Surprisingly, mir-215 and mir-192, which contain the same seed sequence (Additional file 3: Figure S2), are both regulated by p53 [33, 34] and expressed at similar levels after hepatectomy (Figure 1B), followed different patterns of RISC recruitment. This discordance could reflect requirements for sequences outside of the seed region of the two miRNAs for binding to their targets or interaction with RNA binding proteins.

This rapid change in miRNA recruitment to the RISC is unlikely to be due solely to an increase in miRNA transcription, and likely reflects enhanced miRNA processing, stability, changes in miRNA subcellular localization, and/or changes in the abundance of specific mRNA targets or RNA binding proteins. Several miRNAs exhibited changes in RISC recruitment at later time points posthepatectomy. For example, significant RISC enrichment of mir-142-3p occurred 48 h posthepatectomy. This miRNA has been shown recently to inhibit RAC1-mediated colony formation migration and invasion in HCC cell lines [35] and could have a similar function in proliferating hepatocytes. We investigated whether changes in miRNA recruitment to the RISC were associated with corresponding changes in overall miRNA abundance during the regenerating time course (Figure 1A, B). Importantly, with the exception of mir-31, mir-144 and mir-378, we found no correlation between miRNA expression and RISC recruitment and for many of the miRNAs examined, changes in RISC recruitment and overall expression were inversely correlated. Overall, the dramatic changes in the levels of actively engaged miRNAs over a short time frame emphasize the necessity to determine miRNA recruitment to the RISC rather than relying on simple determination of miRNA abundance in order to fully assess the impact of miRNA regulation.

Next, we turned our attention to the mRNAs that were identified in the RISC using our HITS-CLIP assay. We speculated that miRNA regulation would be more significant for those mRNAs that were highly enriched in the RISC relative to their overall abundance in the tissue. Therefore, to identify ‘expression-weighted footprints’, we used RNA-seq to quantify total mRNA expression in quiescent liver, and 1 h, 36 h and 48 h posthepatectomy. Next, we calculated RISC enrichment for all mRNAs relative to their overall abundance (Additional file 4: Table S2). Using the enrichment-weighted values for all mRNAs, we applied k-means cluster analysis to identify groups of genes that exhibited similar changes in RISC recruitment across the posthepatectomy time course. Pathway analysis applied to these gene clusters identified specific networks and biological functional categories that were highly significant (Figure 2A, B and Additional file 5: Table S3 and Additional file 6: Table S4). The overall mRNA expression levels follow the expected expression patterns, with genes involved in basic liver function pathways enriched in the quiescent liver (e.g. clusters 5, 11), and cell cycle genes increasing at 48 hours posthepatectomy (e.g. clusters 3, 8). Cell cycle progression and checkpoint control gene clusters were maximally enriched in the RISC (B) at 36 h (cluster 3: mitosis p < 1.92E-10) and 48 h (cluster 12, DNA replication checkpoint, p < 2.71E-07). In contrast to the pattern of cell cycle gene enrichment in the RISC, genes involved in amino acid metabolism (p < 2.79E-07), lipid metabolism (Cluster 10, synthesis of lipids, p < 3.41E-07), and cell growth (Cluster 10, p <6.45E-08) were gradually decreased in the RISC following hepatectomy. Growth factors expressed in hepatocytes during regeneration including Fgf1 and VegfA were enriched the RISC (Additional file 4: Table S2) suggesting that these factors are also regulated by miRNAs. HGF, which is synthesized by hepatocyte stellate and endothelial cells, was also highly enriched in the RISC in quiescent liver followed by a rapid and sustained fall in enrichment that did not normalize until 48 h after partial hepatectomy (Additional file 4: Table S2). Mobilization of preexisting HGF from extracellular matrix by urokinase plasminogen activator provides a rapidly available source of this growth factor that is required to stimulate hepatocyte proliferation [36]. Our findings suggest that relief of miRNA inhibition of HGF may serve to further enhance HGF mRNA levels during the regenerative response.

Regulation of mRNAs by miRNAs may either “fine-tune” or significantly inhibit target mRNA expression [4]. To investigate the potential regulatory function of miRNAs for their targets in the regenerating liver, we plotted total gene expression and RISC recruitment for two clusters of genes known to be important in the regenerative process after partial hepatectomy [37]. The total regulatory load (TRL) of Gas1, an inhibitor of proliferation, increases coincident with a decrease in its overall expression, suggesting that its inhibition by miRNAs may facilitate the G₁/S transition [38] (Figure 2C). SPHK2 represses transcriptional activation of the immediate-early growth gene Fos via inhibition of histone acetylation [39]. These findings suggest that Fos and SPHK2 inhibition by miRNAs could cooperatively contribute to the resolution of Fos activation posthepatectomy. In contrast, the TRL of mRNAs encoding proteins involved in DNA checkpoint regulation increases at 36 hours posthepatectomy, prior to the observed overall increase in mRNA at 48 hours posthepatectomy (Figure 2D). This temporal relationship suggests that miRNAs may limit and/or delay the increase in expression of replication checkpoint genes during hepatocyte S phase, allowing for completion of DNA replication.

The majority of miRNA target prediction algorithms utilize complementary seed match pairing to the 3^′UTR of candidate mRNAs, although recent reports have identified miRNA regulatory sequences in coding region (CDS) and 5’UTR regions [5–11]. Most of the Ago footprints aligned to sequences in the 3^′UTR of mRNAs, but we also found footprints in the CDS and 5’UTR in a subset of mRNA targets (Figure 3 and Additional file 7: Table S5). Footprints were identified for many genes with known roles in liver growth. For example, we identified a dominant footprint in the 3^′UTR of Ptp4a1 (PRL-1), an immediate-early gene in the regenerating liver associated with cell growth (Figure 3A) [40]. A unique footprint was found in the coding region of the histone deacetylase, HDAC8, a modulator of estrogen related receptor alpha activity (Figure 3B) [41]. An example of a strong 5^′UTR footprint is shown for thyroid receptor beta (THRB), a regulator of hepatic lipogenesis (Figure 3C) [41, 42].

The mRNA fragments identified by HITS-CLIP are predicted to contain or be adjacent to miRNA regulatory elements (MREs), and therefore should enable us to refine the computational predictions of miRNA-mRNA pairs by reducing the ‘search space’ in which to evaluate matches to miRNA seed sequences. For this purpose, we began by compiling all mRNA fragments we obtained from the Argonaute immunoprecipitation, and aligned them to all RefSeq mRNAs. Next, we catalogued all starting positions of these mRNA fragments. These mRNA fragments were then coalesced into ‘footprints’ by the locally strongest accumulation of reads, yielding 472,474 footprints, or sites of RISC occupancy in mRNAs. Finally, these footprints were intersected with all predicted mRNA:miRNA targeting relationships obtained by miRanda, when applied to all known miRNAs and RefSeq mRNAs. This refinement of the computational approach using our experimental data on RISC footprints resulted in 125,949 unique and high confidence miRNA/mRNA pairs, which is a greater than 20-fold refinement of the results obtained from computational target prediction alone (Figure 4).

To assess the accuracy of our predicted miRNA/mRNA pairs, we compared our results to a set of genes whose expression was altered in liver in response to miR-29a-3p antisense oligonucleotide treatment [44]. Hand and colleagues identified 78 genes that were up-regulated upon blocking miR-29a-3p, and an additional 65 genes that were down-regulated. We identified 381 predicted targets of miR-29a-3p in quiescent liver using the methods summarized in Figure 4. Of these predicted targets, 10 (p = 4.16e-06; Fisher exact test, followed by Benjamini-Hochberg correction for multiple testing) were found in the up-regulated set in Hand and colleagues [44] whereas only 3 (p = 0.228) were found in the down-regulated set. Performing the same analysis for all miRNAs, we found the targets of miR-29a-3p to have the second most significant overlap with the up-regulated set obtained experimentally [44]. The most significant overlap was with miR-328-3p that had a similar overlap (10 genes in common with 282 targets and 4 targets in common with miR-29a-3p). Next, we performed a similar comparison of the miR-122-5p targets identified in the quiescent liver [45, 46], and obtained a similar trend. We identified 748 target genes and 866 target mRNAs for miR-122-5p. The overlap between our targets and the microarray data were significant for the genes up-regulated upon interference with miR-122-5p activity (36 of 206 transcripts p = 3.87e-09 for Esau and 42 of 363 p = 3.69e-15 for Krutzfeldt), and not significant for the down-regulated genes (32 of 619 p = 1 for Esau and 16 of 305 p = 1 for Krutzfeldt). For both sets of up-regulated genes, miR-122-5p was the miRNA with the most significant overlap with the microarray data. Based on the differences between the experimental models, our results are likely to underestimate the true extent of mir-29a-3p or miR-122-5p target overlap. Nevertheless, these findings provide strong support for the validity of our target prediction approach.

A prediction of the ceRNA hypothesis is that mRNAs containing the same MREs will compete for a common pool of miRNAs [13]. We postulate that activation of a subset of growth regulatory mRNAs in the regenerating liver may be due in part to relief of inhibitory effects of other highly induced mRNAs that compete for the same pool of miRNAs. At present the extent of these ceRNA networks has been difficult to establish, due to the aforementioned limitations of the computational target prediction algorithms. Therefore, we applied our mRNA:miRNA targeting relationship predictions to build candidate ceRNA networks in the regenerating liver, selecting those mRNAs that had a single, strong miRNA footprint. Several examples of predicted miRNA ceRNA networks identified by our analysis are shown in Figure 5. Two members of the mir-140-3p ceRNA network, Spred2 and Vamp7, exert opposite effects on MAPK signaling pathways via inhibition of Raf phosphorylation [47] and endocytosis of the EGF receptor, respectively [48]. Thus, competition by these two mRNAs for mir-140-3p could balance positive and negative regulation of MAPK signaling in the regenerating liver, providing an example for the potential functional relevance of ceRNA networks.

A significant challenge in the field of microRNA biology is to understand physiologic miRNA:mRNA regulatory relationships in highly complex in vivo systems. We have demonstrated here that the HITS-CLIP assay can be applied to a complex model of growth and proliferation. We demonstrate that dynamic changes in miRNA recruitment to the RISC in most cases did not correlate with overall expression patterns in the regenerating liver, indicating that assessment of expression levels alone does not reflect miRNA activity. This lack of correspondence between RISC recruitment and miRNA expression is likely to be influenced by multiple factors including contributions of RNA binding proteins [49] ratio of target mRNA:targeting miRNA, flanking sequence homology [50] and changes in miRNA subcellular localization [51]. Using miRNA sensory and decoy libraries, Mullokandov and colleagues recently showed that only highly abundant miRNAs showed significant target repression, and in some instances were also less active based on a high mRNA:miRNA ratio or nuclear localization [51]. These findings suggest that the small changes in miRNA expression that occur in the regenerating liver may not contribute significantly to target mRNA regulation. We also show that inclusion of HITS-CLIP data enhances miRNA target predictions in hepatocytes by more than 20-fold over computational approaches alone. These findings provide a framework for understanding the relationship of miRNA regulation to changes in gene expression and will greatly facilitate the determination of candidate ceRNA networks.