Circadian clock-dependent and -independent posttranscriptional regulation underlies temporal mRNA accumulation in mouse liver

Abstract

The mammalian circadian clock coordinates physiology with environmental cycles through the regulation of daily oscillations of gene expression. Thousands of transcripts exhibit rhythmic accumulations across mouse tissues, as determined by the balance of their synthesis and degradation. While diurnally rhythmic transcription regulation is well studied and often thought to be the main factor generating rhythmic mRNA accumulation, the extent of rhythmic posttranscriptional regulation is debated, and the kinetic parameters (e.g., half-lives), as well as the underlying regulators (e.g., mRNA-binding proteins) are relatively unexplored. Here, we developed a quantitative model for cyclic accumulations of pre-mRNA and mRNA from total RNA-seq data, and applied it to mouse liver. This allowed us to identify that about 20% of mRNA rhythms were driven by rhythmic mRNA degradation, and another 15% of mRNAs regulated by both rhythmic transcription and mRNA degradation. The method could also estimate mRNA half-lives and processing times in intact mouse liver. We then showed that, depending on mRNA half-life, rhythmic mRNA degradation can either amplify or tune phases of mRNA rhythms. By comparing mRNA rhythms in wild-type and Bmal1^-/- animals, we found that the rhythmic degradation of many transcripts did not depend on a functional BMAL1. Interestingly clock-dependent and -independent degradation rhythms peaked at distinct times of day. We further predicted mRNA-binding proteins (mRBPs) that were implicated in the posttranscriptional regulation of mRNAs, either through stabilizing or destabilizing activities. Together, our results demonstrate how posttranscriptional regulation temporally shapes rhythmic mRNA accumulation in mouse liver.

Keywords: RNA binding proteins; circadian clock; mRNA half-lives; posttranscriptional regulation.

Fig. 1.

Kinetic model identifies contributions and parameters of rhythmic transcription and rhythmic degradation regulating mRNAs from total RNA-seq. (A) Temporal accumulations of mRNA and pre-mRNA from time-resolved total RNA-seq were used to fit four kinetic models (M1–M4). The rate equation for the temporal accumulation of mRNA $m\left(t\right)$ depends on pre-mRNA $p\left(t\right)$ and the degradation rate $\gamma \left(t\right)$, which are either constant or parameterized with periodic functions (Materials and Methods). The models allow for constant (C) or rhythmic (R) synthesis (S) and degradation (D) in the four combinations: M1 (constant synthesis and constant degradation, CS-CD), M2 (rhythmic synthesis and constant degradation, RS-CD), M3 (constant synthesis and rhythmic degradation, CS-RD), and M4 (rhythmic synthesis and rhythmic degradation, RS-RD). Probabilities for each model are estimated using Schwarz weights (Materials and Methods); the optimal model for each mRNA (one per gene; Materials and Methods) also yields gene-specific parameters (e.g., mRNA half-life, processing time, phases, and amplitudes of rhythmic degradation rates). We applied this approach genome-wide. (B–D) Temporal profiles of mRNA and pre-mRNA of genes assigned to models M2–M4. Data for mRNA (blue) and pre-mRNA (green) with error bars (SE over four biological replicates) are shown as relative expression, that is, the total normalized counts divided by the average value over time. Solid curves are the fitting for the optimal model and estimated parameters (blue, mRNA; green, pre-mRNA; red, degradation). Peak times (phase) and amplitudes are summarized in circle plots. The radial scale of these plots is relative to the largest relative amplitude of mRNA, pre-mRNA, or rhythmic degradation. Absolute half-lives [$\log \left(2\right)/{\gamma }_0$, in hours], if identifiable, are labeled. (B) M2: (Left) Cp had long estimated hl (7.9 h), which damped amplitude of mRNA compared with that of pre-mRNA; (Right) Rorc mRNA was identified in M2 (RS-CD) with estimated constant hl of 2.1 h. (C) M3: Fus mRNA was identified in M3 (CS-RD). The peak time of rhythmic degradation (RD) was ZT18.3 and the relative amplitude of RD was 0.3. Mean half-life was nonidentifiable (Materials and Methods). (D) M4: (Left) Per3 mRNA was identified in M4 (RS-RD). The RD showed a maximum at ∼ZT18, and a relative amplitude of 0.5 mean degradation rate was identifiable with mean hl of 1.6 h; (Right) Cbs mRNA showed a phase delay between mRNA and pre-mRNA > 6 h, which could be explained by M4. Parameters of RD showed a maximum at ∼ZT9 with a mean hl of 8.6 h.

Fig. 2.

Rhythmic mRNA degradation regulates 35% of rhythmically accumulating mRNAs (∼20% in M3 and 15% in M4). (A) Numbers and percentages of rhythmic mRNAs (6,014 with FDR < 0.05, harmonic regression) identified in M2 (3,949), M3 (1,167), and M4 (898). (B, Left) Percentages of mRNAs subjected to RD estimated using the PA test, which were assigned to either M3 or M4 by our method. The analysis is stratified in function of the stringency of the PA test (log₁₀ P values). We used published measurements of half-lives from cell lines as input for the PA test. (Right) Distribution of P values from the PA test according to the classification from our method. (C and D) Distributions of peak times (C) and amplitudes (D) of rhythmic mRNAs in M2–M4.

