MicroRNAs shape circadian hepatic gene expression on a transcriptome-wide scale

To comprehensively uncover miRNA-regulated gene expression in the liver, we generated a mouse model in which miRNA biogenesis could be inactivated in hepatocytes, which constitute around 80% of liver mass (Weibel et al., 1969). We used mice carrying conditional knockout alleles for the Dicer1 gene (termed Dicer^flox in the following), which encodes the ribonuclease that converts pre-miRNAs to mature miRNAs (Harfe et al., 2005), and a Cre-ER^T2 recombinase expressed as an internal ribosome entry site (IRES)-controlled transgene in the 3′ UTR of the endogenous Albumin locus (Alb^Cre-ERT2 in the following) (Schuler et al., 2004). The Albumin locus conferred hepatocyte-specificity and the Cre-ER^T2 fusion allowed for induction of knockout in adult animals using tamoxifen.

Previous publications have reported constitutive hepatocyte-specific knockouts of Dicer (Hand et al., 2009; Sekine et al., 2009a, 2009b) and to our knowledge the alleles we employed have not been used in combination before. We thus first confirmed their suitability for efficient miRNA depletion in the livers of adult animals. To this end, we analysed the kinetics of recombination at the Dicer locus at different timepoints after intraperitoneal tamoxifen (tx) injection. Recombination was undetectable before but highly efficient from 2 to 3 weeks after tx treatment (Figure 1—figure supplement 1A). We assessed the consequences on miRNA expression by northern blot analysis probing for the highly abundant, hepatocyte-specific miR-122, which has been previously noted for its long half-life (Gatfield et al., 2009) and thus served as a proxy for overall miRNA depletion from hepatocytes. One month after tx treatment, miR-122 was virtually undetectable (Figure 1—figure supplement 1B,C). Moreover, loss of mature miR-122 was accompanied by increased levels of its precursor, pre-miR-122 (Figure 1—figure supplement 1B,C), as expected (Harfe et al., 2005). Interestingly, we observed that at later timepoints beyond 2 months after tx treatment, miR-122 expression began to recover (Figure 1—figure supplement 1D,E). The most likely explanation for this observation was that after miRNA loss the liver eventually renewed its hepatocytes from a cell population that had escaped recombination. In summary, we concluded that our experimental system appeared suitable for the efficient knockout of Dicer and, at least judging by miR-122, for the depletion of miRNAs. Approximately 1 month after tx treatment appeared to represent a suitable timepoint for experiments aiming at a comparison of knockout with control animals.

We induced the knockout in more than 40 mice (male, aged 3–6 months), entrained them to light–dark cycles for 1 month (LD 12:12; ad libitum feeding) and sacrificed them at 4-hr intervals around-the-clock (Figure 1A). A littermate control group was treated identically. We combined both wild-type and heterozygous males for controls because a single functional Dicer allele was sufficient to ensure miRNA processing to wild-type levels (Figure 1—figure supplement 2A), as expected from previous studies (Harfe et al., 2005; Kanellopoulou et al., 2005; Murchison et al., 2005; Chen et al., 2008; Frezzetti et al., 2011). Pools of liver total RNA (3–4 animals) were assembled into two independent full time series around-the-clock (ZT0-20) for both knockouts and controls. Northern blot analysis confirmed loss of miR-122 in the knockouts across all timepoints (Figure 1B,C). Microarray (Figure 1—figure supplement 2B) and northern blot analyses (Figure 1—figure supplement 2C) showed that beyond miR-122, miRNA levels were globally decreased. Other hepatocyte-enriched miRNAs (e.g., miR-148a, miR-192 and miR-194 [Sun et al., 2004; Landgraf et al., 2007; Tang et al., 2007; Farid et al., 2012]) were thus virtually undetectable in knockouts. Probably owing to their expression also in non-hepatocyte cells, known ubiquitous miRNAs (e.g., miR-21 [Landgraf et al., 2007], miR-26a, or let-7 family members [Lagos-Quintana et al., 2003]) decreased less strongly, and miRNAs with a particularly strong expression in non-hepatocytes (e.g., miR-126 in hepatic stellate cells [Guo et al., 2013]) were almost unchanged, as expected (Figure 1—figure supplement 2B,C). In all cases, reduced miRNA levels correlated with increased pre-miRNA abundance (Figure 1B,D, Figure 1—figure supplement 1C). Moreover, even in the knockout, pre-miR-122 still showed rhythmic accumulation with a trough at ZT12 (Figure 1B,D), which was a consequence of clock-driven rhythmic transcription at the mir122 gene locus (Gatfield et al., 2009; Menet et al., 2012). Other clock-controlled genes such as Dbp (Ripperger et al., 2000) still oscillated as well (Figure 1E). We thus concluded that cyclic gene expression per se was still functional in Dicer knockouts and that the RNA time series collected from the two genotypes was suitable to uncover the effects of miRNA loss on the rhythmic expression of core clock and clock output genes in hepatocytes.

MicroRNAs inhibit the translation and promote the decay of their mRNA targets. How these events are elicited, and in which order they occur, is still subject of ongoing debate. However, given that in most cases mRNA degradation appears to be an endpoint in miRNA activity, changes in mRNA levels are the most commonly used readout for regulation by miRNAs (reviewed in Fabian and Sonenberg, 2012).

To analyse how rhythmic RNA accumulation was affected by the loss of miRNAs, we employed RNA high-throughput sequencing (RNA-seq). We chose a protocol for library generation that was based on the random priming of rRNA-depleted (‘Ribo-Zero’) total RNA, which permitted the simultaneous quantification of mRNAs (exon-mapping reads) and pre-mRNAs (intron-mapping reads). Given that splicing occurs co-transcriptionally and is relatively fast, the intron-containing pre-mRNAs are short-lived and their abundance can serve as a proxy for gene transcription rates (Koike et al., 2012). This strategy thus identifies post-transcriptional (mRNA) gene expression changes that occur in the Dicer knockouts and likely represent bona fide miRNA targets, as opposed to transcriptional changes (mRNA and pre-mRNA affected), which likely represent secondary effects of miRNA loss. Conceivably, any transcription factor that is a direct miRNA target would engender numerous such secondary effects.

