Critical role of deadenylation in regulating poly(A) rhythms and circadian gene expression

The mammalian circadian clock is deeply rooted in rhythmic regulation of gene expression. Rhythmic transcriptional control mediated by the circadian transcription factors is thought to be the main driver of mammalian circadian gene expression. However, mounting evidence has demonstrated the importance of rhythmic post-transcriptional controls, and it remains unclear how the transcriptional and post-transcriptional mechanisms collectively control rhythmic gene expression. In mouse liver, hundreds of genes were found to exhibit rhythmicity in poly(A) tail length, and the poly(A) rhythms are strongly correlated with the protein expression rhythms. To understand the role of rhythmic poly(A) regulation in circadian gene expression, we constructed a parsimonious model that depicts rhythmic control imposed upon basic mRNA expression and poly(A) regulation processes, including transcription, deadenylation, polyadenylation, and degradation. The model results reveal the rhythmicity in deadenylation as the strongest contributor to the rhythmicity in poly(A) tail length and the rhythmicity in the abundance of the mRNA subpopulation with long poly(A) tails (a rough proxy for mRNA translatability). In line with this finding, the model further shows that the experimentally observed distinct peak phases in the expression of deadenylases, regardless of other rhythmic controls, can robustly cluster the rhythmic mRNAs by their peak phases in poly(A) tail length and abundance of the long-tailed subpopulation. This provides a potential mechanism to synchronize the phases of target gene expression regulated by the same deadenylases. Our findings highlight the critical role of rhythmic deadenylation in regulating poly(A) rhythms and circadian gene expression.

. Overview of the study. (A) Schematic diagram of the model. The model describes four processes that control the poly(A) tail length and mRNA abundance: transcription, degradation, cytoplasmic deadenylation and polyadenylation. The rhythmicities of the four processes, i.e., amplitude and phase, are presumably controlled by the core clock mechanism (shaded molecular circuit), which is not explicitly included in the model.  3) It is impressive that some rich information can already deduced from the model given by equations 1 and 2. It is justified to start with such simple models, which are convenient for thorough parameter space examination. There are, however, certain limitations on coarsegraining the poly-A tail length as two discrete states. Mathematically one can treat the process as a reaction-diffusion process, and can write down either a discrete chemical master equation or continuous 1D process (described by Langevin equations or Fokker-Planck equations ). The authors may discuss limitation of their formulation and alternative future studies.

Response:
We thank the reviewer for pointing out the limitation of coarse-graining the poly(A) tail length as two states and suggesting future study. In the revised manuscript, we have added the descriptions of this limitation and potential future study on p. 21: "In our current model, the poly(A) regulation has been coarse-grained as one-step conversions between a long-tailed and a short-tailed mRNA subpopulation. Such coarse-graining retains the most essential kinetic features of the poly(A) regulation processes, while allowing for significant reduction of the model and efficient global parameter sensitivity analysis. From such analysis we identified the critical role of deadenylation in rhythmic regulation. In reality, both deadenylation and polyadenylation act sequentially, i.e., adding or subtracting one adenosine at a time. Unlike one-step chemical reactions, the kinetics of sequential processes are often non-exponential (Deneke et al., 2013;Eisen et al., 2019). To evaluate the rhythmicities of poly(A) regulation and gene expression more accurately, we will include a linear reaction chain in the model to account for sequential steps of deadenylation and polyadenylation in our future work." The manuscript titled "Critical role of deadenylation in regulating poly(A) rhythms and circadian gene expression" investigates the role of rhythmic deadenylation in regulating circadian gene expression by performing the variance-based global sensitivity analysis. Furthermore, the authors provided the potential mechanism to explain why mRNAs are grouped into three classes based on the presence of rhythms in pre-mRNA and total mRNA. Considering the impact of the result and novelty in the method, I strongly believe that the manuscript fits well with PLOS Computational Biology.

Response:
We thank the reviewer for taking the time to review our manuscript and for the valuable and constructive comments and suggestions. We truly appreciate that the reviewer points out the significance and novelty of our work, as well as its fit with PLOS Computational Biology.

