Critical role of rhythmic poly(A) tail regulation in circadian gene expression

The circadian rhythmicity is deeply rooted in rhythmic regulation of gene expression. The core clock pathway that generates circadian oscillation has been well studied, but it remains unclear how this core rhythm controls rhythmic gene expression. Besides that several core clock components are transcriptional factors and mediate genome-wide rhythmic transcriptional control, mounting evidence has demonstrated the importance of rhythmic posttranscriptional controls. A recent study particularly highlighted rhythmic control of poly(A) tail lengths in hundreds of genes in mouse liver and a strong correlation of poly(A) rhythm with protein expression rhythm. In this work we constructed a simplistic model to study the effect of rhythmic poly(A) tail regulation on circadian mRNA expression. The model depicted rhythmic control imposed upon basic mRNA expression processes, including transcription, polyadenylation, deadenylation and degradation. The model results revealed rhythmicity in deadenylation as the strongest contributor to the rhythmicity in poly(A) tail length, and the rhythmicity in abundance of the mRNA subpopulation with long poly(A) tails, which serves a rough proxy for mRNA translatability. In line with this finding, our model further demonstrated that the rhythmic patterns found in the expression of deadenylases could funnel the peak phases of poly(A) tail length and long-tailed mRNA abundance, respectively, into three distinct groups, which could allow genes within each group to coordinate their functions around the clock. Last, our model suggested factors that contribute to the experimentally observed rhythmicity in poly(A) tail length and total mRNA abundance.