We sequenced knockout and control time series with comparably high coverage. Reads mapped to different RNA classes, as expected (Figure 2A, Figure 2—figure supplement 1, Figure 2—source data 1). We did not detect gross distortions in RNA populations that would compromise the comparability of knockout and control samples, although some differences reached statistical significance (Figure 2—figure supplement 1, Figure 2—source data 1). We then examined transcripts from protein-coding loci, which contributed around half of all obtained sequences (51.6 ± 1.0% in knockout, 47.9 ± 1.3% in control; Figure 2A,B, Figure 2—figure supplement 1, Figure 2—source data 1). Reads mapping to exonic sequences were categorised as ‘mRNA’ and those mapping to introns or intron–exon boundaries as ‘pre-mRNA’ (‘Materials and methods’). Pre-mRNAs are of generally low abundance, but genome-wide there is an estimated 20-fold more transcribed intronic than exonic sequence space (Sakharkar et al., 2005). In accordance, pre-mRNAs were readily detectable transcriptome-wide, representing around a quarter of protein-coding reads (25.5 ± 4.1% in knockout, 21.5 ± 3.1% in control; Figure 2A,B, Figure 2—source data 1). Finally, the RNA-seq data also provided an additional confirmation of the high efficiency of the Dicer knockout (Figure 2—figure supplement 2).

In the complete data set, we detected 13,408 individual protein-coding loci (Figure 2C). The intersecting set (mRNA and pre-mRNA detected in knockout and control) contained 11,990 genes (89.4%); another 361 transcripts (2.7%) common to both genotypes corresponded to intronless genes. To assess the changes engendered by miRNA loss, we first compared mRNA and pre-mRNA levels separately between the genotypes (averaging over timepoints). It was apparent that Dicer deficiency caused large changes in mRNA and pre-mRNA abundance (Figure 2D,E). Applying a 1.5-fold cut-off on differential expression, we found that 30.3% of mRNAs and 15.5% of pre-mRNAs were significantly changed in the Dicer knockout. These changes were skewed towards increased levels for mRNAs (19.8% increased, 10.5% decreased) but slightly less so for pre-mRNAs (9.8% increased, 5.7% decreased). In summary, these findings confirmed our expectation that global miRNA loss would lead to the widespread de-repression of targets that affected a substantial part of the transcriptome. Inevitably, these primary effects would provoke numerous secondary responses involving altered transcription. Striking transcriptional effects were detected for several imprinted loci and for a group of known development-specific and sexually dimorphic genes (see Figure 2—figure supplement 3 for examples), as previously reported from Dicer knockouts (Hand et al., 2009; Sekine et al., 2009a). We excluded the group of transcriptionally most altered genes in subsequent circadian analyses, as vastly different transcription rates would likely confound comparative analyses of rhythmic properties, in particular regarding the precise post-transcriptional contributions (marked in red in Figure 2D,E).

We next calculated transcriptome-wide mRNA/pre-mRNA ratios, which are a measure of mRNA stability and thus an indicator of post-transcriptional regulation (Zeisel et al., 2011). The distribution of mRNA/pre-mRNA ratios was globally shifted to higher values in the knockout (Figure 2F) and we thus assessed whether an increased mRNA/pre-mRNA ratio in the Dicer knockout could be predictive of miRNA activity. Using experimentally validated targets (Hsu et al., 2011) of miR-122 as a positive control and of miR-124 as a negative control (this miRNA is not expressed in liver, [Landgraf et al., 2007]), we indeed observed that the majority of miR-122 targets showed increased mRNA/pre-mRNA ratios in Dicer knockouts, whereas most miR-124 targets did not (Figure 2G). Ribosomal protein mRNAs, whose 3′ UTRs are particularly short (Caldarola et al., 2009) and thus probably not regulated by miRNAs, did not show increased mRNA/pre-mRNA ratios either (Figure 2G). We thus concluded that incorporating mRNA/pre-mRNA ratios in our analyses could enrich for likely direct miRNA targets in the dataset.

We next assessed the expression of transcripts encoding core clock components and clock regulators (listed in Figure 3—source data 1). With regard to mRNA/pre-mRNA ratios, the transcripts of the Period gene family (Per1, Per2, Per3) showed a moderate (average 1.6-fold), but significant increase in Dicer knockouts (Figure 3A). The time-resolved RNA-seq data (Figure 3B) confirmed that Per1 and Per2 mRNAs accumulated to higher levels throughout the day despite similar rhythmic pre-mRNA abundance. For Per3, the higher mRNA/pre-mRNA ratio mainly resulted from decreased transcription (Figure 3—figure supplement 1). Moreover, we noticed that Cry2 showed post-transcriptional up-regulation in the Dicer knockout and was thus a potential miRNA target, although the case was less clear than for Per1 or Per2 since Cry2 transcription was increased at certain timepoints as well. Beyond these genes, the expression of most other core clock components, including Cry1, Clock, Bmal1 (Figure 3B) and others (Figure 3—figure supplement 1) only showed minor changes, if at all. Consistent with the mRNA profiles, western blot analysis confirmed an estimated twofold increase in the levels of PER2 protein and 1.5-fold higher CRY2 (Figure 3C). Surprisingly, PER1 levels were barely increased (1.1-fold higher). A number of mechanisms could account for the observed uncoupling of protein from mRNA levels, including secondary effects of miRNA loss acting on Per1 mRNA translation or the protein degradation machinery. While we cannot pinpoint the exact mechanism that is operative in counter-regulating Per1 mRNA increase, it is interesting to note that the RNA binding protein Syncrip (also known as hnrnp q) is a positive regulator of Per1 translation (Lee et al., 2012) and showed decreased expression in Dicer knockouts (Figure 3—figure supplement 1). Importantly, however, it was clear that core clock functionality was overall unimpaired, as also exemplified by the virtually unchanged expression profiles of genes that are directly regulated by CLOCK:BMAL1, PERs and CRYs, such as Dbp (Figure 3—figure supplement 1, Figure 1E).