Response:
We fully agree with the reviewer that slow mRNA degradation tends to result in Class II. However, degradation rate is not the only factor distinguishing between Class I and Class II, as the distribution of degradation rate significantly overlap between the two classes ( Fig. 4C). For degradation rates that are shared by both classes, other factors, such as the phase difference between transcription and degradation ( Fig. 4D), and the amplitudes of transcription and degradation ( Fig. 4E-F), also play significant roles.
Then one might ask, why are Class I parameter sets significantly enriched in certain ranges for the factors in Fig. 4D-F, whereas Class II parameter sets show nearly even distributions in those factors? This question can be more easily understood by "reconstructing" the parameter distributions in the global parameter space from the 1D distributions ( Fig. R1 below). In the global parameter space, the Class I parameter sets are distributed in a confined region ( Fig. R1, purple), whereas Class II parameter sets are broadly distributed in the remaining region not taken by Class I. Because the Class II region is broad, the 1D projections of this region appear nearly even. This led to our original interpretation of Fig. 4D-F (currently on p. 16): "generation of significant rhythmicity in L+S (Class I) requires sufficient phase difference between transcription and degradation, and sufficiently high amplitudes of transcription and degradation, simultaneously. If any of these conditions are not satisfied, total mRNA abundance would not have significant rhythmicity (Class II)." From an alternative perspective, the three factors in Fig. 4D-F essentially compensate each other in producing "no rhythm" in total mRNA. A similar argument applies to Class III/IV for the interpretation of   4. Case (ii): The parameter sets in case (i) that satisfy −1.15 ≤ log 10 dgrd ≤ 0. Case (iii): The parameter sets in case (i) that satisfy −2 ≤ log 10 dgrd ≤ −1. 15. As the mean degradation rate, dgrd , decreases, there are fewer Class I parameter sets, and they are distributed in a more confined region.  5. Case (ii): The parameter sets in case (i) that satisfy −1 ≤ log 10 dgrd ≤ 1. Case (iii): The parameter sets in case (i) that satisfy −1.5 ≤ log 10 dgrd ≤ −1. As the mean degradation rate, dgrd , decreases, there are fewer Class IV parameter sets, and they are distributed in a more confined region. Response: We thank the reviewer for pointing out the transferability of the Sobol's method to other chronobiology studies. We have added a discussion about the advantage of the Sobol's method in analyzing chronobiology models on pp. 22-23:

The method used in the manuscript (Sobol index) would be useful to investigate
"In addition, the Sobol's method serves as a particularly powerful tool for parameter sensitivity analysis for models in chronobiology. In chronobiology models, oscillation phases are often important quantities of interest. As circular variables, i.e., ZT 0 = ZT 24, phases are intrinsically nonlinear and non-monotonic. Analyzing nonlinear and non-monotonic variables using classic correlation and dependency analyses, such as Pearson correlation and Spearman correlation, could lead to misleading conclusions, because these methods are based on assumptions about linear and monotonic relations between the analyzed data. For example, in our model, Pearson and Spearman correlation analyses demonstrate strong negative correlation between the phases of deadenylation and L/S ratio, a spurious conclusion due to the circular nature of phases (Fig. S6); other pairs of input and output phases suffer different levels of distortion in their Pearson and Spearman correlations (Fig. S6). Based on variance decomposition (see Methods), the Sobol's method circumvents these problems and can effectively analyze nonlinear and non-monotonic variables (Iooss and Lemaitre, 2015). The method can be used widely in chronobiology models to identify the key factors that drive phases of target quantities, such as the phase difference between PER2 and TP53, whose interaction is critical for the crosstalk between the circadian clock and cell cycle ( We also added Fig. S6 to compare the Sobol's method with Pearson and Spearman analyses, using our model as an example.

More descriptive subtitle would be better instead of "Model"
Response: The subtitle has been changed to "Model for rhythmic mRNA and poly(A) tail regulation" on p. 6.

More detailed explanation for Eq. (4) and (5) would be helpful for reader.
Response: The Sobol's method is now explained in more details on pp. 26-27 to facilitate readers' understanding. With one equation added to illustrate the conceptual basis of the Sobol's method, Eqs. (4) and (5) is the contribution from the -th parameter alone. Here , is the contribution from the interactions between the -th and -th parameters, where ~ denotes the combined parameter set except for the -th and -th parameters. Contributions from higher-order interactions between parameters are defined similarly as ( ).
The Sobol indices are then defined as fractions of the decomposed terms in Eq.
(2) out of the total variance, Var( ). In practice, only the single (Eq.(2)) and total-effect (Eq.(3)) indices are calculated because relatively simple algorithm as described below can be designed. Specifically, the single Sobol index, , characterizes the contribution of variance in alone to the total variance in (Eq. (2)). The total-effect, or total 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.(3)). Conveniently, the total-effect contribution equals  ( ) in Equation (7) seems ,1 not ,2 .