Introduction
Rhythmic control of gene expression is a hallmark of the circadian system. The daily rhythms in biochemistry, physiology and behavior ultimately stem from rhythmic gene expression in each cell (Takahashi, Hong et al. 2008, Mohawk, Green et al. 2012. In mammals, approximately 3-15% of mRNAs are rhythmically expressed with a ~24 hr period in any given tissue (Duffield 2003, Zhang, Lahens et al. 2014, Kojima and Green 2015. The rhythmicity originates from a cell-autonomous circadian clock machinery, which consists of a set of core clock genes interlocked by transcription-translation feedback loops (Zhang and Kay 2010, Partch, Green et al. 2014, Takahashi 2015. Many core clock genes encode transcription factors, and interact with their respective target enhancers to exert rhythmic transcriptional controls over circadian mRNA expression (Ueda, Hayashi et al. 2005, Takahashi 2015. While rhythmic transcriptional control has been extensively studied, rhythmic control of gene expression also occurs beyond the transcription process (Kojima, Shingle et al. 2011, Hirano, Fu et al. 2016, Gobet and Naef 2017. Recent genome-wide analyses and mathematical modeling particularly highlighted the role of post-transcriptional regulations in driving rhythmic mRNA expressions , Menet, Rodriguez et al. 2012, Luck, Thurley et al. 2014, Trott and Menet 2018, Wang, Symul et al. 2018. Post-transcriptional regulation targets various processes, such as RNA splicing, nuclear export, cellular translocation, dormancy and degradation of mRNA (Keene 2010). Many post-transcriptional processes are under circadian control, which, in turn, affect the phase and amplitude of mRNA abundance (Kim, Kim et al. 2005, Kojima, Gatfield et al. 2010, Kojima, Shingle et al. 2011, Du, Arpat et al. 2014).
Ultimately, rhythmicity in the transcription and post-transcriptional processes couple with each other and jointly determine the rhythmicity in gene expression. For example, circadian rhythmicity in RNA transcription and turnover (Luck, Thurley et al. 2014) jointly determine the rhythmicity in mRNA level. As yet, it remains poorly understood how rhythmicity in other posttranscriptional processes affect the rhythmicity in gene expression.
One of the post-transcriptional mechanisms that impact rhythmic gene expression is the regulation of poly(A) tail length. The tracts of adenosines at the 3' end of nearly all eukaryotic mRNAs are critical for controlling stability and translatability of the mRNA (Beilharz and Preiss 2007, Eckmann, Rammelt et al. 2011, Jalkanen, Coleman et al. 2014. Hundreds of mRNAs are recently discovered to exhibit robust circadian rhythms in their poly(A) tail lengths in mouse liver (Kojima, Sher-Chen et al. 2012). Interestingly, the rhythmicity in poly(A) tail length is closely correlated with rhythmicity in the level of the corresponding proteins, indicating that the rhythmic poly(A) regulation plays an important role in driving rhythmic protein expression (Kojima, Sher-Chen et al. 2012). Similar daily fluctuations in poly(A) tail lengths were also reported in the mouse brain (Robinson, Frim et al. 1988, Gerstner, Vanderheyden et al. 2012). In addition, the amplitude of mRNA rhythmicity increased in the absence of Nocturnin, a deadenylase (enzyme that removes poly(A) tails from mRNAs) which is rhythmically expressed in various mouse tissues (Wang, Osterbur et al. 2001, Stubblefield, Gao et al. 2018. These data collectively support the importance of poly(A) tail rhythmicity in regulating circadian gene expression.
In this work, we built a mathematical model that describes mRNA dynamics under the regulation of rhythmic transcription, polyadenylation, deadenylation and degradation. We used the model to systematically examine how coupling of rhythmic transcription and poly(A) tail regulation generates rhythmicity in poly(A) tail length and mRNA abundance. Our model results highlighted rhythmicity in deadenylation as the strongest determinant for rhythmicity in poly(A) tail length and in abundance of long-tailed mRNA. The latter can be regarded as a rough proxy for mRNA translatability because sufficiently long poly(A) tail is necessary for mRNA translation (Jacobson and Favreau 1983, Gallie 1991, Gallie 1998, Wells, Hillner et al. 1998, Weill, Belloc et al. 2012. As a corollary of this general finding, we demonstrated that the experimentally observed distinct phases in expression of the deadenylases (Kojima, Sher-Chen et al. 2012) can robustly group the rhythmic mRNAs by their peak phases in mRNA translatability (~long-tailed mRNA level), regardless of other rhythmic controls. This result implies that rhythmic deadenylation could provide a potential mechanism to orchestrate rhythmicity of gene expression across the genome. Finally, we used the model to predict the factors contributing to different classes of combinatorial patterns in the rhythmicity of transcription rate, poly(A) tail length and mRNA abundance (Kojima, Sher-Chen et al. 2012), which provides information for future experimental tests.

Model
In a typical RNA expression process, RNAs are first transcribed in the nucleus and acquire long poly(A) tails as a result of nuclear polyadenylation. After being exported into the cytoplasm, mRNAs undergo deadenylation and are ultimately degraded. Cytoplasmic polyadenylation, as another important post-transcriptional regulation, elongates the poly(A) tail to promote mRNA stability and translatability (Mendez and Richter 2001). Although cytoplasmic polyadenylation is typically associated with translational control in oocyte maturation, early embryo development and synaptic plasticity (Mendez and Richter 2001, Kwak, Drier et al. 2008, Charlesworth, Meijer et al. 2013, Cui, Sartain et al. 2013, it is recently suggested to play a role in circadian rhythmicity of mRNA (Kojima, Sher-Chen et al. 2012). The expression level of Gld2, a poly(A) polymerase responsible for cytoplasmic polyadenylation, was found to exhibit circadian rhythmicity (Kojima, Sher-Chen et al. 2012). Hence, in the model we incorporate polyadenylation, together with transcription, deadenylation and degradation, to capture the major processes that dynamically regulate poly(A) tail length and mRNA abundance ( Fig. 1). Instead of explicitly tracking the exact length of poly(A) tails, the model divides the mRNA population into a long-tailed fraction and a short-tailed fraction, which compares with the same division in the experiment (long-tailed >~100nt, short-tailed <~100nt, (Kojima, Sher-Chen et al. 2012)).
For the sake of simplicity, we made the following assumptions in the model based on experimental evidence. First, degradation only occurs to the short-tailed mRNAs, because mRNA degradation typically requires the poly(A) tail length to be shortened to 10~15 nt (Schoenberg andMaquat 2012, Wu andBrewer 2012). Second, we lump together the nuclear events, i.e., transcription and nuclear polyadenylation, because nuclear polyadenylation is rather ubiquitous (Proudfoot 2011) and the poly(A) polymerases responsible for nuclear polyadenylation have not been found to be rhythmically expressed (Kojima, Sher-Chen et al. 2012). In our model, therefore, the transcription step directly leads to a long-tailed mRNA, and the downstream cytoplasmic deadenylation and polyadenylation are explicitly described. The ordinary differential equations that govern the temporal dynamics of long-tailed ( ) and shorttailed ( ) mRNAs read as Eqs. (1) and (2).
Long-tailed mRNA: Short-tailed mRNA: To capture the circadian regulation of the four processes in Eqs.
where denotes the mean rate, the relative amplitude, and the peak phase, of the process labeled by the subscript. The angular frequency, , equals 2 (24hr) ⁄ . is a fixed parameter, and the other parameters take different values.
In this paper we focus on how rhythmicity in the four modeled processes affects the rhythmicity in the abundance of mRNAs and the ratio between long-tailed and short-tailed mRNAs, because these are quantities measured in the previous circadian transcriptome studies , Kojima, Sher-Chen et al. 2012. Furthermore, we are interested in the rhythmicity of abundance of long-tailed mRNAs, because elongation of the poly(A) tail activates translation (Jacobson and Favreau 1983, Gallie 1991, Gallie 1998, Wells, Hillner et al. 1998, Weill, Belloc et al. 2012. The abundance of the long-tailed subpopulation thus provides a rough proxy for mRNA translatability. From each simulation result based on Eqs. (1) and (2), i.e., time trajectories ( ) and ( ), we extracted the peak phases, relative amplitudes and mean levels of total mRNA abundance ( ( ) + ( )), long/short ratio ( = ( ) ( ) ⁄ ), and long-tailed mRNA abundance, ( ( )) (see Methods). These quantities were subject to further analysis, as elaborated in the following Results sections. For the rest of the paper, we will refer to these quantities, e.g., peak phase of LSR, generally as the "output quantities", unless any specific quantity is referred to.

Rhythmic deadenylation is the dominant contributor to rhythmicity in poly(A) tail length and long-tailed mRNA abundance
To investigate the effects of coupled rhythmicity in transcription, degradation, polyadenylation and deadenylation on the output quantities, we ran 100,000 numeric simulations of the model (Eqs. (1) and (2)) with random parameter values for the mean rates, relative amplitudes, and phases of each process: transcription, deadenylation, degradation, and polyadenylation ( Table 1 and Fig. S1), except for the mean rate of transcription, which only affected the overall abundance of mRNAs without any effect on the output quantities (see

Supplementary Materials).
Our model results revealed rhythmicity in deadenylation as the strongest contributor to the rhythmicity of poly(A) tail length, followed by rhythmicity in polyadenylation. Over the sampled parameter sets, the scatter plots of the simulation results demonstrated that the peak phases of LSR were strongly concentrated within a ~3 hr-wide stripe centered around 10 hrs before the corresponding input peak phases of deadenylation ( Fig. 2A), indicating a strong dependence of the LSR phase on the deadenylation phase. With a much more diffuse distribution of data points, the phase of LSR depended on the phase of polyadenylation much more weakly ( Fig. 2A). The phases of transcription and degradation had very little correlation with the phase of LSR ( Fig. 2A). We further performed variance-based sensitivity analysis (Sobol indices, (Sobol 2001, Iooss andLemaitre 2015)) to systematically quantify the impacts of each model parameter on each output quantity. Based on simulation results from a large number of random parameter sets, each Sobol index quantifies the impact of a particular model parameter on a particular output quantity, reporting the percentage of variance in the output quantity contributed by the variance in the model parameter. The estimated Sobol indices (Fig. 2B) corroborated with the dependencies indicated by the scatter plots. For example, the peak phase of deadenylation was suggested by the Sobol indices to be the strongest contributor to the peak phase of LSR.
Variance in the peak phase of deadenylation alone explained ~50% of variance in the peak phase of LSR (longest bar in Fig. 2B). Consistently, the relative amplitude of deadenylation was also found to be the strongest contributor to the relative amplitude of LSR (Fig. S2), even though the mean level of LSR depended nearly equally on the mean rates of deadenylation and polyadenylation (Fig. S2).
Our model results also showed a moderate impact of rhythmic deadenylation and polyadenylation on the rhythmicity of total mRNA abundance. Although rhythmicity of total mRNA abundance, as expected, was most correlated with rhythmicity of transcription and degradation, variances in the peak phases of deadenylation and polyadenylation still made a moderate contribution to the variance in the peak phase of total mRNA (Figs. 2C, D). This contribution stemmed from the regulation of mRNA stability by poly(A) tail length. In the model, this regulation is depicted by the simplistic assumption that degradation is restricted to the short-tailed mRNAs (Eqs. (1) and (2)). For the same reason, relative amplitudes of deadenylation and polyadenylation contributed to the relative amplitudes of mRNA to an extent right next to the mean rate of degradation (Fig. S4).
We further used the model to examine the effects of the four processes on the rhythmicity of mRNA translatability, using the abundance of long-tailed mRNA as a proxy. Although it reflects the multiplication of total mRNA level and LSR, the rhythmicity of long-tailed mRNA abundance relied most heavily on rhythmicity of deadenylation, followed by that of polyadenylation (Fig. 2F). The peak phases of long-tailed mRNAs were strongly concentrated in the ~3 hr-wide strip centered around 10 hrs before the peak phases of deadenylation (Fig. 2E), highly similar to that observed in the peak phases of LSR ( Fig. 2A). Consistently, the relative amplitudes of deadenylation and polyadenylation were also among the strongest contributors to the relative amplitude of long-tailed mRNA (Fig. S3). Hence, deadenylation and polyadenylation were found to be the strongest contributors to the rhythmicity of mRNA translatability.
Taken together, these results underscore the importance of deadenylation and polyadenylation in regulating the rhythmicity of poly(A) tail length and mRNA abundance (especially abundance of the long-tailed subpopulation). We particularly found deadenylation the dominant contributor to the rhythmicity of poly(A) tail length and long-tailed mRNA abundance.

Rhythmic deadenylation could robustly group genes by their poly(A) tail rhythms
Rhythmicity in transcription, deadenylation, polyadenylation and degradation of mRNA are ultimately controlled by rhythmicity in the abundance or activity of molecules mediating these processes, e.g., transcription factors, deadenylases and poly ( Intrigued by the above observation, in this section we use our model to explore the potential functional role of having diverse peak phases in deadenylases. First, we plugged in the model the observed distribution of transcription peak phases centered around ZT15 . In three separate in silico experiments, we chose degradation, deadenylation or polyadenylation, and let the process peak at three narrow windows centered around ZT0, 8, 16, and let the peak phases of the other two processes evenly distribute around the clock (2nd, 3rd and 4th columns in Fig. 3). The case with three deadenylation peak windows significantly grouped the peak phases of LSR and long-tailed mRNA abundance into three windows (3rd column in Fig. 3), whereas the cases with degradation or polyadenylation peak windows generated results were comparable to the case where all three processes have evenly distributed peak phases (1st column in Fig. 3), with wide distributions in all peak phases. To test the effect of the actual rhythmic patterns observed in deadenylases and polyadenylases, we further plugged in narrow peak phase windows centered around ZT3.5 for polyadenylation, and around ZT2, ZT5 and ZT13 for deadenylation. The simulations resulted in three highly distinct groups of mRNAs in their peak phases of poly(A) tail length and long-tailed mRNA abundance (5th column in Fig. 3).
These results corroborate with the strong sensitivity of rhythmicity in poly(A) tail length and long-tailed mRNA abundance to the rhythmicity in deadenylation (Fig. 2B, F). Note that the mean rates and relative amplitudes of all the processes assumed random variables as in the previous section (Table 1 and Fig. S1). Therefore, our results showed that multiple peak phases in deadenylation, but not other processes, could robustly group the rhythmicity of poly(A) tail length and mRNA translatability (represented by long-tailed mRNA abundance) into distinct peak phase windows, regardless of variations in the mean activity or rhythmicity of other processes.

mRNA half-life and other factors explain difference between Class I and Class II PAR mRNAs
In the previous transcriptome-wide study of rhythmicity in poly(A) tail lengths (Kojima, Sher-Chen et al. 2012), the mRNAs with rhythmic changes in the poly(A) tail length were grouped into three classes according to their rhythmic characteristics, i.e., whether 24-hour rhythmicity was detected in the levels of pre-RNA and total mRNA. The former is an indicator for transcription rate. The Class I mRNAs are rhythmic not only in poly(A) tail length, but also in pre-RNA and total mRNA (Fig. 4A). The Class II mRNAs are rhythmic in poly(A) tail length and pre-RNA, but not in total mRNA (Fig. 4A). The Class III mRNA are only rhythmic in poly(A) tail length, but not the other two (Fig. 5A). Factors like differences in mRNA half-lives were suggested in the original study to explain why the rhythmic patterns were different in these classes (Kojima, Sher-Chen et al. 2012). Here we leveraged our model to systematically identify factors that can lead to the combinatorial rhythmic patterns in these classes.
We first attempted to identify the model parameters that contributed most to the distinction between Class I and Class II. Because the only difference between Class I and II are the rhythmicity of total mRNAs, we focused on identifying model parameters that contribute most to the relative amplitude of total mRNA abundance. The Sobol indices revealed the mean rate constant of mRNA degradation as the strongest contributor to the amplitude of total mRNA (Fig. 4B). Indeed, the mRNA half-lives found in the cases with Class II behaviors exhibited a much wider distribution than those with Class I behaviors, extending towards larger values (Fig.   4C). This finding corroborates with the experimental observation that the average mRNA halflife is longer in Class II than in Class I (Kojima, Sher-Chen et al. 2012). After the mean rate constant of mRNA degradation, the Sobol indices also revealed the relative amplitudes of transcription and degradation, and the mean rate constants of deadenylation and polyadenylation, as contributors to the amplitude of total mRNA (Fig. 4B). Higher mean rate constant of deadenylation also contributed to Class I (Fig. 4D). This is related to the above finding about mRNA half-lives because deadenylation promotes degradation and hence increasing mean rate constant of deadenylation has a similar effect on mRNA turnover as increasing the mean rate constant of degradation. The Class I behaviors are also associated with stronger amplitudes in transcription and degradation (Fig. 4F, G). These factors were indeed shown to significantly affect amplitude of mRNA abundance in a previous theoretical study on coupling of rhythmic transcription and degradation (Luck, Thurley et al. 2014). Furthermore, the Class I behaviors were found to associate with antiphasic rhythms in transcription and degradation (Fig. 4H).
These additional factors differentiating Class I and Class II awaits future experimental testing.

Class III PAR mRNAs likely have longer mRNA half-life and non-rhythmic degradation
Compared with Class I and Class II, Class III does not have rhythmic transcription.
Because rhythmicity of transcription serves as input to our model, we cannot use the model to investigate the origin of lack of transcriptional rhythmicity. However, we were interested in understanding why all PAR mRNAs without transcriptional rhythmicity also lack rhythmicity in total mRNA level, and only exhibit rhythmicity in LSR (Kojima, Sher-Chen et al. 2012). For the sake of discussion, we use "Class IV" to refer to the non-existing group of PAR mRNAs that would exhibit rhythmicity in total mRNAs and poly(A) tails, but not in pre-RNA (Fig. 5A). We tried to identify model parameters that could make the difference between Class III and the nonexisting Class IV. Because both Class III and IV do not have rhythmic transcription, we ran 100,000 model simulations with non-rhythmic transcription (i.e., relative amplitude of transcription = 0, while keeping the other parameters sampled from the same distributions as before (Table 1 and Fig. S1). Without rhythmic transcription, the resulting Sobol indices again revealed the mean rate of degradation as the most sensitive contributor to the relative amplitude of mRNA abundance, which explained the difference between Class III and IV (Fig. 5B). The Class IV behaviors required much shorter half-life of mRNA (Fig. 5C). Therefore, the absence of Class IV indicates that the group of mRNAs without rhythmic control at the transcriptional level likely have long mRNA half-life. Indeed, Class III has the longest average mRNA half-life measured among all mRNAs that are rhythmically expressed (Fig. 2D in (Kojima, Sher-Chen et al. 2012)). We also found the relative amplitude of degradation as the next strongest factor that differentiates Class III from IV (Fig. 5B). The rhythmic characteristics of Class IV would require markedly higher amplitude of mRNA degradation (Fig. 5D). According to this result, the complete missing of Class IV PAR mRNAs implies that the group of mRNAs without transcriptional rhythmicity might also lack rhythmic degradation control, and hence all the mRNAs in this group exhibit the Class III, but not the Class IV behavior. Such RNAs likely depend on rhythmic poly(A) tail regulation to mediate downstream rhythmicity in protein translation. This conjecture awaits future experimental testing.

Discussion
Many posttranscriptional steps are shown to participate in circadian control of gene expression (Kojima, Shingle et al. 2011), among which rhythmic regulation of poly(A) tail length was recently discovered and characterized at the transcriptome level (Kojima, Sher-Chen et al. 2012). Notably, the rhythmicity in poly(A) tail is both the cause and result of rhythmic gene expression. On one hand, rhythmicity in poly(A) tail length controls rhythmicity in mRNA degradation (Luck, Thurley et al. 2014), because mRNA stability relies on poly(A) tail length (Jacobson and Favreau 1983, Wells, Hillner et al. 1998, Eckmann, Rammelt et al. 2011, Weill, Belloc et al. 2012). On the other hand, rhythmicity in poly(A) tail length itself stems from rhythmicity in different gene expression processes, such as polyadenylation and deadenylation.
Down to the mechanistic basis, it is the dynamic coupling between multiple rhythmic processes, i.e., transcription, polyadenylation, deadenylation and degradation, that generates the concomitant rhythmic patterns in poly(A) tail lengths and mRNA abundance. Such dynamic coupling strongly influences the rhythmic patterns in the gene product, as shown previously for coupling between rhythmic transcription and rhythmic degradation (Luck, Thurley et al. 2014).
In this work, we developed a simplistic mathematical model to study how rhythmic inputs from the transcription, degradation, polyadenylation and deadenylation processes control the output rhythmicity in mRNA abundance, poly(A) tail length and long-tailed mRNA abundance (as a proxy for mRNA translatability). Our model results highlighted rhythmic deadenylation as the strongest factor in controlling the peak phases and amplitudes of rhythmic poly(A) tail length and abundance of long-tailed mRNA (Fig. 2). This finding revealed a crucial role of rhythmic poly(A) regulation, especially rhythmic deadenylation, in rhythmic control of mRNA.
Following the identification of rhythmic deadenylation as the most important factor in controlling poly(A) tail length and long-tailed mRNA abundance, we used the model to test the potential effects of having three distinct peak phases in deadenylases as those suggested by previous experiments (Kojima, Sher-Chen et al. 2012). According to the model results, these three groups of deadenylases, coupled with the experimentally observed single peak phase of poly(A) polymerase (Kojima, Sher-Chen et al. 2012) and concentration of transcription peak phases around ZT15 , can robustly divide peak phases in poly(A) tail length and abundance of long-tailed mRNAs into three distinct groups, regardless of further rhythmic control (e.g., by rhythmic degradation) (Fig. 3). This finding indicates that rhythmic deadenylation may not only diversify the rhythmic patterns in mRNA functionality, but also organize the rhythmic genes into groups with close phases. Such grouping could potentially provide an effective mechanism to coordinate expression of genes that synergistically mediate physiological functions at distinct times of the day. Genes related to the same cellular or physiological functions could potentially be regulated by the same deadenylase, leading to similar phases in their expression and subsequent synergistic functioning. Our model further confirmed that such grouping effect is not significant if the same three peak phases are associated with the polyadenylation or degradation processes, instead of the deadenylation process. In other words, the diversity of rhythmicity in deadenylases is more effective than diversity in other rhythmic processes in imposing rhythmic grouping on mRNAs.
Our model suggested factors that likely contribute to the rhythmic characteristics of three observed classes of PAR mRNAs. Generally, multiple factors could affect the rhythmic characteristics of an mRNA, reflecting complexity and flexibility of the mechanisms controlling circadian gene expression. In principle, four classes of PAR mRNAs should exist, depending on whether total mRNA and preRNA are rhythmic or not (Fig. 5A). Only three classes, however, were reported previously (Kojima, Sher-Chen et al. 2012). For PAR mRNAs without rhythmicity in preRNA, they ubiquitously lack rhythmicity in total abundance of mRNA. In addition to the previous suggestion of longer half-life (Kojima, Sher-Chen et al. 2012), our model analysis suggested non-rhythmic mRNA degradation as another possible contributing factor to this phenomenon (Fig. 5B, D). It is up to future experiments to test the effects of factors predicted by the model that lead to different classes.
The finding that rhythmic deadenylation plays a dominant role in regulating rhythmicity of poly(A) tail and mRNA translatability poses two interesting questions. First, could deadenylation help bridge the circadian rhythm with other rhythmicity? For example, genes related to poly(A) tail regulation, such as Parn, one of the rhythmically expressed deadenylases, are suggested to gain rhythmicity from clock-independent signals (Wang, Symul et al. 2018).
Equipped with a strong impact on mRNA rhythmicity, such deadenylation process could potentially help couple the circadian rhythm to other rhythmic signals, e.g., feeding, to optimize physiological functions around the clock. Second, how do different deadenylases interact with each other in the regulation poly(A) tail rhythmicity? In other words, could deadenylases that peak at different phases regulate the same mRNA and impose a more complex rhythmic pattern?
These questions need future studies to resolve.

Model simulation and extraction of phase, amplitude and mean from simulation results
For any given parameter set, Eqs. (1) and (2) were simulated using the ordinary differential equation (

Parameter sampling
We performed global parameter sensitivity analysis (Saltelli 2008) on the model to analyze the general contribution of each parameter to any specific output quantity (i.e., peak phase, relative amplitude and mean of ( ) + ( ), ( ) ( ) ⁄ and ( )). To perform such global sensitivity analysis, one needs to simulate the model with randomly chosen parameter values that maximally represent the parameter space. In this study we drew random parameter values from the distributions listed in Table 1 and plotted in Fig. S1. The peak phases and relative amplitudes were sampled form uniform distributions of their possible ranges by definition (Table). The mean reaction rates were sampled from log-normal distributions suggested by previous genomic scale measurements (see sources indicated in Table 1). It is noteworthy that we set the mean transcription rate as constant, as it only causes proportional changes in ( ) and ( ), and does not affect the rhythmic patterns of any quantity (see Supplementary Materials). To improve the accuracy of the global sensitivity analysis for models with many parameters, one needs sampled parameter values that well represent the parameter space. To this end, we used the sampling method of Sobol sequence, which was known to ensure good representation of the high-dimensional parameter space (Iooss and Lemaitre 2015).

Sensitivity analysis
To evaluate the impact of each model parameter (e.g., phase of deadenylation) on each model output (e.g., relative amplitude of LSR), we used a variance-based global parameter sensitivity analysis method, the Sobol indices (Iooss and Lemaitre 2015). For each pair of parameter and output in the model, the first-order Sobol index, , characterizes the contribution of variance in alone to the total variance in (Eq.(4)). The total-effect index, , characterizes the contribution of variance in , as well as the variance caused by its coupling with other parameters, to the total variance in (Eq. (5)). The larger these indices are, the more sensitive is to , or the more impact has on .
where the subscript ~ indicates all indices except for .
The calculation of variance was slightly modified for the peak phases, the values of which live on a circular domain. This results in a sudden jump of the numeric values of the phase at the ends of the circular domain, e.g., 24hr to 0hr, even though it reflects a continuous change. This problem can artificially increase the estimated variance. To correct the problem, the difference between phases, 1 − 2 , was replaced by its sine value, sin ( 1 − 2 ), in the evaluation of the Sobol indices.