In our knockout model, miRNA loss occurred in the hepatocytes of a liver that was under constant entrainment by a genetically intact SCN clock. It was thus conceivable that a more severe phenotype was masked and would only become apparent under non-entrained, free-running conditions. It has indeed recently been reported that the period length of free-running rhythms in Dicer–deficient MEFs is dramatically shortened by ≈2 hr (Chen et al., 2013). Moreover, the authors identified the lack of miRNA-mediated repression of Per1 and Per2 as the underlying mechanism. We thus measured free-running rhythms in liver explants from Dicer knockout and control mice carrying the mPer2^Luc reporter allele (Yoo et al., 2004). Interestingly, our experiments showed a trend towards period lengthening in knockout liver explants with a period that was on average 41 min longer than in controls. The effect, however, did not reach statistical significance (p=0.197) (Figure 3D,E). Overall, we concluded that Per1, Per2 and Cry2 were likely miRNA targets in liver, but that the hepatocyte clock was relatively resilient to miRNA loss. The difference in phenotypes of liver and MEF Dicer knockouts could point to tissue-specific miRNA functions in clock regulation (‘Discussion’).

We next explored how circadian gene expression was affected beyond the core clock. We first wished to identify transcripts whose rhythms were post-transcriptionally driven by miRNAs and assess whether this mechanism could underlie the discrepancies between cyclic transcription and cyclic mRNA accumulation that have been reported (Koike et al., 2012; Le Martelot et al., 2012; Menet et al., 2012). To identify oscillating gene expression transcriptome-wide and as comprehensively as possible, we adapted previously published methods for rhythmicity detection that we used at relatively low stringency (‘Materials and methods’). Especially for pre-mRNAs, which are naturally noisier due to their low abundance, we wanted to avoid applying overly stringent criteria, as this could lead to the under-detection of transcriptionally and the overestimation of post-transcriptionally generated rhythms. Indeed, a previous study had noted that many genes that had failed to pass the imposed rhythmicity criteria on the transcriptional level visibly still exhibited transcriptional patterns that resembled those of the mRNA, albeit with much higher variability (Menet et al., 2012).

In control and Dicer knockout mice, 1630 and 1902 mRNAs, respectively, cycled with at least 1.5-fold peak-to-trough amplitudes (≈13% and 14%, respectively, of detected mRNAs; around half were shared between both groups; Figure 4A, Figure 4—source data 1), consistent with previous estimates that up to 15% of mRNAs are rhythmic in wild-type mice (Lowrey and Takahashi, 2004; Vollmers et al., 2009). On the pre-mRNA level, around 17% of transcripts were rhythmic (2083 in control and 2101 in knockouts; >50% shared). In both knockouts and controls, >60% of oscillating mRNAs were also rhythmic at the pre-mRNA level. These values are within the range of what has been reported previously from wild-type mice (Koike et al., 2012; Le Martelot et al., 2012; Menet et al., 2012). Interestingly, a sizeable number of transcripts were rhythmic specifically in either controls or Dicer knockouts (Figure 4A, Figure 4—figure supplement 1). We were particularly intrigued by the relatively high number of transcripts that appeared to become rhythmic in the Dicer knockouts (Figure 4A, sectors B–D) and investigated whether this was caused as a side effect of low stringency rhythm detection. However, the overlap of rhythmic transcripts between genotypes did not increase with higher thresholds despite reducing the overall number of detected rhythmic events (Figure 4—figure supplement 2). Further analyses revealed that Dicer KO-specific rhythmic transcripts were strongly enriched for certain Gene Ontology (GO) terms, of which cell cycle-specific genes (especially those involved in DNA replication, such as the members of the minichromosome maintenance complex, MCM) were particularly noteworthy (Figure 4—source data 2, see Figure 4—figure supplement 3 for examples). These findings are in line with increased liver regeneration after miRNA loss (Figure 1—figure supplement 1D,E) and the circadian gating of the cell cycle in the regenerating liver (Matsuo et al., 2003).

