{'fraction_non_alphanumeric': 0.038778539981570596, 'fraction_numerical': 0.023369509573052113, 'mean_word_length': 4.65234375, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Redox rhythm reinforces the circadian clock to gate immune response

Recent studies have shown that in addition to the transcriptional circadian clock, many organisms, including Arabidopsis, have a circadian redox rhythm driven by the organism’s metabolic activities1–3. It has been hypothesized that the redox rhythm is linked to the circadian clock, but the mechanism and the biological significance of this link have only begun to be investigated4–7. Here we report that the master immune regulator NPR1 (non-expressor of pathogenesis-related gene 1) of Arabidopsis is a sensor of the plant’s redox state and regulates transcription of core circadian clock genes even in the absence of pathogen challenge. Surprisingly, acute perturbation in the redox status triggered by the immune signal salicylic acid (SA) does not compromise the circadian clock but rather leads to its reinforcement. Mathematical modelling and subsequent experiments show that NPR1 reinforces the circadian clock without changing the period by regulating both the morning and the evening clock genes. This balanced network architecture helps plants gate their immune responses towards the morning and minimize costs on growth at night. Our study demonstrates how a sensitive redox rhythm interacts with a robust circadian clock to ensure proper responsiveness to environmental stimuli without compromising fitness of the organism.

Life on Earth has evolved the circadian clock to anticipate diurnal and seasonal changes 8 . This "scheduling" mechanism coordinates biological processes to reduce random energy expenditures and increase fitness. In Arabidopsis, the daily time-keeping is driven by three interlocked transcription-translation feedback loops (TTFLs): the core loop, the morning loop, and the evening loop. The core loop consists of three transcription factors (TFs): two partially redundant morning-phased CIRCADIAN CLOCK ASSOCIATED 1 (CCA1) and LATE ELONGATED HYPOCOTYL (LHY), and the evening-phased TIMING OF CAB2 EXPRESSION 1 (TOC1). CCA1/LHY and TOC1 are repressors of each other's expression 9,10 . Besides the TTFL circadian clock, non-transcriptional redox oscillations exist in all domains of life, including Arabidopsis 3 . Even though redox rhythm was shown to influence the TTFL clock 3 , how these two oscillatory systems are linked molecularly and what the biological significance of having two oscillatory systems is remain largely unknown.
To begin addressing these questions, we examined the daily changes in the reductionoxidation coenzymes NADPH and NADP + in Arabidopsis under constant light and found them to display circadian rhythms (P < 10 −4 ), with NADPH peaking before subjective dawn and NADP + peaking before subjective dusk (Fig. 1a, b). Moreover, their ratio also oscillated in a circadian manner (Extended Data Fig. 1). These data support the existence of widespread metabolic and redox rhythms in plants beyond the previously reported oscillations of oxidized peroxiredoxin, H 2 O 2 and catalases 3,7,11 . It is known that the plant immune-inducing signal salicylic acid (SA) can alter the cellular redox to trigger defence gene expression 12 . We found that under constant light, treating plants with SA could significantly perturb NADPH and NADP + rhythms as well as their ratio (Fig. 1a, b and Extended Data Fig. 1), indicating that the redox rhythm is sensitive to external perturbations.
We next examined whether this SA-triggered redox rhythm perturbation could be transduced to the circadian clock by first focusing on the evening-phased TOC1, which is responsive to many environmental factors 13 . Using quantitative PCR (qPCR), we observed significant increases in amplitude and average expression of TOC1 upon SA treatment (Fig. 1c). Similar results were observed using a transgenic line carrying a reporter of the TOC1 promoter fused to luciferase (TOC1p:LUC) 14 (Fig. 1d, e and Extended Data Fig. 2a). Strikingly, the period of the TOC1p:LUC expression rhythm did not change regardless of whether SA was applied at subjective dawn ( Fig. 1d) or dusk (Fig. 1e).
To study the effect of endogenous SA, which oscillates in a circadian manner 15 , on the clock, we crossed the TOC1p:LUC reporter into the SA biosynthesis mutant, sid2 (SA induction-deficient 2) 16 . We found that the amplitude and the average expression of TOC1 were significantly reduced in sid2 and this phenotype was rescued upon exogenous SA treatment (Extended Data Fig. 2b). Our results indicate that endogenous SA plays a part in the redox rhythm that modulates the amplitude and average expression of the circadian clock.
SA-induced redox changes can lead to reduction of the master immune regulator, NPR1, the release of NPR1 monomer for nuclear translocation, defence gene induction 12 and subsequent degradation mediated by the nuclear SA receptors NPR3 and NPR4 17 . To test whether the SA-mediated regulation of TOC1 is through NPR1, we crossed TOC1p:LUC into the npr1 mutant 18 . We found that the mutation not only dampened the basal expression of TOC1 but also abolished the SA-triggered increases in expression regardless of the time of treatment ( Fig. 2a and Extended Data Fig. 3a-c).
We hypothesized that NPR1 is an intrinsic regulator of TOC1 in response to the rhythmic accumulation of the endogenous SA 15 . Through western blotting, we indeed found a circadian oscillatory pattern for the NPR1 monomer (P < 0.01) with a peak at night (Fig. 2b and Extended Data Fig. 4a). Therefore, oscillation in the endogenous SA level may drive the rhythmic nuclear translocation of NPR1 to regulate the circadian clock genes. To test this hypothesis, we used mutants of cytoplasmic-localized thioredoxins (TRX), trx-h3 and trx-h5, in which NPR1 nuclear translocation is largely impaired 19 . We found that both the basal rhythm of TOC1p:LUC and its responsiveness to SA were diminished in trx-h3 trx-h5 ( Fig.  2c and Extended Data Fig. 5a), suggesting the requirement of NPR1 nuclear translocation in regulating TOC1 expression. Besides SA, glutathione reduced ethyl ester (GSHmee), a redox-altering reagent 20 , could also enhance TOC1 expression in an NPR1-dependent manner (Extended Data Fig. 5b), suggesting that NPR1 is a general redox sensor in modulating this clock gene.
NPR1 is a transcription cofactor of the TGA class of TFs in SA-induced defence gene expression 21 . Using yeast one-hybrid assay, six Arabidopsis TGAs were found to have strong binding affinities to the TOC1 promoter at the two TGA-binding sites (TBS) (Fig.  2d). To confirm this in planta, we mutated TBS in the TOC1p:LUC reporter (TOC1p (TBSm):LUC) and transformed it into Arabidopsis. We found that these mutations significantly inhibited transcription of the reporter (P < 0.001), indicating that TGAs are transcription activators of TOC1 (Fig. 2e). A direct role that NPR1/TGA plays in regulating TOC1 expression was further confirmed through chromatin immunoprecipitation (ChIP) in which association of NPR1 to TBS in the TOC1 promoter was significantly enhanced upon SA induction (Fig. 2f).
TOC1 is unlikely the only clock gene regulated by NPR1, because lowering the TOC1 level shortens the clock period while elevating the level lengthens the period 22,23 . However, no such perturbation was observed in npr1 (Extended Data Fig. 3d) or after SA treatment (Fig.  1d, e). Moreover, SA treatment at dawn should have caused an immediate induction in TOC1 expression instead of a 12-hour delay (Fig. 1d). To systematically search for other NPR1-targeted clock genes, we performed mathematical modelling using the P2012 circadian model 13 under the assumption that NPR1 is also a transcriptional activator of other clock genes (X and Y in Fig. 3a) (see Methods and Extended Data Figs. 6 and 9 for details).
We first optimized the P2012 model to fit the TOC1 expression in npr1 (Fig. 2a), which was a single parameter fit (i.e., basal expression in the absence of functional NPR1). The heat map of the best least-squares fit showed a characteristic "crosshair" pattern centred on PSEUDO-RESPONSE REGULATOR 7 (PRR7) (Fig. 3b), indicating that the basal regulation of PRR7 by NPR1 best explains the unchanged TOC1 period in npr1 (Extended Data Fig.  6a, b). This prediction was verified using qPCR in which PRR7 transcript levels in npr1 were found to be significantly lower compared to wild type (Fig. 3c). The second fitting for SA-induction data involved multiple parameters. We used our fixed basal expression parameter and NPR1 western data ( Fig. 2b and Extended Data Fig. 4) to fit the TOC1 expression from Fig. 2a. The resulting heat map showed a "crosshair" pattern for LHY/CCA1 (Fig. 3d), suggesting that either one or both of these genes is responsive to SA through the function of NPR1. Using qPCR, we found that while CCA1 and EARLY FLOWERING 3 (a negative control) did not respond to SA, the amplitude of LHY expression was significantly elevated by SA (P < 0.05) (Extended Data Fig. 6c-e) as predicted by our model. This result was further confirmed using the LHYp:LUC reporter (Fig. 3e). Consistently, the amplitude of basal LHY expression was reduced in npr1 (P < 0.05) while that of CCA1 remained unchanged (Extended Data Fig. 6f, g). Because LHY is an antagonist of TOC1 in the clock, induction of LHY by SA explains the delayed increase in TOC1 after SA treatment at dawn (Extended Data Fig. 6i, j) when LHY has its highest expression. This balanced network architecture of NPR1 regulating both the morning-phased LHY and the evening-phased TOC1 (Fig. 3a) strengthens the clock when the redox rhythm is perturbed.
To investigate the effect of reinforced circadian clock on plant immunity, we examined SAinduced resistance against bacteria in a toc1 mutant and found it to be more sensitive to induction than wild type and npr1 (Fig. 4a). While TOC1 negatively regulates resistance against bacteria, CCA1 and LHY have been reported to positively regulate resistance against bacteria and oomycetes 24,25 , timing immunity for the morning when temperature and humidity are the most favourable for infection 25 . We hypothesized that SA/NPR1-mediated induction of both morning and evening components of the circadian clock plays a major role in maintaining this diurnal difference in plant sensitivity to pathogen challenge.
To test this, we first examined induction of WRKY40, a direct target of TOC1 26 , and PR1, a direct target of NPR1 27 , 3 hours after SA application either in the subjective morning (ZT24) or evening (ZT36) under constant light. Both defence genes had higher induction after the morning treatment (Fig. 4b). We next performed microarray to investigate this time-of-day-specific sensitivity globally (GSE61059) (Patterns of representative genes verified by qPCR shown in Extended Data Fig. 7a, b). We found more genes showing higher induction by morning SA treatment than the evening treatment (like PR1) (Fig. 4c). They are mainly defence-related genes (Fig. 4d). In contrast, a larger number of genes appeared to be more repressed after the evening treatment and they are enriched in plant growth and development (Fig. 4d). Furthermore, promoter analysis of the differentially induced genes in the morning showed significant enrichments for both cis-elements bound by CCA1/LHY and TOC1, and those of the differentially repressed genes in the evening had significantly enriched cis-elements bound by CCA1/LHY (Extended Data Fig. 7c, d). Collectively, these data strongly support our hypothesis that acute perturbation in redox rhythm caused by SA treatment leads to increased expression of both positive and negative regulators of defence, but with the former in the morning and the latter in the evening ( Fig. 3a and Extended Data Fig. 7c, d). This may increase the diurnal differences in sensitivity to pathogen challenge in plants.
Gating defence towards the morning may also be a mechanism for plants to minimize interference on growth at night 28,29 . We observed that SA treatment of dark-grown plants could help sustain a more robust circadian rhythm (P < 10 −15 ) than mock treatment (P < 0.001) (Fig. 4e), but induced a severe loss in fresh weight (Fig. 4f). This is consistent with our hypothesis that untimely induction of immunity at night is detrimental to plant growth. Besides gating immune response, reinforcement of the clock may also help increase photosynthesis and negate the redox perturbation through enhanced expression of CHLOROPHYLL BINDING PROTEIN 2 (CAB2) and evening-phased CATALASE 3 (CAT3) but not the morning-phased CAT2 (Fig. 4g).
We propose that in Arabidopsis, the daily redox rhythm is intrinsically linked with the basal expression of the circadian clock through NPR1 ( Fig. 4h and Extended Data Fig. 8). Perturbation in redox rhythms caused by SA during pathogen challenge is sensed by NPR1 to trigger defence gene expression and to reinforce the circadian clock. The wiring of NPR1 to defence genes as well as to the clock shows how plants gate their immune responses towards the morning to anticipate infection while minimizing fitness costs on plant growth which occurs mainly at night 30 .

Plant materials
The TOC1p:LUC (Col-0), LHYp:LUC (Col-0) 14  , wild-type CAT3, CAT2 and TOC1 promoters and mutated TOC1 promoter (amplified using QuikChange Lighting Multi Site-directed mutagenesis kit, Agilent Technologies) were cloned into the pDONR207 vector (Invitrogen) through the Gateway BP reaction (Invitrogen) and then transferred to the destination vector pGWB235 32 through the Gateway LR recombination reaction (Invitrogen). Agrobacterium-mediated transformation of Arabidopsis was performed as previously described using wild-type plants 33 . Homozygous T3 lines of CAT2p:LUC, CAT3p:LUC and different T1 lines of TOC1p:LUC and TOC1p(TBSm):LUC were selected and used for the luciferase imaging experiment. All primer sequences used for making the transgenic constructs are listed in Extended Data Table 1.

NADP + and NADPH measurement
Three-week-old wild-type (Col-0) plants grown in soil under diurnal condition (12 h light/12 h dark) were treated with water or 1 mM SA at subjective dawn and samples were collected every 4 h for two days under constant light conditions. NADP + and NADPH were measured according to Queval and Noctor 34 with modifications. Briefly, 50 mg of 3-week-old leaves were pulverized in liquid nitrogen using Genogrinder and extracted using 10 mM Tris-HCl (pH 8.0, 1 ml Tris-HCL per 100 mg tissue). The homogenate was centrifuged at 16,000 g for 10 min at 4 °C. The supernatant was separated into two 0.2 ml aliquots. To extract NADP + , 50 μl 1 M HCl was added to one 0.2 ml aliquot. The mixture was heated in boiling water for 1 min. Then 25 μl MES (pH 5.6) was added and the pH of the extract was adjusted to 5-6 using 0.2 M NaOH. To extract NADPH, 50 μl 1 M NaOH was added to the other 0.2 ml aliquot. The mixture was heated in boiling water for 1 min. Then 25 μl MES (pH 5.6) was added and the pH of the extract was adjusted to 7-8 using 0.2 M HCl. Three 20 μl aliquots of the NADP + and the NADPH extracts were used as technical replicates. Samples containing only the extraction buffer were used as blank. The measurement of the samples and the derivation of the standard curves were performed according to Queval and Noctor 34 .

RNA extraction and quantitative PCR
Three-week-old plants grown in soil under diurnal conditions (12 h light/12 h dark) were treated with water or 1 mM SA at subjective dawn and samples were collected every 4 h for two days under constant light conditions. RNA extraction was performed as previously described 18 . cDNA synthesis (SuperScript III, Invitrogen) and quantitative PCR (SYBR Green, Roche) were performed according to the manufacturer's protocols. All primer sequences used for qPCR are listed in Extended Data Table 1.

Luciferase activity measurement
Plants grown in soil with under 12 h light/12 h dark cycles for three weeks were sprayed with 2.5 mM luciferin (Gold Biotechnology) in 0.02% Triton X-100 (Sigma) to activate and deplete pre-existing luciferase due to its instability in the presence of the substrate. Then the plants were transferred into constant light condition. 1 mM SA (Sigma) or water (as control) was sprayed 24 hours later (ZT24). The fifth and sixth leaves were harvested and rinsed three times in 50 ml water. Luciferase was extracted and relative activity was measured according to the manufacturer's protocol (Luciferase Assay System, Promega).

Luciferase imaging
Plants grown in soil under 12 h light/12 h dark cycles for three weeks were sprayed with 2.5 mM luciferin (Gold Biotechnology) in 0.02% Triton X-100 (Sigma) one day before luciferase imaging. Plants were then placed into the imaging system (Nightshade LB985) under either constant light or dark conditions and assayed for bioluminescence by acquiring images with exposure time of 20 min. To test the effect of SA or GSHmee, 1 mM SA (Sigma)/3 mM GSHmee (Sigma) or water (as control) was sprayed at different indicated times. Subsequent quantifications of bioluminescence intensity were performed using Image J.

Analysis of circadian rhythms
The quantified time course bioluminescence data were decomposed into a line and a sine wave with exponentially decaying amplitude: using GraphPad Prism 6. The intercept of the line at the y-axis ("b") was considered as the average expression level. The period and amplitude were inferred from the sine wave. The exponential decay was used to account for the dampening of bioluminescence over time. The best fitted value, standard error and degrees of freedom were used for statistical analysis.

Western blot
Three-week-old wild-type (Col-0) plants grown in soil under diurnal condition (12 h light/12 h dark) were treated with water or 1 mM SA at subjective dawn and samples were collected every 4 h for two days under constant light. Detection of the NPR1 monomer protein on a non-reducing SDS-PAGE gel was performed as previously described using an antibody against NPR1 12 .

Yeast one-hybrid assay
The TOC1 promoter was first cloned into the pDONR P4-P1R vector (Invitrogen) through the Gateway BP reaction. The entry clones were recombined into destination vectors pMW#2 (Invitrogen) and pMW#3 (Invitrogen). Mutagenesis of the TOC1 promoter was performed using QuikChange Lighting Multi Site-directed mutagenesis kit (Agilent Technologies) according to the instruction manual. TOC1p_Mut1 (the TOC1 promoter mutated in the 1 st TGA-binding site), TOC1p_Mut2 (the TOC1 promoter mutated in the 2 nd TGA-binding site) and TOC1p_Mut1+2 (the TOC1 promoter mutated in both TGA-binding sites) were cloned into destination vectors pMW#2 and pMW#3 through the Gateway cloning kit (Invitrogen). The coding sequences of TGAs were cloned into pDONR207 and subsequently transferred into pDEST-AD by the Gateway LR reactions. Transformation of constructs into the yeast strain YM4271 was performed as previously described 35 . Betagalactosidase reporter activities were measured using ONPG as the substrate 36 and normalized to the control with an empty vector pDEST-AD. All primer sequences used for yeast one-hybrid assay (Y1H) are listed in Extended Data Table 1.

Chromatin immunoprecipitation
Three-week-old soil-grown 35S:NPR1-GFP (in npr1-1) 12 plants were treated with either water (CK) or 1 mM SA at dusk and samples were collected 3 hours after treatment. Chromatin immunoprecipitation (ChIP) was performed as described previously 37 . Immunoprecipitation was performed using a polyclonal antibody against GFP (Ab290, Abcam) and Dynabeads Protein G (Invitrogen). The purified ChIP samples were subject to qPCR using primer pairs for the promoter region (−639 to −589 bp upstream of the start codon) and the coding region (+753 to +803 bp downstream of the start codon) of TOC1. Fold of enrichment was calculated using the comparative Ct method 38 using the input samples as normalizers. All primer sequences used for ChIP are listed in Extended Data Table 1.

The mathematical model of the Arabidopsis circadian clock
We applied the P2012 plant circadian model from the Plant Systems Modelling portal 13 to elucidate new connections between SA signalling through NPR1 and known plant circadian genes. This numerical ordinary differential equation (ODE) model in MATLAB consists of 32 ODEs and includes transcription terms for 10 genes, which are LHY/CCA1, PSEUDO-RESPONSE REGULATOR 9, 7, 5 (PRR9, PRR7, PRR5), TOC1, EARLY FLOWERING 4, 3 (ELF4, ELF3), LUX ARRHYTHMO (LUX), GIGANTEA (GI) and ABA receptor (ABAR). The 133 parameters in P2012 were previously fitted to multiple datasets in various lightdark photoperiods, different genetic backgrounds, and ABA signalling. It is important to note that the P2012 model was designed to understand and predict changes in period and phase when perturbed by genetic or environmental variations. The model does not aim to reflect the exact transcriptional profiles or the absolute protein concentrations.
It was recently shown that plant circadian models exhibit a 'period overshoot' when transitioning from LD to LL cycles 39 . This period overshoot introduced a constant phase delay in the LL data relative to the LD data. Since this phenomenon was not observed in the experimental data, it is an artefact of the mathematical model. For example, the P2012 model predicts that TOC1 mRNA peaks at ZT18 (i.e., 6 h after subjective light-to-dark transition) under LL conditions 13 instead of the real peak time at ZT12. Moreover, our experiments indirectly measured TOC1 expression via the luciferase reporter, which is known to exhibit delays 40,41 . This delay was deduced to be 2 h because the luciferase reporter peaked at ZT14 (Fig. 1d) whereas the TOC1 mRNA peaked at ZT12 (Fig. 1c).
To take this 2-h delay and the 'period overshoot' in the model into consideration, we empirically measured a 4-h delay between simulation TOC1 mRNA levels and our luminescence data. This total 4-h delay was inferred by aligning the second peak after the LD to LL transition in our luciferase experiments (38 h) and TOC1 mRNA in the model (42 h). We subsequently used this 4-h delay to correctly align and fit the P2012 simulation TOC1 mRNA to our experimental luciferase data.

Addition of NPR1 regulation to the circadian clock model
While keeping the original P2012 parameters fixed, we added NPR1 as a transcriptional activator of TOC1, as it has been shown experimentally. We also added NPR1 as a transcriptional activator of two additional clock genes ("query pair"). Our goal was to systematically determine which query pair best fitted our measured TOC1p:LUC expression in wild type and npr1-3 in mock-and SA-treated plants. We multiplied the P2012 transcriptional synthesis term of TOC1 and each gene in the query pair by their own NPR1dependent regulatory function F(t) (i.e., non-competitive activation). Each regulatory function F has the form: where [NPR1(t)] is the NPR1 monomer concentration over time, n b is the basal, NPR1independent transcription level of the gene of interest, n a is the maximum NPR1-activated transcription level of the gene and K d is the effective DNA-binding dissociation constant for the gene. The [NPR1(t)] monomer levels for mock-treated and SA-treated plants were taken from western blot data in Fig. 2b and Extended Data Fig. 4. The NPR1 data for the mocktreated and SA-treated plants were then averaged, normalized, and linearly interpolated to serve as an input function for modelling (Extended Data Fig. 9a, b).

Least squares fitting of the TOC1p:LUC data
For every query pair and the TOC1 gene, we optimized n b , n a , K d parameters (among the 9 parameters, 6 are independent; see below) to give the best least squares fit of the TOC1 mRNA in the model to the patterns of TOC1p:LUC expression in wild-type and npr1-3 with mock-or SA-treated plants over several circadian cycles (Fig. 2a).
Because the time of sampling and the waveforms were different between our experiments and the P2012 model, our luciferase data could not be fitted directly to the model. To solve the sampling time discrepancy, data points from our experiments and the P2012 datasets were interpolated (via cubic spline) to a time resolution of 0.1 h. To circumvent the waveform issue, we first calculated the ratio (R) of npr1-3/WT and SA-treated/mock-treated WT in experimental TOC1p:LUC data: where TOC1 WT , TOC1 npr1 and TOC1 SA are interpolated experimental data for TOC1p:LUC in (1) mock-treated wild type, (2) mock-treated npr1 mutant, (3) SA-treated wild type, respectively. We then created target P2012 simulation of TOC1 npr1 or TOC1 SA mRNA data by multiplying the TOC1 WT mRNA data in the simulation by R npr1 or R SA. We optimized the regulatory function parameters (n b , n a , K d ) for each combination of query genes and TOC1. The parameter optimization used nonlinear least squares fitting to minimize the sum of squared residual of TOC1 expression (i.e., squared difference between the model and the target TOC1 mRNA profile). To account for the 4-h delay inherent to 'period-overshoot' in the P2012 model and the use of the reporter, we started fitting at 28 h, which corresponds to 24 h in our experiments. We fitted 3-day-long npr1-3 mutant data and 2-day-long WT and SA treatment data.

Fitting of the npr1-3 mutant data
For the first part of the fitting process, we set [NPR1] = 0, which resulted in a single parameter (n b ) fit for each candidate gene. Because we coupled NPR1-activation to TOC1 and two other query genes, 3 parameters in total were fitted for each query pair. The optimal parameters were restricted to a value between 0 and 1, where 0 represents no transcription in absence of NPR1 and 1 represents the absence of any regulation by NPR1. We used the function fmincon in MATLAB (2013b, MathWorks, Natick, MA, USA) with sequential programming algorithm without restrictions. We found the 3 parameters that minimized the least squares residual of TOC1 model output to target P2012 TOC1 npr1 data. Because nonlinear least squares fitting uses a deterministic algorithm that can become trapped in a local minimum, we ran the simulations from 15 different random starting points to find the global minimum. We confirmed that nonlinear least squares fitting of most query pairs converged to the same global minimum when started from random parameters (Extended Data Fig. 9c, d).
A plot of the best least squares fit for each query gene showed a cross-hair pattern centred on PRR7 in Fig. 3b. The npr1 data (reduced TOC1 amplitude, no change in phase or period) are mostly explained by including PRR7 as an additional target of NPR1 regulation (Extended Data Fig. 6a, b), such that PRR7 levels are reduced in the npr1 mutant (Fig. 3c).
A role for NPR1 in regulating basal expression of PRR7 is consistent with previous genetic data showing that lowering TOC1 expression shortens the circadian period 23 , whereas lowering PRR7 expression lengthens the circadian period 42 . Analytical work has previously shown how mutations in opposing components in the clock can lead to unchanged period 43 . We used our best-fit P2012 model to verify that lowered PRR7 expression lengthens the period by ~2 h, whereas lowered TOC1 expression shortens the period by ~2 h in npr1. Thus, the simultaneous, balanced reduction in expression of two opposing nodes (i.e., TOC1 and PRR7) explains why the period of TOC1p:LUC expression is not altered in npr1.

Mock-treated WT data constrain the parameters
During our fitting procedure for mock-treated WT data, we discovered that optimal always exhibited a simple relationship (Extended Data Fig. 9e). The constraint that explains this empirical relationship is that our final choice of should not alter clock expression in mock-treated WT. Mathematically, should have no effect, on average, on gene expression in mock-treated WT. This condition restricts to depend on because the time-averaged <[NPR1]> in WT is normalized to 1, such that Thus, we recovered the simple empirical relationship observed in Extended Data Fig. 9e as

Fitting of the SA treatment data
By using the best-fit value for from Step 1 and the constraint from Step 2 to fix , we only needed to fit a single parameter (k d ). We restricted k d to lie between 0 (i.e., always maximum expression, insensitive to SA treatment and NPR1 monomer levels) and 5. We verified that allowing k d to be as high as 50 did not significantly improve the least squares fit. For each query pair, we ran nonlinear least squares fit from 15 random starting points. We found that again, in general, they converged to the same global optimum (Extended Data Fig. 9f, g). Fig. 3d shows a cross-hair pattern centred on LHY/CCA1 suggesting that LHY/ CCA1 activation by the induced, arrhythmic NPR1 levels during SA treatment should counteract the effect of TOC1 induction. We noted that the best-fit solution from a TOC1only case shows immediate TOC1 induction after SA treatment (Extended Data Fig. 6i) instead of the observed delay until dusk. However, if we added LHY as an NPR1 target in addition to TOC1, the modelling results best fit our experimental data (Extended Data Fig.  6j). Consistently, LHY expression was found to be induced by SA (Fig. 3e), but reduced in the npr1 mutant (Extended Data Fig. 6g).

Limitations of the model
There are limitations to our model. First, our model was only fitted to the expression of one gene (TOC1p:LUC) under 3 conditions (i.e., wild-type with and without SA treatment, and npr1). Second, our model combines LHY/CCA1 into one gene and cannot resolve the experimental differences that we observed in those genes ( Fig. 3e and Extended Data Fig.  6c, d). Third, our model for the SA-induction data pre-sets PRR7 at maximum expression even without SA treatment (i.e., K d = 0) (Extended Data Fig. 9g). This is unlikely an accurate reflection of a real physiological state. Last, all our experiments and modelling were done under constant light conditions and the additive (acute) light activation terms are effectively 0 because the hypothetical protein responsible for light activation has decayed to 0. These additive light terms should not be affected by our assumption of non-competitive activation by NPR1. Uncovering the proper relationship between the light-dependent terms and the NPR1-dependent terms (i.e., competitive versus non-competitive activation) would require experiments under diurnal conditions. This is outside of the scope of the current paper.
Even though the modelling approach correctly predicted PRR7 and LHY as direct targets of NPR1, it had mixed results with ELF3, which is not a direct target. Our model predicted ELF3 to decrease after SA treatment and to increase in the npr1 mutant. The lack of induction in ELF3 by SA (Extended Data Fig. 6e) was consistent with the model prediction. However, the significantly decreased expression of ELF3 observed in npr1 (Extended Data Fig. 6h) was not in agreement with the model. This discrepancy suggests there are other links between NPR1 and the circadian clock that the current model cannot capture. A future model should be fitted to ELF3p:LUC, which would be an informative constraint.

Bacterial infection
Three-week-old plants grown in soil were pre-treated with water or 1 mM SA at ZT12 and infiltrated 24 hours later with Pseudomonas syringae pv. maculicola (Psm) ES4326 (OD 600nm = 0.001) as previously described 12 . Briefly, 8 plants/genotype/treatment were inoculated with Psm ES4326 and bacterial growth was measured 3 days post inoculation.

Microarray analysis
To test time-of-day-specific sensitivity to SA, 3-week-old soil-grown plants were transferred to constant light condition one day before treatment. Water or 1 mM SA were applied in the subjective morning (ZT24) or evening (ZT36). The fifth and sixth leaves were sampled 0 and 3 hours after treatment. RNA was extracted, amplified, labelled and hybridised to ATH1 GeneChip (Affymetrix) as previously reported 25 . The arrays were normalized with RMA algorithm and centred to median. Two-way ANOVA (P < 0.05) and Student's t-test with multiple comparison correction (P < 0.05, Fold changes > 2) were used to identify genes that were significantly more induced or more repressed by water or SA when treated at ZT24 or at ZT36. Athena program (http://www.bioinformatics2.wsu.edu/Athena) was used to identify cis-elements bound by CCA1/LHY including evening element, CCA1-binding site, CCA1 motif1 BS in CAB1 and CCA1 motif2 BS in CAB1. Hypergeometric distribution was used to determine statistical significance. Enriched gene ontology categories were identified using BiNGO (http://www.psb.ugent.be/cbd/papers/BiNGO/ Home.html).

Fresh/dry weight measurement
Three-week-old soil-grown wild-type (Col-0) plants were transferred into constant dark condition or normal diurnal condition at dusk. After 36 h, water (control) or 1 mM SA was applied. Two days after treatment, pictures were taken and fresh/dry weight was measured.

Statistical analysis
Statistical analysis was performed using GraphPad Prism 6 (GraphPad Software). All the centre values shown in the figures are means of technical (Figs. 1a, 1b, 2d and Extended Data Fig. 1) or biological (all other figures where applicable) replicates. Experiments were repeated twice for Figs. 1a, 1b and CAB2p:LUC in 4g. All other experiments were repeated three times where applicable. Harmonic regression ( ) followed by ANOVA test was used to identify statistically significant oscillation. The null hypothesis is that all data across different time points are sampled from the same normal distribution. Student's t-test with multiple comparison correction was performed to identify statistically significant differences between mock and treated samples. Two-way ANOVA was used to assess significant interactions between genotype and treatment or between time of treatment and treatment. Significant interactions suggest the effect of the treatment is dependent on genotype or time of treatment. Tukey's multiple comparisons test was performed to identify the orders of samples that are significantly different from each other. All statistical tests are two-sided tests where applicable.

Extended Data
Extended Data Figure 1. Circadian oscillation of the NAPDH/NADP + ratio NADPH/NADP + ratios in 3-week-old soil-grown plants derived from Fig, 1a, b. Water (CK) or 1 mM SA was applied at 0 h. Data are mean ± s.e.m. (n = 3). White bars represent subjective days and grey bars represent subjective nights. Harmonic regression analysis suggests significant circadian oscillation of water-treated NADPH/NADP + ratio (P < 0.0001).
Extended Data Figure 2. The effects of exogenous and endogenous SA on TOC1 expression a, Luciferase activity measurements using the TOC1p:LUC plant extracts. Relative luciferase activity of the fifth and sixth leaves from 3-week-old soil-grown TOC1p:LUC plants. Water (CK) or 1 mM SA was applied at ZT24. LL, constant light. a.u., arbitrary unit. Data are mean ± s.e.m. (n = 6 biological replicates; t-test; ***, P < 0.001). b, TOC1p:LUC activity rhythms in 3-week-old soil-grown wild-type (WT) and sid2 plants treated with water (CK) or 1 mM SA at subjective dusk (black arrow) (mean ± s.e.m., n = 8 plants). White bars represent subjective days and grey bars represent subjective nights. The bar graphs represent the estimates of amplitude and average expression of TOC1p:LUC, respectively (mean ± s.e.m.). The letters above the bars indicate statistically significant differences between groups at P < 0.05 (Tukey's multiple comparisons test). n.s., nonsignificant (two-way ANOVA, non-significant interaction between genotype and treatment). This experiment has been repeated three times with similar results. (trx-h3 h5) (mean ± s.e.m., n = 6 plants). LL, constant light. a.u., arbitrary unit. White bars represent subjective days and grey bars represent subjective nights. The bar graphs show the estimates of amplitude and average expression (mean ± s.e.m.; t-test; ****, P < 0.0001). b, TOC1p:LUC activity rhythms in 3-week-old soil-grown wild-type (WT) and npr1 plants treated with water (CK) or 3 mM GSHmee at subjective dusk (black arrow) (mean ± s.e.m., n = 8 plants). The bar graphs represent the estimates of amplitude and average expression of TOC1p:LUC, respectively (mean ± s.e.m.). The letters above the bars indicate statistically significant differences between groups at P < 0.01 (Tukey's multiple comparisons test). **, P < 0.01; ****, P < 0.0001 (two-way ANOVA). These experiments have been repeated three times with similar results.
Yellow represents a high value or a target of CCA1/LHY or TOC1. Blue represents low a value or not a target of CCA1/LHY or TOC1. 'X' represents a gene that was more induced by SA when treated at ZT24 (c) or more repressed by SA when treated at ZT36 (d). Arrows represent activation. Block arrows represent repression. P values were determined based on hypergeometric distribution. Extended Data Figure 8. NPR1 senses and transduces redox signals to trigger transcriptional reprogramming SA-triggered redox changes induce the oligomer-to-monomer switch of NPR1. The monomer then enters the nucleus and up-regulates both defence genes and clock genes through interaction with TGA TFs. Extended Data Figure 9. Technical details for model fitting a, Normalized NPR1 monomer abundance in mock-treated samples. The blue line presents the mean values from Fig. 2b, where the value at 48 h (marked with an open star) was inferred to be the same as 0 h. The red line represents the smoothened values used for modelling by averaging over two days to create a one-day trace, which was then repeated over two days. The smoothened data were normalized, such that the time average of NPR1 was equal to 1. LL, constant light. a.u., arbitrary unit. White bars represent subjective days and grey bars represent subjective nights. b, SA-treated NPR1 monomer abundance. NPR1 monomer abundance after SA treatment from Extended Data Fig. 4b was normalized so that 0 h has the same value as the corresponding mock-treated NPR1 monomer level. Based on the assumption that the SA induction lasted for two days, the value of the last time point was inferred to be equal to the basal level (marked with an open star). c, Coefficient of variation (CV) of least squares residual Σ for 15 different, random initial parameters for the model fitting of npr1 data. d, CV of n b * for 15 different, random initial parameters for the model fitting of npr1 data. e, Optimal exhibit a linear relationship. log(Σ) was plotted as a function of n a and K d for mock-treated TOC1-only coupling (no query pairs). A "low", linear Σ region is evident and is described by a simple analytical linear relationship. n b * =0.5689 h −1 . f, CV of Σ for 15 different, random initial parameters for the model fitting of SA-treated data. g, The CV of K d * for 15 different, random initial parameters for the model fitting of SA-treated data.