Fig. 3.

Distributions of estimated half-lives and processing times of mRNA in mouse liver. (A) Distribution of estimated mRNA half-lives [$\log \left(2\right)/{\gamma }_0$, in hours] in M2 and M4 for transcripts with identifiable half-lives. Vertical lines indicate the medians of the distributions. (B) Distribution of pre-mRNA processing times ($1/k$, in minutes) in models M2 and M4. Vertical lines indicate the medians of the distributions for transcripts with identifiable parameters. (C) Half-lives and processing times of core clock and clock output mRNAs. We plotted standard errors (SEs) of the estimates when available (SI Materials and Methods and Dataset S1).

Fig. 4.

Phases and amplitudes of rhythmic mRNA degradation. (A and B) Distribution of estimated phases and peak-to-trough amplitudes (in log₂) of rhythmic degradation in M3 (A) and M4 (B) for transcripts with identifiable and high-confidence degradation parameters (coefficient of variation of estimated relative amplitude <0.4 and SE of estimated phase <1 h; SI Materials and Methods and Dataset S1).

Fig. 5.

Rhythmic mRNA degradation serves as an amplifier for long-lived mRNAs and phase tuner for short-lived mRNAs. (A, Left) Histogram of phase delays in M4 (peak time of RD minus peak time of RS). (Right) Phase delay vs. mRNA half-lives suggests three classes of transcripts. C1, Long-lived mRNAs (hl > 5 h; green); C2, short-lived mRNAs with short delay (hl < 5 h, delay < 12 h; orange); and C3, short-lived mRNAs with long delay (hl < 5 h, delay > 12 h; dark red). (B, Left) Histogram of relative amplitude ratio $\left({\epsilon }_{\gamma }/{\epsilon }_s\right)$ between RD and RS. (Right) Scatterplot of relative amplitude ratio vs. mRNA half-life color-coded for C1–C3. (C) Scatterplot of relative amplitudes [(maximal expression − minimal expression)/2 × (mean of expression)] of mRNAs in C1 vs. predicted under assumption of constant degradation (Left) and vs. relative amplitudes of pre-mRNAs (Right) in log scale. (D) Scatterplot of peak times of mRNAs in C1 vs. predicted peak times under assumption of constant degradation (Left) and vs. peak times of pre-mRNAs (Right). (E, Left) Temporal profiles of pre-mRNA (green) and mRNA (blue) of C1 gene Smagp (values are normalized read counts divided by the temporal average; error bars show SE over four biological replicates). The gray line shows predicted mRNA profile under assumption of constant degradation (relative amplitude of mRNA degradation ${\epsilon }_{\gamma }$ set to zero). (Right) Arrows in the circle plot depict the phase (angular coordinate) and relative amplitude (radial coordinate) of pre-mRNA (green), mRNA (blue), and predicted mRNA under the assumption of constant degradation. (F–H) Idem as C–E for genes in C2. Rhythmic degradation advances peak times of mRNAs in C2 without affecting relative amplitude. Pre-mRNA and mRNA of Wee1 oscillates with comparable phases, whereas assuming constant degradation, its mRNA would show larger phase delays. (I–K) Idem for genes in C3. Rhythmic degradation delays peak times of mRNAs in C3 without affecting relative amplitude. Pre-mRNA and mRNA of Slc4a4 show large delays, whereas assuming constant degradation, its mRNA would show smaller phase delays (K). (A–K) Transcripts were selected as for Fig. 4, with, in addition, a threshold on the rhythmicity of their pre-mRNA and mRNA to ensure reliable phase estimation (FDR < 0.05 rhythmicity test, and relative amplitude >0.1; Materials and Methods).

Fig. 6.

BMAL1-dependent and -independent rhythmic mRNA degradation. (A) Heat map of mRNA relative expression of genes classified in M3 in mouse liver in the following conditions (biological replicates are averaged): WT ad libitum (AL WT, Right; this condition was used for the main analysis, 12 times points, 48 samples), WT restricted feeding (RF WT, Middle; 6 time points, 12 samples), and Bmal1^−/− restricted feeding (RF Bmal1^−/−; 6 time points, 12 samples). All data are from ref. and provided in Dataset S5, together with the parameter estimates. Genes above the thickest horizontal separator have similar rhythmic parameters in AL WT and RF WT. The two thinner horizontal separators split the genes into those that keep the same rhythmic parameters in RF Bmal1^−/− (Top, 569 mRNAs, group 1), those that are expressed in a constant manner (Bottom, 292 mRNAs, group 2), and those that exhibit 24-h oscillations with altered rhythmic parameters (Middle, 59 mRNAs). (B) mRNA peak time distributions for groups 1 and 2. (C) Phases of inferred activities of mRBP motifs in groups 1 (blue) and 2 (red). Activities are inferred using a penalized linear model that integrates mRBP binding site at the 3′-UTR of mRNA with phase and amplitude of mRNA degradation (Materials and Methods). Angles on the circle indicate peak activity times (ZT times). (D) Inferred activities and protein abundance of candidate mRBPs over time (50). mRBPs are predicted to regulate stability of transcripts from group 1 (plotted in blue), from group 2 (in red), or in some cases from both groups (in purple). Rhythmic activities that are antiphasic to rhythmic protein accumulation predict a stabilizing effect.