To explore the relationship between transcript rhythmicity and miRNA regulation, we first consulted the mRNA/pre-mRNA ratios as an indicator of post-transcriptional regulation and of mRNA stability changes that could be indicative of direct miRNA activity (Figure 4B,C). For 30% of rhythmic transcripts in control mice (483 out of 1630), mRNA/pre-mRNA ratios were significantly increased in Dicer knockouts (Figure 4—figure supplement 1A, sectors I–P). Globally, rhythmic transcripts showed a slightly stronger shift to higher mRNA/pre-mRNA ratios than average (Figure 4C; p=0.001; Welch's two-sample t test). We concluded that several hundred circadian mRNAs were potentially miRNA-regulated but that there was no striking enrichment or depletion for miRNA regulation among the circadian transcriptome.

To identify candidates for post-transcriptionally driven rhythmicity involving miRNAs, we investigated the group of 290 transcripts whose cyclic accumulation occurred exclusively at the mRNA level in control animals (sector I in Figure 4A, see Figure 4—source data 1). A heatmap representation confirmed pronounced rhythmicity mainly at the control mRNA level, as expected, but also suggested that a considerable number of transcripts possibly still possessed underlying, yet noisier rhythms at the level of control pre-mRNA and/or Dicer knockout mRNA (Figure 4D). We thus examined all profiles individually, also taking into account mRNA/pre-mRNA ratio changes occurring in the Dicer knockouts (in total 93 of the 290 transcripts in sector I thus had increased mRNA/pre-mRNA ratios in knockouts, see Figure 4—figure supplement 1). The visual inspection confirmed that only a minority of transcripts showed the features expected for true miRNA-driven post-transcriptional rhythmicity and that, instead, many of the non-rhythmic assignments were probably a result of noisier data (Supplementary file 1 in Du et al., 2014 and Figure 4—source data 1 for whole transcriptome gene expression plots and rhythmicity parameters). Visual selection resulted in no more than 20 transcripts for which we would confidently postulate that their rhythmic accumulation could indeed be miRNA-driven. As shown for Ddx17, Slc1a5, Stx2, Uba6, Zfp697 (Figure 4E–I) and others (Figure 4—figure supplement 4), these transcripts were all characterised by relatively low amplitude rhythms. As a validation with an independent technique, we confirmed the circadian expression detected by RNA-seq for two of the genes, Stx2 and Uba6, also by quantitative real-time PCR (qPCR) (Figure 4—figure supplement 5A). In summary, these findings suggested that miRNAs are only of marginal importance as drivers of circadian rhythmicity. Although we were able to identify several transcripts whose cyclic accumulation was Dicer-dependent and occurred post-transcriptionally, miRNA-mediated mechanisms can probably not account for the previously observed large discrepancy between mRNA and transcription rhythms (Koike et al., 2012; Le Martelot et al., 2012; Menet et al., 2012).

Hundreds of rhythmic mRNAs displayed increased stability (as judged from mRNA/pre-mRNA ratios) in Dicer knockouts (Figure 4B,C) and many of these transcripts likely represented direct miRNA targets. Circadian transcription paired with miRNA-mediated post-transcriptional control is thus probably a widespread phenomenon. We therefore investigated how such dual regulation added up at the level of mRNA cycling. Kinetic models of how mRNA stability influences rhythmic mRNA accumulation show that for a cyclically transcribed gene, the more stable the transcript, the later the phase and the lower the amplitude of its cycling (Le Martelot et al., 2012), as schematically shown in Figure 5A.

We examined the circadian properties of transcripts that were rhythmic across all conditions, that is in knockouts and controls at the pre-mRNA (>1.5-fold amplitude) and mRNA levels (without amplitude cut-off). Moreover, we grouped the transcripts according to whether their mRNA/pre-mRNA ratios were significantly increased in the Dicer knockout (i.e., transcripts more likely to represent direct miRNA targets; N = 167) or not (i.e., transcripts less likely to represent direct miRNA targets; N = 505). A heatmap representation of the transcripts considered is shown in Figure 5B (see also Figure 5—source data 1 for gene list).

We first determined how miRNA loss affected the phases of rhythmic mRNA accumulation (Figure 5C). We thus calculated the time lag between the transcripts' pre-mRNA and mRNA peaks for the Dicer knockouts, and subtracted the equivalent values calculated from the control mice. Positive and negative ‘phase delay differences’ thus signified that in the Dicer knockout, mRNA accumulation peaked later and earlier, respectively, than in controls (an additional analysis of absolute pre-mRNA and mRNA phases can be found in Figure 5—figure supplement 1A,B). Phase delay differences were plotted against the mRNA/pre-mRNA ratio change (Figure 5C; the 167 transcripts with a significantly higher ratio in the knockout are marked in red). This analysis revealed a shift to later phases of mRNA accumulation in the Dicer knockout, which correlated with increased mRNA/pre-mRNA ratios (blue line in Figure 5C; Pearson's r(670) = 0.17; p=9.24e−06). In the absence of miRNAs, mRNAs with significantly increased mRNA/pre-mRNA ratios thus peaked on average 48 min later (vertical red line in Figure 5C), whereas the average for all other rhythmic transcripts (vertical black line) lay at 5 min. Both groups of transcripts, that is those that showed higher mRNA/pre-mRNA ratios and those that did not, were significantly different (p-value 1.72e−07; Wilcoxon rank sum test; cumulative density plot at the top of Figure 5C). Although the globally detectable effect of miRNA loss on the phases of mRNA rhythms was relatively modest (less than 1 hr), individual transcripts were affected more dramatically (Figure 5E). For genes such as Aconitase 2 (Aco2), F-box protein 21 (Fbxo21), Nicotinamide phosphoribosyltransferase (Nampt) or Nerve growth factor receptor (Ngfr), the phase-shifts amounted to several hours (Figure 5E). We also confirmed the observed effects for selected genes by qPCR (Figure 4—figure supplement 5B).

Using a similar approach, we next assessed how the loss of miRNAs influenced the amplitudes of mRNA rhythms. For each transcript, we calculated the ratio of the mRNA amplitude to the pre-mRNA amplitude in the Dicer knockout and divided this value by the corresponding ratio in control mice. We thus found that miRNA loss led to globally shallower amplitudes in the Dicer knockout (Figure 5D), although the mean effect was small and amounted only to a 7% decrease in fold-amplitude (see difference between red and black vertical lines in Figure 5D; p-value for difference between transcripts with and without increased mRNA/pre-mRNA ratio: 9.39e−03; Wilcoxon rank sum test). Nevertheless, individual examples clearly suggested that miRNAs assume important functions in ensuring that rhythmically transcribed genes give rise to rhythmic mRNAs that oscillate with the desired amplitudes and magnitudes, as shown in Figure 5F (qPCR validations in Figure 4—figure supplement 5C,D) for the expression of DNA methyltransferases 3A (Dnmt3a), 3B (Dnmt3b) and Neuroepithelial cell transforming gene (Net1), as well as for several other genes in Figure 5E. In some cases, such as Lysosomal acid lipase A (Lipa), rhythmicity generated at the pre-mRNA level was completely lost at the mRNA level when miRNAs were absent (Figure 5G). Here, miRNA-mediated decay could thus represent the mechanism ensuring that after cyclic transcription, mRNAs are sufficiently unstable to show rhythmic accumulation at all. Interestingly, the protein encoded by Lipa, lysosomal acid lipase (LAL), is the key enzyme hydrolysing cholesteryl esters and triglycerides stored in lysosomes after LDL receptor-mediated endocytosis (Fouchier and Defesche, 2013) and has been reported to be rhythmic in liver, but non-rhythmic in other organs (Tanaka et al., 1985). Conceivably, miRNAs could be involved in rendering Lipa rhythms tissue-specific.

Altogether the effects of miRNA loss on amplitudes were less uniform than those on phases (compare Figure 5C,D). Several reasons could account for this difference. First, we noticed that amplitude estimations appeared technically more error-prone than phase estimations. Both parameters were read from cosine curves fitted to the data, but the peak-trough symmetry imposed by the cosine function frequently resulted in the underestimation of amplitudes; this was especially the case for high amplitude rhythms with pronounced ‘spiky’ appearance (i.e., with rapid rising and declining phases; see cosine fits of Nr1d1/Rev-erbα and Dbp in Figure 3—figure supplement 1). Second, miRNAs likely only represent one of several mechanisms operative in amplitude modulation; direct consequences of miRNA loss on target mRNA amplitudes might thus be partially masked by secondary effects occurring in the Dicer knockouts. Third, the model of how phases and amplitudes should ideally correlate (Figure 5A) is almost certainly an oversimplification. Indeed, for transcripts such as Cpeb1 (Figure 5H), later phases were even accompanied by higher amplitudes. Globally, there was indeed no significant correlation between phase delays and amplitude decreases (Figure 5—figure supplement 1E). An uncoupling of amplitude and phase effects could occur, for example, through miRNA activity that is not constant over the day but confined to specific timepoints (‘Discussion’). Overall, we concluded that the phase delays seemed to be the dominant, consistent consequence of miRNA loss in liver. Nevertheless, both the amplitude decrease and the phase delay were individually highly significant, as also shown by permutation tests in which we randomly reshuffled either the timepoints of the original data or the Dicer KO/control assignments on a genewise basis, followed by analyses of rhythm parameters. Among >10,000 datasets with permuted timepoints, the likelihood of finding similar or stronger delays in phase difference (in both size and significance) or correlations (between mRNA/pre-mRNA ratios and phase difference) than those in the original data were very low (Figure 5—figure supplement 2A). The observed shift in phase differences was therefore not an intrinsic data set property. By contrast, because the expression level differences (Dicer KO vs control) were untouched by the reshuffling of timepoints, a sizeable number of permutations still showed amplitude effects comparable to the original data (Figure 5—figure supplement 2B). However, when reshuffling occurred at the level of the genotypes, the permuted data no longer retained any trend reminiscent of the original data, neither for phases nor for amplitudes (Figure 5—figure supplement 2C,D). We concluded that the correlations uncovered between miRNA loss and rhythmic gene expression parameters were not likely to have arisen by chance, for example as a result of a specific predisposition of the data structure or distribution that could have occurred in our particular time series and the downstream analyses.