Response:
The wrong subscript has been corrected. Figure S1.

Labels (i.e. A, B, C, …) are missing in
Response: Panel labels have been added in Figure S1.

Is it true that deadenylation is the strongest contributor for S and S/L?
Response: Yes, deadenylation is the strongest contributor for rhythmicities in S and S/L. We have calculated the Sobol indices for S and S/L (Fig. R2 below). The amplitudes and phases of the two quantities strongly depend on the rhythmicities of deadenylation and polyadenylation, similar to L and L/S. Note that S/L is just the reciprocal of L/S, Hence, the corresponding Sobol indices for S/L are expectedly almost identical to those for L/S. Interestingly, the mean level of S is almost exclusively determined by the mean degradation rate (Fig. R2). This can be understood from the steady-state level of S when all four input processes are non-rhythmic. In this non-rhythmic case, = trsc dgrd ⁄ only depends on the rate constants for transcription and degradation. Hence, in the rhythmic model, the mean S is expected to mainly depend on the mean rates of transcription and degradation. Because the mean transcription rate is set to be 1 in our study (as it does not affect the rhythmic patterns in any case), the mean S strongly depends only on the mean degradation rate. Following a similar reasoning, mean L depends on mean deadenylation and polyadenylation rates (Fig. S2), because in the non-rhythmic case, = trsc ( polyA + dgrd ) deA dgrd ⁄ depends on polyA and deA . However, in the model without cytoplasmic polyadenylation ( polyA = 0), L depends only on mean deadenylation rate (Fig. S3), because the corresponding non-rhythmic case has = trsc deA ⁄ .

Page 5: …recently …. 2012???
Response: The word "recently" has been removed on ln. 80. It has also been removed or replaced on ln. 20 and 117, when referring to the same paper.

Response:
We chose numerical solutions, because approximate analytical solutions were too complex and intractable, as the reviewer suspected. At the beginning, we did use variation of constants to find the approximate analytic solution for ( ) and ( ). Specifically, we plugged the following ansatzes into Eqs. (1) and (2), and solved for the coefficients.
Unfortunately, the analytical solutions are rather complex and intractable. For example, the peak phase of long-tailed mRNA, Since these analytic solutions did not provide much insight about the behaviors of the system, we resorted to numeric solutions.
We have added the following clarification for the choice of method on pp. 8-9.
"Because the parameters of the model are largely unknown and likely vary significantly from gene to gene, we need to investigate the dependency of output rhythmicities on input rhythmicities in the global parameter space (i.e., the entire possible range of parameter values).
In previous studies, such dependency was analyzed by deriving approximate analytic solutions to models with up to two rhythmic input processes ( Motivated by the reviewer's comment, we also added a brief statement in the model setup section to clarify the expected high variations of model parameters across different genes. This is on p. 7: "Note that the four processes can assume different amplitudes and phases for different genes, because these regulations can be mediated by different

Response:
We suspect that the reviewer meant L+S instead of S, as Fig. 2 did not show any data for S. To answer the reviewer's question about the effect of deadenylation on S, we calculated the Sobol indices for S ( Fig. R3 below), and found that deadenylation also strongly controls the rhythmicity of S, as the reviewer expected. In the manuscript, we only presented the results for L, because poly(A) tail facilitates translation initiation and in this work we use L as a rough proxy for mRNA translatability.  Figure 1, I am not surprised that Deadenylation plays a major role. It has a central position and modulates the decay of L and the production of S. If this is the correct interpretation, I would expect for constant DeA a major control of Polyadenylation on the L/S ratio since this reaction is also at the center of the chain and effects L and S.

Response:
We agree that polyadenylation lies at the center of the reaction chain, affecting both long-tailed and short-tailed mRNA like deadenylation. We hereby calculated the Sobol indices when deadenylation is non-rhythmic. In this case, the rhythmicity of L/S ratio is indeed mostly controlled by polyadenylation, as the reviewer expected (Fig. R4 below). Our main results in the paper, interestingly, showed that when both polyadenylation and deadenylation are rhythmic, deadenylation is stronger than polyadenylation in regulating the rhythmicities of L/S ratio (Figs. 2A, 2B, 3(iv), 3(v)).