Statistical characterization of circadian transcripts

To assess whether genes are regulated by the circadian oscillator, most methods take advantage that circadian regulation of transcript abundance resemble a sinusoid. To estimate circadian period of the regulation of a particular transcript, the main idea is to find the sinusoid that most closely matches its abundance over time [32,33]. However, in NAM-treated plants the changes in abundance of circadian-regulated transcripts have a considerable number of non-sinusoidal profiles (S1 Fig). To overcome this problem, we devised a learning approach based on pseudo-sinusoidal functions to properly assess the rhythmicity and the corresponding circadian period of signals from gcRMA normalized microarray data of NAM treated plants. To infer period, phase and amplitude, linear trends are eliminated by removing the best straight-line fit and pseudo-sinusoidal functions are fitted to each signal to minimize the 2-norm error. Pseudo-sinusoidal functions account for many signals that are periodic but not sinusoidal. Pseudo-sinusoidal functions are constructed by joining together two sinusoids with different periods. Hence, a complete oscillation of a pseudo-sinusoidal function consists of the first sinusoid (of period p₁) in the first half-oscillation, and the second sinusoid (of period p₂) in the second half-oscillation (Fig 2A). The resulting period of the pseudo-sinusoidal function is defined as . This can be expressed by: where A is a scaling factor that accounts for the amplitude of the signal and φ₁ is the phase of the signal. The algorithm searches possible combinations of p₁ and p₂ to minimize the least square distance between pseudo-sinusoidal functions and the data. We allowed periods p₁ and p₂ to vary between 12 and 36 hours. A perfect sinusoid gave a high fit for the wild-type background dataset. We found that three periodic signals were highly represented in the dataset. In particular, those with p₁, p₂ equal to: p/2, p/2 (pure sinusoid); p/2+3.8, p/2–3.8 (p1 is greater than p2); and p/2–7.3, p/2+7.3 (p1 is smaller than p2) (Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The effects of NAM on the circadian regulation of the transcriptome.

(A) Illustration of the shape of S. The first panel shows two period of a perfect sinusoidal shape, whereas the second panel displays the segmentation of the period P into p1 and p2, where p1 is greater than p2. p1 and p2 follows the formula: P = (p1 + p2)/2. The last panel displays the case were p1 is smaller than p2. (B) Number of periodic transcripts that have been identified in untreated and NAM-treated plants, as well as the intersection. (C) Circadian period of untreated and NAM-treated transcripts plus minus standard deviation. The mean increase of period following the NAM treatment is of 3.3h. (D) Amplitude analysis (normalized) for the same transcripts.

https://doi.org/10.1371/journal.pcbi.1006674.g002

We used a logistic regression framework to generate a probabilistic discriminative model that estimates the probability of a gene to be rhythmic given its time course data. In this case, the classification problem only contains two classes: rhythmic (C₁) and arrhythmic (C₂). For each transcript, a set of 8 features x = {X₁,X₂,…,X₈} is computed and empirically believed to be crucial to distinguish between rhythmic and arrhythmic transcripts.

The features were computed from 2 signals: the first signal (A) corresponds to the average of replicates and (B) being a single replicate for which the L2-norm error with the best fitted pseudo-sinusoidal function is lower than for the other replicate. The following features were computed: ratio of power in the 18–32 hours frequency range (of (A) and (B)), L2-norm of the error to the best fit of pseudo sinusoidal function (of (A) and (B)), the variance of the power spectrum (of (A) and (B)) and the amplitude of the best fitted pseudo-sinusoidal function (of (A) and (B)).

The log of the ratio of probabilities between the two classes, also known as the log odds, is given by [34]:

The goal of the logistic regression is to estimate σ for a linear combination of the X_n features such that:

The weights b_i of the independent variables X_i were estimated using the mnrfit function in MATLAB. The algorithm is initially trained with a mix of 100 rhythmic and 100 arrhythmic transcripts randomly chosen from the dataset and visually inspected to show clear (ar)rhythmicity. Finally, the decision boundary was set so that if p(C₁|x)>0.5, the gene was classified as rhythmic, and vice versa. Our approach, therefore, is inspired by the patterns observed in the dataset but not strictly constrained to pure cosine shapes. With the inclusion of the S function, we allow the search for asymmetric signals, which represent a large part of the transcriptome. A main distinction with the previously introduced algorithms, therefore, is the data-specific, learning approach devised to allow for a wider range of periodic signals. However, this offers additional advantages such as a dedicated way to handle noise between replicates, or the information in the frequency domain of the signal, which are both learned from the data. Comparison of performances with standard periodicity assessment tools is shown on S2 Fig.

Network inference and analysis by DyDE

Like most biological systems, circadian clocks have a large number of feedback loops. Hence, a perturbation anywhere in the network typically affects all nodes (in this case, their molecular concentration and time profiles), which makes the problem of inferring the entry point of a perturbation hard using DE analysis. We proposed, instead, that key mechanisms involved in NAM-induced perturbations in the circadian system of Arabidopsis can be captured by identification and comparison of regulatory dynamics before and after the perturbation occurred. Assume that a perturbation, such as NAM, changes the regulatory dynamics between two genes (e.g. by binding to a transcription factor) while leaving intact the rest of the system. Due to feedback interconnections, all the clock genes would change their expression, which, in turn, would change the expression of all circadian genes. While thousands of genes change their expression, only one regulatory link was actually affected. Our goal is to find this link (or links, in case of multiple perturbation entry points). To achieve this, we developed DyDE that looks for changes in links, instead of nodes. DyDE uses Linear Time-Invariant (LTI) models, a type of black box model, to systematically capture the dynamics underlying the biochemical mechanisms of circadian gene regulation, without relying on a priori knowledge of the system or extensive database. Such models benefit from a rich theory and a well-established collection of tools that makes the analysis of its dynamical properties straightforward, as contrast to detailed mechanistic models. In addition, the estimation of the parameters of such models is reliable and computationally efficient. The description of biological mechanisms of the Arabidopsis circadian clock from time-series data by LTI models has been studied in [35]. More recently, the performances of such linear modelling approach to reverse engineer the clock topology were compared for two extensively used Arabidopsis oscillator models [36] and confirmed that the majority of oscillator links can be represented by simple linear dynamics. However, the use of LTI models to detect dynamical perturbation in the gene regulatory network resulting from chemical treatments is novel.

The first step of DyDE consists of uncovering dependencies and quantifying dynamics between genes with LTI models. Our mathematical framework estimates a collection of Single Input-Single Output (SISO) models between pairs of genes to characterize the system dynamics. The limited number of available time points restricted the modelling of SISO systems to first and second order models. Overall, second order systems did not improve significantly the fitness of models and resulted in a considerable increase of false positives (overfitting). Hence, in this analysis of the circadian system, we restrict the model order to one. Mathematically, the dynamics between two genes can be represented as: where u(t) and y(t) represent the time series of the regulatory gene and the regulated gene, respectively. In addition, b y(t) corresponds to the degradation rate of gene y, a u(t) corresponds to the influence of u(t) on the rate of y(t) and c is a constant offset. System identification is performed using the function ‘pem’ implemented in MATLAB to minimize the prediction error [37]. The model has a total of three parameters (a, b, and c), leading to efficient solutions. We chose a subspace initialization algorithm since it performed similarly as randomizing initial conditions–for the vast majority of models (99%), the final solution was identical with either method. This suggests that the chances of being trapped into a local minimum are negligible.

The estimation of parameters requires low computational time: a single system between a pair of genes is typically identified within few seconds (Intel Core i5). This modeling is independently repeated for all available pairwise genes, where each gene takes its turn as being an input and then an output to another gene. This modeling approach, therefore, generates a large amount of SISO LTI models (n²−n models, where n corresponds to the amount of genes, and self-regulation is not considered) to describe the system. Each potential link between two genes is validated if the corresponding model reproduces the dynamics involved with a sufficient degree of precision, which is characterized by a high goodness of fit, defined as: where y is the validation data, is the average value of the validation data, and is the estimated output. MATLAB function compare can be used to compute the fitness of the model. A fitness equal to 100% corresponds to a perfect identification. The choice of such metric is motivated by the dependency of noise towards the abundance of gene expression. When the distance of the true data points towards the mean is large (represented by the denominator in the above equation), the fitness conveniently penalizes less the error term, which lies in regions where the intrinsic noise involved in the gene expression is potentially the largest.

The second step consists in identifying the effect of a treatment, NAM in our case, on the biological network. While a treatment might affect the abundance of many transcripts, only a few links are affected, as depicted in red in Fig 1. Hence, checking whether links are affected before and after perturbation can potentially lead to the finding of the entry point of the treatment. For this purpose, two cases are of particular interest. First, a link between two genes is validated in the untreated system alone (i.e. it is not possible to find a combination of a, b and c so that the model in the treated system provides a good match with the data anymore). Second, a link is validated in both systems, but the way one gene regulates the other may change; this is a much subtler change in the dynamics of the link. The latter case requires us to compare the dynamics between both links. Here, we use a rigorous and well-established tool from engineering known as the nu-gap [38]. Originally developed to address the stability properties of closed loops systems defined in the same feedback loop, the nu-gap essentially measures the distance, from a perturbation point of view, between linear models. This property is particularly relevant in the context of circadian clock networks, which consist in regulatory networks with several feedback loops. This then facilitates us to determine the significance of the dynamical change of a link between experimental conditions. The nu-gap returns a value between 0 to 1, quantifying whether the models are similar or very different, respectively. [39] have suggested that values above ~0.2 could be used to infer the main target of a perturbation. The nu-gap is computed using the gapmetric function in MATLAB. It should be applied to all models that have been estimated in both networks. If the signals are concentrated around a particular range of frequencies (such as oscillating signals), the gap should be measured ‘locally’ around that range of frequencies only, since they dominated the model estimation in Step 1.

Next, we explain the key ideas behind DyDE through a small number of genes in the Arabidopsis circadian oscillator. For example, the following model considers TOC1 as an input and PRR9 as an output. where b represents the strength of activation or repression induced by TOC1 on the expression rate of PRR9, and a corresponds to the degradation rate of PRR9. These parameters are estimated by minimizing the prediction error from the untreated time-series for both TOC1 and PRR9.

In this case, we found a model in good agreement with the data (57% fitness), suggesting that indeed TOC1 regulates PRR9 (Fig 3A). Moreover, the model demonstrates that the rate of change of the concentration of PRR9 is proportional to the concentration of TOC1. Note that the other way around (i.e., PRR9 regulating TOC1) could not be established since the respective model has a low goodness of fit (16%, Fig 3B). These results are consistent with the literature [40]. Hence, we would then establish a link from TOC1 to PRR9, but not the other way around (Fig 3C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Network inference and analysis by dynamical differential expression (DyDE).

(A) Ordinary Differential Equations (ODEs) capture the dependence of the rate of the concentration of a transcript on the concentration of another transcript. First order linear models are used to represent the dynamics between two genes. Here, a good agreement (plain line) with the data (dotted line) was found (57% goodness of fit). (B) The inverse regulation is considered. In this case, it is not possible to find a combination of parameters so that a first order linear model captures the dynamics involved. For this inverse regulation the model that best described the data obtained a goodness of fit of only 16%. (C) A threshold by which each model is (in)validated is applied on the goodness of fit of the models. As an example, a threshold of 46% would consider a link from TOC1 to PRR9 but not the other way around. The same threshold is applied to all models. (D) A first order linear model is evaluated in the presence of nicotinamide between the same species. The nu gap is then applied to compare models (A) and (D) to quantify whether the models are similar, or significantly affected by NAM.

https://doi.org/10.1371/journal.pcbi.1006674.g003

Then, a model is estimated between TOC1 and PRR9 from the NAM-treated time-series. From the untreated and treated time-series alone, it is unclear whether the link dynamics have changed (Fig 3D). The optimal model parameters, however, have significantly changed. A nu-gap of ~0.5 confirms that indeed the link has been affected. This result indicates that there is large perturbation in the regulatory dynamics that links TOC1 to PRR9, which, therefore, should be considered as a strong candidate for being an entry point for NAM in the system. If true, knocking down either TOC1 or PPR9 would therefore lead to NAM no longer affecting the clock. This analysis is then repeated for all common links between untreated and treated plants.

DyDE applied to the Arabidopsis circadian clock genes

We considered a total of 17 known clock genes: CCA1, LHY, PRR9, PRR7, PRR5, RVE8, GI, TOC1, ZTL, ELF4, ELF3, PHYA, PHYB, CRY1, CRY2, CHE and PRR3. However, the core oscillator genes ZTL, ELF3, PHYB, CRY1, PRR3 and CHE were identified as non-rhythmic in the presence of NAM, which was confirmed by visual inspection (S1 Fig). Hence, these genes are excluded from the modeling of NAM targets as they cannot be contributing to the rhythmic dynamics of the remaining oscillator components that are measured in the presence of NAM.

As a first step, we computed models for all available pairs of the clock genes for both conditions, totaling 220 SISO models (110 in untreated and 110 in NAM). We kept only those models with good agreement with the data, i.e. above a fitness threshold. On one hand, the user-defined threshold has to be set large enough to reliably capture the dynamics involved between genes, and provide the nu-gap analysis with comparable models. On the other hand, the threshold has to be set sufficiently low to consider enough gene-to-gene relationships to detect a dynamical perturbation in the network. Here, the fitness threshold was set to 46% as we noted that below this threshold, the amount of unknown regulations dramatically raised (S3 Fig; S3 Table).

In total, 70 regulatory links were retained for untreated plants and 55 links for NAM-treated plants between the 11 clock genes. The untreated models describe 70% of the known regulatory pathways among these 11 genes (S3 Table; S3 and S4 Figs [40]). 64% of which, had the expected activation or inhibition effect. These numbers are remarkable, taking into account the model simplicity, and confirms that the majority of clock links can be represented by simple linear dynamics [35,41,42].

In particular, 28 links were present in the untreated samples but not in the NAM-treated samples. These 28 links form a network from now on referred to as “regulation loss” network, which captures the links abolished by NAM. In addition, 42 links are present in both conditions which form a network, so called “common” network that is common to both treated and untreated plants (S3 Table).

We used the nu-gap to identify those links among the common network whose dynamics were significantly affected by NAM. Fig 4A and S4 Table depict the comparison of the dynamics of each link with the nu-gap. All regulatory interactions are somehow affected by the treatment, which is expected from the interconnected circadian network. Let us then consider the highest nu-gap values, which are associated with the following links: TOC1 to PRR9 (0.5), those originating from CRY2 to ELF4 (0.47), LHY (0.42) and RVE8 (0.37) and PRR9 to CRY2 (0.35). Interestingly, the only inferred interaction originating from CRY2 that does not seem affected connects to TOC1 (nu-gap of 0.06). These results suggest that a major dynamical change is induced to CRY2 in the dynamical response of the circadian clock to NAM. In addition, the largest nu-gap value suggests that the causality within the time course data of TOC1 and PRR9 has changed significantly differently towards the treatment, as compared to the other parts of the circadian network.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. DyDE applied to the Arabidopsis circadian oscillator genes.

(A) Common network and nu-gap analysis. The common network displays the models that have been validated in both untreated and NAM-treated plants. A directed arrow from gene a to gene b (blue circles), therefore, represents a dynamical model that captures the dependency of b on a. Red arrows represent the models associated with the top five highest nu-gap values. (B) Bar plot comparing the connectivity loss (%) associated to each gene. For a particular gene, the connectivity loss corresponds to the total amount of incoming and outgoing links that were validated in untreated plants but not in NAM-treated plants. Error bars represent the standard deviation of connectivity loss for ± 5% change in fitness threshold selection.

https://doi.org/10.1371/journal.pcbi.1006674.g004

We then used a standard network topology metric to identify the genes that are central to the drastic changes in dynamics captured by the regulation loss network. This topology metric accounts for the connectivity of a gene, i.e. the number of its incoming and outgoing links. This measure is estimated for each gene of the regulation loss network. As an example, PRR7 has six incoming links and nine outgoing links for untreated plants. The connectivity of PRR7 in untreated plants is then equal to 15. Among those, only six of were present in NAM-treated plants. PRR7, therefore, has a connectivity of nine in the regulation loss network, which correspond to a loss of 60% of its connectivity from untreated to NAM treated plants. As a result, CCA1 (61%), PRR7 (60%), TOC1 (57%) exhibit the highest connectivity drop (Fig 4B; S5 Table). This result identifies the biological functions of CCA1, PRR7 and TOC1 as being highly affected by NAM in the regulation of the circadian clock.

DyDE, therefore, identifies the regulatory dynamics of TOC1-CRY2-CCA1-PRR7 as being predominately affected by NAM as a result of both nu-gap and connectivity analysis. Accordingly, the strong emergence of the blue light receptor CRY2 in the nu-gap analysis suggests that nicotinamide alters the regulation of the interactions between light signaling and the circadian oscillator. These findings are further examined through mutant analysis and single wavelength light experiments.

PRR7/PRR9 inter-regulation together with TOC1 are targets of nicotinamide

To test the predictions that TOC1, CRY2, CCA1 and PRR7 are associated with the effect of NAM on the circadian oscillator, we experimentally investigated the sensitivity of circadian mutants to NAM. All mutants responded to NAM with increased circadian periods, with the exception of two independent lines of the same T-DNA insertion allele of PRR7, which were insensitive (prr7-3 p > 0.95; prr7-11 p > 0.95 Fig 5A and 5B; S5 Fig; S6 Table). The insensitivity of prr7-11 to NAM was confirmed by measuring circadian rhythms of leaf movement (S6 Fig). prr7-11 was not affected by NAM at any tested concentration, contrasting with a dose-dependent effect of NAM on circadian period in other prr mutants and associated backgrounds (R² > 0.9; Fig 5C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Reverse genetic analysis validates the prediction that TOC1 and PRR7 are associated with the effect of nicotinamide on circadian period.

(A) The change in circadian period caused by 20 mM NAM (period difference) in circadian mutants measured using CCA1:LUC, TOC1:LUC, and CAB2:LUC. (B) CCA1:LUC activity for Col-0 and prr7-3 in the presence (yellow) and absence (white) of 20 mM NAM. (C) Dose response of circadian period to NAM for prr7-11, prr9-10, prr7-3prr9-10, prr3-1, prr5-1 and toc1-2. Bars ±SD. N >16 from > 2 technical replicates. Open symbols indicate minus NAM. Mean time courses for these data are shown in S5 Fig. Statistical analysis is detailed in S6 Table.

https://doi.org/10.1371/journal.pcbi.1006674.g005

In contrast, toc1-2 and TOC1-ox had significantly greater responses to NAM than wild type (Fig 5A; S6 Table). These results support our predictions that NAM induces dynamical changes specifically to PRR7 and TOC1. No dramatic changes of period, however, were observed for cry2-1 and cca1-11, suggesting that these might not contribute directly to the response to NAM.

Finally, derived from the nu-gap analysis, the possible change in the dynamical behavior of PRR9 in mediating the effect of NAM on the clock was evaluated with a prr7-3 and prr9-10 double mutant. prr7-3 and prr9-10 had an epistatic interaction, with the double mutant responding to NAM by a 5.3 ± 1.6 h increase of period, more than either the insensitive prr7-3 or the oversensitive prr9-10 alone (Fig 5A). The epistasis of prr9-10 to prr7-3 was confirmed at all concentrations of NAM tested (Fig 5C).

Nicotinamide-induced changes in period are associated with a blue light signaling pathway

The mutant analysis did not confirm the modeling dynamical perturbation of CRY2 in the response to NAM. However, CRY2 is one of a pair of cryptochrome blue light photoreceptors and so mutant analysis might not be the most appropriate tool to investigate the role of the blue light photoreceptor. To investigate further we also investigated the role of blue light in the response to NAM using monochromatic light conditions. High frequency measurements of the circadian promoter fusions PRR9:LUC, PRR7:LUC, TOC1:LUC, CCA1:LUC, LHY:LUC and GI:LUC were collected in the presence or absence of 20 mM nicotinamide under constant blue or red light (S7 Fig).

In the absence of blue light, NAM was without effect on the circadian period or amplitude of CCA1:LUC (Fig 6A) and other promoter:luciferase fusions (Fig 6B). This demonstrates that input pathways associated with blue light are sensitive to NAM. Under blue light exposure, all promoter:luciferase fusions considered had an increase in period in the presence of NAM (Fig 6B). Under red light exposure, the period response was either negligible (PRR9:LUC, CCA1:LUC, LHY:LUC, GI:LUC) or negative (PRR7:LUC, TOC1:LUC). These results suggest that blue light increase the response of circadian period to NAM, while red light decrease its responsiveness.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. The effect of nicotinamide on circadian period requires blue light.

(A) CCA1:LUC rhythms in monochromatic blue or red light ± 20 mM NAM. (B) The change in circadian period ± 20 mM NAM in constant white (grey) and monochromatic red or blue light for a range of reporters. (C) 20 mM NAM abolishes circadian rhythms of [Ca²⁺]_cyt in both Col-0, prr7-11 and prr7-prr9. (D) Robust circadian rhythms of [Ca²⁺]_cyt in monochromatic blue light are abolished by 20 mM NAM. (E) Red light induced elevations of [Ca²⁺]_cyt early in the photoperiod are not abolished by nicotinamide. (F) PHYB-ox enhances circadian rhythms of [Ca²⁺]_cyt in monochromatic red light. Light conditions are indicated by the colored boxes on the X axes. White is red/blue mix, black is dark, monochromatic light is represented by the appropriate color with subjective night shaded darker than subjective day. NAM indicated by yellow. Bars are SD. n > 7.

https://doi.org/10.1371/journal.pcbi.1006674.g006

Having previously proposed that the effects of NAM on the circadian system are due to the inhibition of the production of the Ca²⁺-agonist cADPR [12], we tested if the response to NAM of prr7-11 is due to altered Ca²⁺ signaling. We investigated, therefore, the inhibitory effects of NAM on circadian [Ca²⁺]_cyt oscillations in prr7-11 and in light signaling mutants in red and blue light. 20 mM NAM was equally effective in abolishing circadian rhythms of [Ca²⁺]_cyt in both Col-0, prr7-11 and prr7-3 prr9-10 (Fig 6C). This suggests either that there are multiple sites of action of NAM or that PRR7 is downstream of the effects of NAM on [Ca²⁺]_cyt.

In constant blue light, there were robust oscillations of [Ca²⁺]_cyt in plants with functional CRY1 photoreceptors, being abolished in cry1 and, cry1cry2 but unaffected by cry2, phototropins and Phy loss-of-function mutants (Fig 6D, S8 Fig). Under blue light, NAM abolished [Ca²⁺]_cyt oscillations but did not reduce oscillations further in cry1 or cry1cry2 (Fig 6D, S9 Fig). High amplitude oscillations of [Ca²⁺]_cyt were dependent on blue light because in constant red light, [Ca²⁺]_cyt increased early in each cycle but without a subsequent decrease (Fig 6E; S8 Fig). This red light-induced increase in [Ca²⁺]_cyt was dependent on PHYB (S8 Fig).

To examine the role of PHYB further we measured [Ca²⁺]_cyt in PhyB-ox (Fig 6F, S8 Fig) and determined that in these plants [Ca²⁺]_cyt was rhythmic with a sinusoidal period of 25.0 ± 0.5 h in constant red light (Fig 6F). NAM was without effect on [Ca²⁺]_cyt in constant red light, even in the PHYB-ox background (Fig 6E and 6F, S6 and S7 Figs) demonstrating that blue light regulates circadian oscillations of [Ca²⁺]_cyt through a NAM-sensitive pathway. This pathway appears to be required for the major oscillatory dynamics of [Ca²⁺]_cyt.

Extension of DyDE to the rhythmic transcriptome

DyDE was further adapted to explore the rhythmic genome for additional targets for NAM and novel clock genes. For this purpose, models were computed between each pair of the 988 genes that were scored rhythmic in both untreated and NAM treated conditions, resulting in 2 million models corresponding to potential interactions.

We selected the models that exhibit the highest goodness of fit (over 80%) in both untreated and NAM-treated plants to minimize the identification of erroneous interactions and computed their nu-gap value to investigate dynamics affected by NAM. As a result, out of ten, two models only were retained with a nu-gap > 0.2. These models identified the regulation of AT5G35970 (P-loop containing nucleoside triphosphate hydrolases superfamily protein) by AT2G21860 (violaxanthin de-epoxidase-like protein) and the regulation of ATG21660 (GRP7/CCR2) by AT1G78600 (LZF1/BBX22) as being altered by NAM. The regulation of AT5G35970 by AT2G21860 may be important as AT5G35970 is identified by DyDE as being a hub regulated by four circadian oscillator genes (S8 and S9 Tables). The second link is easier to explain because GRP7 along with GRP8 forms a slave oscillator driven by the circadian clock that regulates ABA responses [43]. GRP7 is an RNA binding protein regulated by ADP ribosylation [44]. As the enzymes that perform ADP ribosylation are inhibited by NAD, this could suggest a role for nicotinamide inhibiting ADP ribosylation of an oscillator or slave oscillator component.

Then, the fitness threshold was released to 60% to further investigate novel clock components. For this purpose, we searched for those genes for which models can be computed from/to clock components (S8 Table). Models with a nu-gap value above 0.2 were discarded as a consistency criterion. Finally, candidates were ranked according to their connectivity with the clock. As a result, 20 high potential genes were isolated (S9 Table). The whole genome analysis of clock input and output hubs and the nu-gap analysis suggest interesting roles for previously characterized genes, including AT3G47500 (CYCLING DOF FACTOR3) [45], AT4G38960 (BBX19) [46], AT1G78600 (BBX22) [47,48], AT3G22840 (CRY3) [49], AT1G28330 (DRM1), AT2G33830 (DRM2) [50,51] and uncharacterized genes including AT5G35970.

Experimental procedures