The correct timing and extent of rhythmicity is important for the daily execution of clock-regulated physiological functions. It was thus likely that the observed phase and amplitude changes had an impact on liver functions. To identify pathways that were particularly affected by miRNA loss, we investigated whether the post-transcriptionally up-regulated circadian transcripts were associated with specific GO terms (Figure 6A). Using the same data set of circadian transcripts as before (Figure 5), several over-represented GO terms could be identified with statistical significance (FDR-adjusted p-value<0.05) in the group of rhythmic transcripts with increased mRNA/pre-mRNA ratios in Dicer knockouts, but not among the remaining rhythmic transcripts. However, the enrichment was relatively low, indicating that the function of miRNAs in regulating circadian gene expression was broad rather than confined to specific rhythmic functions. The individual GO terms affected three major fields: cell adhesion, apoptosis/development, and (lipid) metabolism. In particular, the latter caught our attention because serum analyses of Dicer knockout animals indicated metabolic phenotypes as well, especially for cholesterol, triglyceride, and glucose levels, which were all reduced (Figure 6—figure supplement 1). At least in part, these effects were likely mediated by the loss of miR-122, whose inactivation has previously been shown to impact lipid metabolism (Krutzfeldt et al., 2005; Esau et al., 2006; Gatfield et al., 2009; Wen and Friedman, 2012).

The breadth of miRNA loss (including potential secondary effects) poses an obvious limit to the suitability of Dicer knockouts for further functional studies. Ideally, the targeted inactivation of single miRNAs that fulfil prospective functions in circadian regulation would allow for more straightforward investigations. To identify such miRNAs, we analysed whether the circadian transcripts that became more stable in Dicer knockouts were enriched for predicted miRNA binding sites (Lewis et al., 2005). Interestingly, this analysis revealed nine miRNA seed families with more circadian targets than expected by chance (Figure 6B) specifically in the data set with higher mRNA/pre-mRNA ratios. Of these, miR-25, miR-92, miR-122, miR-22, and miR-29 are expressed at readily detectable levels in liver, as shown previously (Landgraf et al., 2007; Gatfield et al., 2009; Vollmers et al., 2009) and confirmed by northern blot (Figure 1—figure supplement 2) and small RNA-seq (data not shown). Given that miR-122 already has a history as a circadian output modulator (Gatfield et al., 2009), it is tempting to speculate that the other miRNAs may represent specialised regulators of the rhythmic transcriptome as well. Interestingly, miR-29 is also one of the miRNAs that was reported to regulate Per1 in MEFs (Chen et al., 2013). In the future, it would thus be exciting to explore the role of these miRNAs in targeted experiment in the context of clock-regulated biological pathways.

The combination of increased mRNA/pre-mRNA ratios and miRNA target predictions is suggestive of bona fide miRNA-mediated regulation. Nevertheless, before embarking on follow-up experiments for individual transcripts, validation experiments to map miRNA binding sites in the targets' 3′ UTRs and to identify the responsible miRNAs will be necessary. To this end, we designed a luciferase reporter assay based on the lentiviral transduction of Dicer knockout and control primary hepatocytes (Figure 6C). We chose the core clock transcripts to demonstrate the utility of the assay, as we were particularly intrigued by the difference in phenotypes between the Dicer knockout in liver (i.e., mild period lengthening, Figure 3E) and the previously reported strong period shortening in MEFs (≈2 hr), which the authors attributed to miRNA-regulation of Per1 (by miR-24 and miR-29) and Per2 (by miR-24 and miR-30) (Chen et al., 2013). We tested whether the firefly luciferase reporter gene carrying the full-length Per1, Per2 or Cry2 3′ UTRs (or fragments thereof) was subject to de-repression in Dicer knockout hepatocytes, as compared to control hepatocytes. Consistent with the RNA-seq and western blot analyses (Figure 3B,C), this assay indeed identified Per2 as the main miRNA-regulated core clock component (Figure 6D). Moreover, truncated versions of the 3′ UTR that contained the miR-24 and miR-30 sites identified by Chen et al. recapitulated some of the miRNA regulation observed for the full-length Per2 3′ UTR sequence. However, additional regulation was also conferred by a portion of the 3′ UTR containing a predicted target site for the miR-25/92 seed family (fragment 1739–2129, Figure 6D) and full regulation may thus require the synergistic activity of several sites. Per1 and Cry2 full-length 3′ UTRs showed 1.2-fold and 1.8-fold up-regulation, respectively (Figure 6D). This was in the same range as what we had observed for the endogenous proteins by western blot (Figure 3C). In particular, the lack of strong Per1 regulation in liver was a striking difference to the MEF data and could be responsible for phenotypic differences, as Chen et al. have proposed that the faster translation and higher accumulation of PERs (approximately twofold higher PER1 and PER2 peak levels in Dicer knockout) shortens the time delay within the oscillator's main feedback loop and speeds up the clock. It is possible that tissue-specific miRNA expression contributes to the observed differences; in line with a previous study (Vollmers et al., 2012), we have thus observed that miR-30 is among the most abundant miRNA species in liver (within the top 10 of miRNAs detected by small RNA-seq), but that miR-29 and miR-24 are expressed at lower levels (not in the top 50) (data not shown), which could thus partially explain the comparably weak effect of the Dicer knockout on PER1 in liver. In summary, in addition to shedding light on miRNA regulation of core clock transcripts, we concluded that the reporter assay in hepatocytes represented a useful tool for future studies aimed at verifying direct miRNA targets and at unravelling the contribution of specific miRNAs.

Animal studies were conducted in accordance with the regulations of the veterinary office of the Canton of Vaud (authorization VD2376). All alleles used in the study have been published before that is, Dicer^flox (Harfe et al., 2005), Alb^Cre-ERT2 (Schuler et al., 2004), and mPer2^Luc (Yoo et al., 2004), and were kindly provided by the Tabin, Metzger, and Takahashi labs, respectively. The genetic background of the animals used in this study was mixed. Male littermate knockout (Dicer^flox/flox;Alb^Cre-ERT2) and control (Dicer^flox/+ or Dicer^+/+;Alb^Cre-ERT2) animals aged 3–6 months received tamoxifen treatment by intraperitoneal injections over 5 days (total 2 mg) essentially as described (Schuler et al., 2004). We combined both wild-type and heterozygous male littermates for the controls because microarray analyses indicated that a single functional Dicer allele was sufficient to ensure miRNA processing to wild-type levels (Figure 1—figure supplement 2A), as expected from previous studies (Harfe et al., 2005; Kanellopoulou et al., 2005; Murchison et al., 2005; Chen et al., 2008; Frezzetti et al., 2011). Haploinsufficiency has indeed so far only been reported for tumour suppressor functions of Dicer in certain non-liver tumorigenesis models (Kumar et al., 2009; Arrate et al., 2010; Lambertz et al., 2010; Nittner et al., 2012; Yoshikawa et al., 2013).

After the last injection, mice were entrained to a 12-hr light:12-hr dark photoperiod with free access to food and water for 1 month, sacrificed at the respective Zeitgeber times (ZT 0, 4, 8, 12, 16, 20), and the liver was snap-frozen in liquid nitrogen. Oligos used to genotype Dicer knockout efficiency from genomic DNA (Figure 1—figure supplement 1A) were DicerF1, DicerR1, and DicerDel (see Supplementary file 2 in Du et al., 2014 for sequences).

Total RNA was prepared essentially as described previously (Gatfield et al., 2009). The efficiency of the Dicer knockout was confirmed for each liver sample by northern blot probing for miR-122. The protocol for miRNA northern blots has been described (Gatfield et al., 2009). Total RNA pools for the two time series in knockouts and controls were assembled using identical amounts of RNA from 3 to 4 animals per pool. All pools were DNase-digested (RQ1 DNase, Promega, Madison, WI) in order to eliminate potential contamination from genomic DNA that would have rendered pre-mRNA quantifications impossible. The absence of genomic DNA contamination was confirmed in all pools by the comparison of RT− and RT+ reactions in quantitative real-time PCR (qPCR) using genomic probes. For qPCR analysis, cDNA was synthesised from 6 μg of DNase-treated total RNA using random hexamers and SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA) according to the supplier's instructions. cDNA was PCR-amplified using FastStart Universal SYBR Green Master (Roche, Basel, Switzerland). Mean levels were calculated from triplicate PCR assays for each sample and normalised to those obtained for the control transcripts GusB and Eef1a1 (Figure 1E) and Csnk1a1, Ctsd, Nudt4, Smg1 and Trip12 (Figure 4—figure supplement 5). All oligonucleotide sequences are listed in Supplementary file 2 (Du et al., 2014).

Total RNAs (100 ng; same samples as for the RNA-seq series) were hybridised to Mouse miRNA Microarray, Release 18.0, 8 × 60K (Agilent Technologies, Santa Clara, CA) using miRNA Complete Labelling and Hyb Kit (Agilent Technologies) according to the supplier's instructions. An invariant normalisation method using five internal spikes (hur_1, hur_2, hur_4, hur_6, mr_1) was used to normalise Dicer KO vs control samples in Figure 1—figure supplement 2B. Quantile normalisation was used for heterozygote knockout vs wild-type samples in Figure 1—figure supplement 2A.

The dual luciferase reporter cassette in which the bidirectional phosphoglycerate kinase 1 (PGK1) promoter (Amendola et al., 2005) drives the expression of firefly luciferase (FL, used for 3′ UTR cloning) and renilla luciferase (RL, for normalisation) was designed in silico and purchased as a synthetic clone (GenScript, Piscataway, NJ) in pUC57 plasmid. After excision from the plasmid via flanking SalI sites, the cassette was used to replace the XhoI cassette in lentiviral plasmid pWPT-GFP (Addgene, Cambridge, MA), resulting in plasmid prLV1 that was used to clone 3′ UTRs of interest via XhoI/NotI restriction sites downstream of the firefly luciferase coding sequence. 3′ UTR sequences were amplified by PCR from liver cDNA or genomic DNA with specific oligonucleotides that carry NotI sites in the forward (F) and XhoI sites in the reverse (R) oligos, as listed in Supplementary file 2 (Du et al., 2014). The identity of the cloned UTRs was verified by sequencing.

From the prLV1 vectors with cloned UTRs, lentiviral particles were produced in 293T cells using envelope vector pMD2.G and packaging plasmid psPAX2 as previously described (Salmon and Trono, 2007). Viral supernatant was spun 2 hr at 25,000 rpm, 4°C using Optima L-90K Ultracentrifuge (SW32Ti rotor; Beckman, Brea, CA), then viral particles were resuspended with primary hepatocyte medium.

Mice were anesthetised and livers were perfused through the inferior vena cava with 50 ml of washing buffer (137 mM NaCl, 2.7 mM KCl, 0.5 mM Na₂HPO₄, 10 mM HEPES, pH 7.65) supplemented with 0.5 mM EDTA, then 50 ml of digestion buffer (washing buffer supplemented with 7 mM CaCl₂, 0.4 mg/ml collagenase [C5138; Sigma-Aldrich, St. Louis, MO]) at flow rate of 5 ml/min and at 37°C. Isolated cells were filtered with a cell strainer (100 μm, BD Falcon, Franklin Lakes, NJ) and washed with 40 ml of primary hepatocyte medium (Medium 199, GlutaMAX Supplement, supplemented with 1% penicillin/streptomycin/glutamine (PSG), 0.1% BSA, 10% FCS [all Gibco/Life Technologies, Carlsbad, CA]). The cells were resuspended with 10 ml of primary hepatocyte medium, counted, and plated in 12-well plates (coated with 0.2% gelatin) at a density of 2 × 10⁵ cells/well. Medium was changed 4 hr later and viral supernatant was added. The cells were maintained at 37°C and 5% CO₂ until harvested. For each lentiviral construct, transductions were performed in triplicates. Each experiment was repeated at least twice with cell preparations from independent mice.

6 days after lentiviral transduction, cells were collected using 5x Passive Lysis Buffer (Promega) and luciferase activity was measured using the Dual-Glo Luciferase Assay System (Promega) according to the manufacturer's protocol. Firefly luciferase signals were normalised to Renilla luciferase, and for each 3′ UTR construct this signal was then normalised to that of lentivector-control plasmid (without cloned 3′ UTR, thus only containing generic vector 3′ UTR). Signals in control mice were set to 1.

Mice were anesthetised with isoflurane and sacrificed by decapitation. The liver was excised and immediately placed in cold Hank's balanced salt solution (Invitrogen). Then, the liver was sliced into small pieces and cultured separately on Millicell culture membranes (PICMORG50; Millipore, Billerica, MA) with 1.1 ml of HEPES-buffered phenol red-free DMEM (Gibco) supplemented with 2% B27 (Invitrogen), 1% PSG (Gibco), 4.2 mM Na₂HCO₃ and 0.1 mM luciferin. Cultures were maintained at 37°C and 5% CO₂ in a light-tight incubator, and bioluminescence was monitored continuously using LumiCycle 32 (Actimetrics, Wilmette, IL).

Proteins from mouse liver nuclei were prepared according to the NUN procedure (Lavery and Schibler, 1993). Each sample analysed by SDS-PAGE was a pool of protein extracts from 3 to 5 mice. SDS-PAGE and immunoblot analysis were performed according to standard protocols. Antibodies used were rabbit CRY1, CRY2, PER1, PER2, BMAL1, and CLOCK (kindly provided by S Brown and J Ripperger) and U2AF65 (Sigma). Western blots were quantified using ImageQuant TL 8.1 software (Life Sciences).

All blood chemistry was performed as described (Le Martelot et al., 2009).

DNase-treated total RNAs were subjected to rRNA depletion using Ribo-Zero Magnetic Kit (Human/Mouse/Rat, Epicentre, Madison, WI) following the supplier's instructions. The resulting rRNA-depleted RNA samples were then used to prepare random-primed cDNA libraries using the Illumina Tru-Seq RNA Sample Preparation Kit (Illumina, San Diego, CA) according to the manufacturer's recommendations. Multiplexed libraries were sequenced following the supplier's protocol (100 bp single-end reads) on the Illumina HiSeq 2500 at the Lausanne Genomics Technologies Facility.

Quality of the sequencing reads was initially assessed based on the quality values produced by the Illumina pipeline Casava 1.8.2. Sequencing runs were only then allowed into our data analysis work flow, if their statistics related to quality of base calling and to preliminary alignment against mouse genome were within three standard deviations of the mean of all runs. The five quality related statistics used were the percentage of clusters passed filtering (%PF clusters), mean quality score (PF clusters), percentage of reads aligned, mean alignment score and percentage of alignment error, and had the following population means and standard deviations, respectively, 93.39 ± 1.85, 36.34 ± 0.5, 72.24 ± 3.0, 290.6 ± 15.13 and 0.98 ± 0.15. Presence of adapter sequences in the datasets was checked with cutadapt utility (Martin, 2011). Less than 0.65% of the reads were estimated to include an adapter sequence longer than 10 bps; and more than 95% of such reads were left unmapped after the alignment of the reads to the mouse genome. Therefore, it was not necessary to trim reads to remove the adapter sequences.

Sequences were first mapped to local databases of mouse rRNA (Ensembl and NCBI) and rodent-specific repeat sequences (RepBase Update version 18.10 [Jurka et al., 2005]) using bowtie2 version 2.1.0 (Langmead and Salzberg, 2012) with default alignment parameters. Then reads were mapped to Genome Reference Consortium GRCm38 (mm10) version of mouse reference genome sequence using tophat version 1.4.1 with known mm10 transcripts provided via the—transcriptome-index option (Trapnell et al., 2009). For each read all mapping outcomes were considered and only non-rRNA, non-repeat reads that mapped uniquely to the mouse genome were selected for further analysis (Figure 2—source data 1).

Expression levels for mRNA and pre-mRNA were estimated per locus rather then per isoform. To this end, gene annotations from Ensembl mouse database release 68 (Flicek et al., 2013) were flattened similarly as described earlier (Anders et al., 2012) with the following modifications. Prior to flattening, transcript isoform models whose expression was not evidenced by count data or splice junction coverage were removed from the annotation database. This was achieved first by estimating the expression level of each isoform contained in Ensembl database by cufflinks software version 2.0.2 (Trapnell et al., 2010) with following parameters: -G Mus_musculus.GRCm38.68.gtf -u -j 0.8 -m 160 -s 50 -N. Isoform models which did not comprise at least 5% of the total expression of the locus or which were only expressed in fewer than three samples were excluded from the gene annotation database. Furthermore, spliced isoform models, whose splice junctions were not covered on average by 6 reads, or 40% or more of whose splice junctions were not covered at least by 6 reads were also removed from the annotation database. In addition to the removal of non-supported isoform models, novel exons overlapping with reference on the opposite strand were identified and added to the annotation database to improve the unambiguity of read counting. Briefly, cufflinks was run in reference annotation based transcript (RABT) mode (Roberts et al., 2011) and isoform models were assembled with cuffmerge software. Novel exon models with a class code ‘X’, which overlap with a reference exon on the opposite strand, were then validated by splice junction coverage as outlined before and added to the annotation database. Finally, the gene models were flattened (Anders et al., 2012).

An in-house Python script was used to count the reads mapped within each annotation feature in a similar way as implemented in htseq-count utility software (HTSeq: analysing high-throughput sequencing data with Python. [http://www-huber.embl.de/users/anders/HTSeq/↗]). Only reads which can be unambiguously identified as either exonic (continuous or spliced) or intronic for a single locus were counted towards the mRNA or pre-mRNA counts of that locus, respectively. Mappable and countable mRNA and pre-mRNA lengths (in bps) for each locus were calculated by means of generating all possible 100-bp long reads in silico (faux reads) for each transcript type and counting the faux reads through identical mapping and counting work flow used for real experimental reads.

Preliminary inspection of the extent of differential expression and the presence of highly expressed genes was carried out by cumulative percentage plots of raw counts; and accordingly read counts of mRNA and pre-mRNA datasets were normalised with upper quantile (Bullard et al., 2010; Dillies et al., 2013) and TMM (Robinson and Oshlack, 2010) normalisation methods, respectively. Prior to normalisation, transcripts which did not have at least 10 counts in at least one third of the samples were removed from the datasets. Differential expression (DE) analysis was performed using the DESeq analysis tool version 1.12.1 (Anders and Huber, 2010). Briefly, normalisation factors calculated earlier were used to create DESeq count datasets for mRNA and pre-mRNA counts annotated with full experimental design (treatment × time) and dispersion values were calculated with default settings. Calculated p-values for DE of all transcripts between control and Dicer KO conditions were then adjusted for false discovery rate (FDR) by Benjamini and Hochberg (BH) method (Benjamini and Hochberg, 1995). A transcript was considered differentially expressed if it had an adjusted p-value smaller than 0.05 and had a fold-change greater than 1.5.

For further filtering of lowly expressed transcripts and for better comparability between datasets, RPKM values were calculated as the number of counted reads per 1000 mappable and countable bases per geometric mean of normalised read counts per million. The geometric mean of normalised read counts was 53,692,658 and 16,285,588 for mRNA and pre-mRNA datasets, respectively. Transcripts which did not have a RPKM value greater than 0.1 or 0.01 in at least one third of the samples were considered as ‘not expressed’ in the mRNA or the pre-mRNA dataset, respectively.

As an estimate of the mRNA stability, we have calculated the log2 mRNA to pre-mRNA ratios of RPKM values separately for each sample. The distribution of log-ratios was inspected to assess normality. To test if a transcript's mean stability over all timepoints was increased in the absence of miRNAs, we applied a one-sided t test between control and KO conditions (n = 12). The p-values were then adjusted for FDR using BH method.

We used the combined power of a parametric test based on harmonic regression similar to those implemented in cosiner based algorithms (Cugini, 1993; Yang and Su, 2010) and a non-parametric test, JTK_CYCLE (Hughes et al., 2010). To adapt the harmonic regression to sequencing count data, several modifications were introduced. First, a negative binomial generalised linear model (GLM), with a log link function was used to regress the rounded normalised counts on a cycle with a 24-hr period. Count data were not detrended or smoothened before the regression. Dispersion values that were already fitted for each gene during DE analysis were used as an estimate of the shape parameter theta. From the fitted parameter values, then, amplitude and phase estimates were calculated for each gene. The goodness-of-fit for the harmonic model (sinusoidal expression) was tested with a likelihood ratio test against the updated null model including only the intercept parameter (constant expression). The analysis was performed using the MASS package (Rindskopf, 1997) in statistical computing environment R.

For both tests, the two biological replicates per timepoint were treated as independent replicates and were not combined into a false 48-hr data-series. The p-values calculated by each test were then adjusted separately for FDR by BH method (Benjamini and Hochberg, 1995). For all further calculations, amplitude and phase estimates obtained from the harmonic fit were used. We categorised mRNAs and pre-mRNAs as rhythmic when they passed at least one of the statistical tests and had an amplitude-ratio (peak/trough) ≥1.5. Combining two complementary tests in this way provided for relatively low stringency in rhythm detection, which we deemed important in particular for the pre-mRNA data, which is naturally noisier due to the low abundance of short-lived pre-mRNA molecules.

To assess the significance of the effects of the treatment on amplitude and phase estimates, we performed two permutation tests (N > 10,000). For both tests, we applied the permutation on the normalised read counts per gene, always maintaining the mRNA and pre-mRNA attributes. Each permuted data set was then analysed by exactly the same steps used on the experimental data. Finally, the distribution of the difference in means, significance of the effects, and correlation coefficients under the null hypothesis were inspected for both amplitude and phase effects. To enable the comparison of the significance of the effects between samples of different sizes, the Wilcoxon test statistics were converted into a standardised z-score.

In the first test, we permuted the timepoints while maintaining the treatment effects (Dicer KO vs control) to test the hypothesis that the observed phase delay differences in the stabilised circadian transcripts was not merely an outcome of altered expression levels (in particular lower amplitudes) in the Dicer KO (in which case it would be observable rather frequently on randomly generated cyclic transcripts). In the second test, the genotype (Dicer KO/control) was permuted to test the hypothesis that the observed effects on both amplitudes and phase differences were not random.

Statistical testing of enrichment of Gene Ontology terms within lists of genes was performed with MetaCore version 6.17 software suite (Thomson Reuters, New York City, NY). For enrichment analysis of miRNA targets within lists of genes, the ‘Conserved Family Info’ table from the TargetScanMouse database version 6.2 (Lewis et al., 2005) was used to create an association database between miRNA families and the genes from the current study. Using an in-house Python script, raw counts of targets and non-targets were extracted from this database for the genes included in list. Background list included all 11,841 genes from Figure 2F. The statistical significance of enrichment of targets for each miRNA family in the test list against the background list was tested via Fisher's exact test on the contingency tables created from the raw counts. Calculated p-values were adjusted for FDR using BH method (Benjamini and Hochberg, 1995). Analysis was carried out in R environment (R Core Team, 2013).

The RNA-seq data set produced in this study has been deposited in the Gene Expression Omnibus (accession number GSE57313↗). Supplementary files 1 and 2 are available at http://www.unil.ch/cig/page102471.html↗ and have been deposited in the Dryad Digital Repository (Du et al., 2014).

The following datasets were generated:

NH Du, AB Arpat, M De Matos, D Gatfield, 2014, Data from: MicroRNAs Shape Circadian Hepatic Gene Expression on a Transcriptome-Wide Scale, http://dx.doi.org/10.5061/dryad.cd726↗, Available at Dryad Digital Repository under a CC0 Public Domain Dedication.

NH Du, AB Arpat, M De Matos, D Gatfield, 2014, MicroRNAs Shape Circadian Hepatic Gene Expression on a Transcriptome-Wide Scale, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE57313↗, Publicly available at NCBI Gene Expression Omnibus.