Circadian protein regulation in the green lineage II. The clock gene circuit controls a phospho-dawn in Arabidopsis thaliana

The 24-hour rhythmicity in the levels of many eukaryotic mRNAs contrasts with the long half-lives of most detected proteins. Such stable molecular species are not expected to reflect the RNA rhythms, emphasizing the need to test the circadian regulation of protein abundance and modification. Here we present circadian proteomic and comparable phosphoproteomic time-series from Arabidopsis thaliana plants, estimating that 0.4% of quantified proteins and a much larger proportion of quantified phospho-sites were rhythmic under constant light conditions. Approximately half of the quantified phospho-sites were most phosphorylated at subjective dawn, when SnRK1 protein kinase was a candidate regulator. A CCA1-over-expressing line that disables the plant clock gene circuit lacked most or all circadian protein phosphorylation, indicating that the phospho-dawn is regulated by the canonical clock mechanism. The few phospho-sites that fluctuated despite CCA1-over-expression still tended to peak in abundance close to subjective dawn, indicating that their temporal regulation was less dependent on the clock gene circuit. To test the potential functional relevance of our datasets, we conduct phosphomimetic experiments using the bifunctional enzyme fructose-6-phosphate-2-kinase / phosphatase (F2KP), as an example. The clock gene circuit controls diverse protein targets in metabolism and physiology via phosphorylation, including where the bulk abundance of the cognate proteins is arrhythmic.


Introduction
Most circadian observations in Arabidopsis can be explained by a genetic repressilator which is composed of a network of mostly negatively interacting transcription factors (Pokhilko et al, 2012;Fogelmark & Troein, 2014). In addition to transcriptional interactions, this transcriptional-translational feedback loop (TTFL) system requires post-translational modification of transcription factor proteins . Phosphorylation of CCA1 by casein kinase (CK) 2 for example is necessary for the circadian clock function (Daniel et al, 2004).
Protein phosphorylation is involved in the circadian clock mechanism not only in plants but also in fungi, animals and cyanobacteria (Nishiwaki et al, 2004;. While the transcription factors of TTFLs in animals, fungi and plants are not evolutionarily conserved , many kinases that play a role in circadian timekeeping across eukaryotes, are similar. For instance CK2 is also important for circadian timing in mammals (Maier et al, 2009) and fungi (Yang et al, 2003). Protein phosphorylation is also involved in the output of the circadian clock (Kusakina & Dodd, 2012). However, only one study has so far addressed the question of how pervasive circadian protein phosphorylation is (Choudhary et al, 2015).
In addition to TTFLs, circadian oscillations can be driven by mechanisms that are independent of rhythmic transcription, non-transcriptional oscillators (NTOs). The cyanobacterial circadian clock is based on rhythmic autophosphorylation of the KaiC protein together with the KaiA and KaiB protein and this mechanism does not even require a living cell to create oscillations (Nakajima et al, 2005).
Evidence for NTOs also exists in eukaryotes.  found that in Ostreococcus tauri cells, which do not transcribe in extended darkness, the protein peroxiredoxin (PRX) is rhythmically overoxidized during 3 days of constant darkness. The same overoxidation rhythm was found in red blood cells which are enucleated and therefore do not possess a TTFL . Rhythmic PRX over-oxidation was also found in transcriptionally arrhythmic Arabidopsis TOC1-Ox plants as well as a Neurospora line that lacks a functional TTFL. The PRX over-oxidation rhythm even exists in cyanobacteria and archaea. Feeney et al, 2016 show that circadian magnesium ion transport can occur in transcriptionally inactive Ostreococcus tauri cells. Therefore, at least some eukaroytes possess an NTO which appears to have the same redox modification output and therefore may be evolutionarily ancient and conserved (Edgar et al, 2012).
In this study, we used mass spectrometry based proteomics and phosphoproteomics on circadian time courses to address the following questions: (1) How pervasive are rhythms in protein abundance and phosphorylation as a clock output in a normally functioning circadian clock system? and (2) can protein abundance or phosphorylation be rhythmic in a plant with a disabled transcriptional oscillator?
To investigate (1), we used a time course of WT plants, and for addressing (2), we generated time courses of CCA1 overexpressing plants (CCA1-Ox), which are impaired in their TTFL. We generated global proteomics and phosphoproteomics data in parallel from the same protein extracts. Our analysis revealed that the transcriptional oscillator is required for most if not all rhythmic protein phosphorylation, and that most rhythmic phosphopeptides peak at subjective dawn. We also found this 'phospho-dawn' trend in the time courses with fluctuating phosphopeptides in the CCA1-Ox. Lastly, we selected one rhythmically changing phosphosite of the bifunctional enzyme fructose-6-phosphate-2kinase / phosphatase (F2KP) and demonstrate that phospho mimetic mutation of the rhythmic site can have functional relevance.

Sample preparation
Protein extraction and precipitation was carried out according to method 'IGEPAL-TCA' in our previous publication (Krahmer et al, 2015). Briefly, protein was extracted and precipitated with TCA and phase separation, then washed with methanol and acetone. 500μg re-suspended protein was digested using a standard in-solution protocol and peptides were desalted. Before drying, the eluted peptides were separated into two parts: 490μg of the digest was used for phosphopeptide enrichment, with decoy search and an ion-cutoff of 20. In all but one cases, these parameters resulted in a false-discovery rate (FDR), of less than 5% (phospho I: 3.5%, phospho II: 3.2, global I: 6.8%, global II: 4.5%, calculated using the formula 2*d/(n+d) (Elias & Gygi, 2010), n and d being the number of hits in the normal and decoy databases, respectively, using an ion score cutoff of 20). Peptides were quantified by their peak area by Progenesis, and proteins were quantified by using the sum of the quantitation of all their unique peptides. Where peptides matched very similar proteins, multiple accession numbers are shown in exported results from Progenesis (supplementary data S1, S2).
In order to remove duplicates of phosphopeptides due to alternative modifications other than phosphorylation or missed cleavages, we used the qpMerge software following the Progenesis analysis (Hindle et al, 2016). The data are publicly available in the pep2pro database (Baerenfaller et al, 2011) at www.pep2pro.ethz.ch (Assembly names 'ed.ac.uk Global I', 'ed.ac.uk Global II', 'ed.ac.uk Phospho I', 'ed.ac.uk Phospho II') and have been deposited to the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE partner repository (Vizcaíno et al, 2014) with the dataset identifier PXD009230. Exported .csv files from Progenesis with all peptide and protein quantifications can be found in the online supplementary material (Data S1 and S2).

General statistics and outlier removal
Results of statistical analyses are summarised in data S3. Outlier analysis, statistics and Venn diagrams were done using R version 3.2.1 (R Core Team, 2015). Zero values for the quantification were exchanged for 'NA'. For outlier analysis and parametric tests such as ANOVA, arcsinh transformed data was used to obtain a normal distribution, while untransformed data was used for plotting time courses and for the non-parametric JTK_CYCLE analysis. For phosphoproteomics analysis all replicates in which the Pearson correlation coefficient among replicates of the same time point was lower than 0.8 were regarded as outliers (supplementary methods, supplementary data S4). In global dataset II, the first run replicate of each time point had to be excluded as an outlier due to an apparent drift (figure S13, S15). To generate heat maps, the abundance of each peptide or protein was normalized by the time course mean of the peptide or protein, followed by taking the log2 to centre values around 0, and the heatmap.2 function from the pvclust R package was applied (Suzuki & Shimodaira, 2015).

JTK_CYCLE analysis
We used the R-based tool JTK_CYCLE (Hughes et al, 2010) to determine rhythmicity, with the following modifications: (1) we ran the JTK_CYCLE algorithm for each phosphosite or protein separately rather than the entire list, to allow handling of missing quantification values for some replicates. Benjamini Hochberg (BH) (Benjamini & Hochberg, 1995) correction of p-values was carried out after application of the JTK_CYCLE tool.
(2) Since our time course durations are close to the periods of rhythms we are searching for, some identifications were assessed as rhythmic by the original JTK_CYCLE tool that were increasing or decreasing continually over the entire time course.
We excluded these from the group of rhythmic identifications. These continually rising or falling identifications were marked with 'excl.' in front of the p-value (supplementary data S5).
Kinase prediction using GPS3.0 Kinases for each site were predicted using the Arabidopsis thaliana specific GPS 3.0 prediction tool (Liu et al, 2015). For each phosphorylation site, an amino acid sequence was generated that contained 50 amino acids on either side of the phosphorylated residue using a python script. This resulted in 101 amino acid long sequences, unless the phosphosite was closer than 50 residues to the C or T terminus in which case the missing positions were filled by 'X'. The phosphopeptides used as foreground were the significantly rhythmic phosphosites (JTK_CYCLE p-value < 0.05) peaking at a given time point, all other phosphosites identified in the same experiment were used as background. The high threshold setting was used to minimize false-positive predictions and searches were done for S and T residues. In order to reduce the complexity of the dataset, we used a simplification: where kinases from different families were predicted, only the one with the highest difference between score and cutoff was used.
For foreground and background the numbers of predictions for each occurring kinase group were counted and the Fisher's Exact test was used to determine predicted kinases that were significantly enriched in the foreground group (p<0.05).

GO analysis
For GO analysis, foreground and background were chosen as in the kinase prediction analysis. With these groups, GO analysis was conducted using the topGO script (Alexa et al, 2006), followed by Fisher's exact test to determine enrichment of terms (supplementary data S6).

Generation of F2KP point mutations and expression constructs
The F2KP coding sequence in the pDONR221 vector, lacking a stop codon, was kindly provided by Dr Sebastian Streb, ETH Zürich. The QuikChange Lightning Site-Directed Mutagenesis kit (Agilent Technologies) was used to introduce point mutations using primer pair AspF and AspR for mutation of S276 to aspartic acid or AlaF and AlaR for mutation to alanine (Table S 3). The F2KP coding sequence (the mutated constructs as well as the WT construct) was amplified by PCR using primers F2KP-F and F2KP-R (Table S 3), introducing restriction sites for AflII (3' end) and XbaI (5' end) and a stop codon.
Using these restriction sites, digested PCR products were ligated into the pmcnEAVNG expression vector which allows expression in E.coli with an N-terminal GST tag and a T7 promoter for IPTG inducibility. Plasmids were transformed into Rosetta™2(DE3)pLysS Competent expression strain (Novagen).

Expression of GST-F2KP constructs in E. coli cells
The three constructs, WT, S276E and S276A were expressed in E.coli and purified using the GST tag.
Two independent expression experiments were performed (experiment 1 and experiment 2), each in triplicates. 200ml E.coli cultures were grown at 37°C with 100μg/ml ampicillin and 25μg/ml chloramphenicol and induced with 1mM IPTG at OD600 values between 0.6 and 0.8. Cells were harvested by centrifugation after 2.5h of expression at 37°C. Each pellet was from 50ml bacterial culture. For purification of GST-F2KP, pellets were lysed in 2.5ml PBS with Complete protease inhibitor cocktail (Roche) with a probe sonicator. After clearing of the lysates, AP was carried out using GSH-agarose beads, with 167µl GSH-agarose bead suspension (Protino glutathione agarose 4B, Macherey-Nagel), a binding incubation of 30min at room temperature, 4 washes with 10x the volume of the bead suspension and elution in PBS with 100mM reduced glutatione (pH 8.0) 3 times 30min at room temperature.

F2KP activity assay
F2KP activity of purified F2KP was measured as described in Kalt-Torres & Huber, 1987 (F-2,6-BP producing reaction) and Van Schaftingen et al, 1982 (measurement of generated F-2,6-BP by its activation of PFP and subsequent production of NADH produced from glycolytic enzymes).

Western blot quantification of F2KP concentration in AP eluates
To test whether the differences in F2KP activity of eluates could be caused by differences in abundance in the eluate, we quantified the amount of F2KP in equal volumes of eluates by western blotting..
Samples were prepared for SDS-PAGE with 25% 4xLDS (Life Technologies, NP0008) and 20mM DTT and were incubated at 70°C for 10min. Two concentration series of an equal mic of all three eluates were loaded on each gel for relative quantitation. Individual eluates were examined in duplicates (expression 1) or triplicates (expression 2). 4-12% Bis-Tris Gels (Life Technologies) were run and protein was blotted onto nitrocellulose membrane using the iBlot ® system. Two primary antibodies were used in parallel: anti-F2KP from rabbit, raised against amino acids 566-651 (Villadsen & Nielsen, 2001) and anti-GST from mouse (Thermo Scientific), both at a dilution of 1:1,000 overnight. Secondary antibodies were goat-anti-rabbit (IRDye®800CW, LI-COR) and goat-anti-mouse (IRDye®680RD, LI-COR). Bands were quantified using the ImageStudioLite (version 2) software.

Circadian protein phosphorylation requires the canonical transcriptional oscillator
We generated global proteomics and phosphoproteomics datasets for two independent circadian time courses and for two genotypes each -WT and the CCA1-Ox line which has an impaired circadian oscillator (Wang & Tobin, 1998). This resulted in four datasets: dataset I (Zeitgeber time (ZT) 12 to ZT32) and dataset II (ZT24 to ZT48 (CCA1-Ox) or ZT52 (WT)) ( Figure 1). We removed outliers before conducting further stastical analysis as detailed in the methods section.
We identified 2287 phosphopeptides in dataset I, 1164 in dataset II, which condensed to between 1000 and 1500 in each dataset after applying qpMerge (Hindle et al, 2016) to remove duplicate phosphopeptides (Table 1 a). These were from several hundred proteins in each dataset and over 1000 in both datasets together (Table 1b). In the global protein analysis, we identified 1896 and 1340 proteins in dataset I and II, respectively, adding up to a total of 2501 for both datasets combined (Table   1a,b). To assess the circadian rhythmicity of each phosphopeptide or protein, we employed the nonparametric JTK_CYCLE tool (Hughes et al, 2010) as it can be applied to time courses of only one cycle, taking the curve shape into account. Unless otherwise stated we considered periods of 22 to 26h and excluded continuously increasing or decreasing profiles from the group of rhythmic phosphopeptides or proteins. In dataset I, which started 12h earlier than dataset II, 606 (40%) phosphopeptides were rhythmic in the WT, in dataset II 100 (8.8%) when using p-values; 338 (23%) (dataset I) and 26 (2.3%) (dataset II) after carrying out p-value correction for multiple testing using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995). The fraction of rhythmic proteins in the global proteomics analysis was smaller than in the case of phosphopeptides: 171 (9.0%) in dataset I, 45 (3.4%) in dataset II had JTK_CYCLE p-values < 0.05, 6 (0.32%) and (0.45%) with q-value < 0.05.
Much fewer phosphopeptides and a similar number of global proteins were rhythmic in the CCA1-Ox, however after adjusting for multiple testing, no significant identifications remained apart from two phosphopeptides in dataset I (Table 1a).
As expected, the global analysis did not cover all proteins identified by phosphoproteomics ( Figure   2a,c). We only found very few accessions with rhythms in protein as well as rhythmic phosphopeptides. About half of the rhythmic phosphopeptides or proteins also had rhythmic transcripts 'Phospho-dawn': most rhythmic phosphopeptides peak in the subjective morning Analysis of the number of phosphopeptides that peak at each time point revealed that 45% (dataset I) and 73% (dataset II) of rhythmic phosphopeptides peak at subjective dawn in the WT (Figure 3, Figure   4a-b, Figure S 1a -b). This is in agreement with previous observations in Ostreococcus and Arabidopsis (Noordally et al, 2018;Choudhary et al, 2015). This trend to peak at dawn was not observed in the global proteomics dataset ( Since phospho dawn may not be completely abolished by impairment of the TTFL, we reasoned that there may be a kinase activity that is active in the subjective morning which is either very robustly dawn-timed, even when the TTFL is barely operating, or an alternative oscillator, possibly an NTO, contributes to the dawn phased kinase activity. Therefore, characterisation of these dawn-phased phosphopeptides could help to identify either a very robust dawn-phase TTFL output, or potentially even consequences of an NTO. We therefore focused on this group of phosphopeptides for further data analysis: GO analysis to search for common enriched functions of these phosphoproteins ; kinase prediction enrichment to obtain candidate kinases that may be involved in phospho dawn; and phosphosite motif analysis. In all these analyses we used the ZT24 or ZT48 peaking rhythmic (JTK_CYCLE p-value <0.05) phosphopeptides as foreground and all other identified phosphopeptides as background. As we used p-values instead of q-values in these analyses, we regard only those results that are supported by both datasets as reliable.

GO analysis
We analysed GO term enrichment on the phospho dawn group, and other time points where the foreground was sufficiently large (Table 3a,b, supplementary data S3). Only one GO term ('cotyledon development') was shared between the two phosphoproteomics datasets and the ZTs of enrichment are 16h apart (Table 3b). Several terms were shared between the WT and the CCA1-Ox (allowing shorter periods) in dataset I, most of them related to energy metabolism or ion homeostasis, and all of them were enriched at ZT24 or ZT28. Apart from overlap of exact GO IDs, we found enrichment of terms related to translation in the WT at 24h in both phosphoproteomics datasets, which is consistent with rhythmic phosphorylation of RPS6 (supplementary data S3, Choudhary et al, 2015)).

Kinase prediction and phosphorylation motif analysis
For predicting candidate kinases phosphorylating the subjective-dawn phased phosphopeptides, we used the software GPS3.0 (Liu et al, 2015) (Xue et al, 2008) because it offers a species-specific setting for Arabidopsis thaliana. We searched for enrichment of kinase groups among the rhythmic (JTK_CYCLE p-value <0.05,  (Hrabak et al, 2003). Among the phospho dawn peptides we found known (direct or indirect) SnRK1 targets ( Figure 5, Table S 2). To complement the GPS3.0 analysis, we searched for the SnRK1 target pattern, firstly with hydrophobic amino acids at the -5 and +4 positions  (supplementary methods) and found significant enrichment in the phosphopeptides peaking at ZT24 in dataset II, allowing periods between 22 and 26h in the CCA1-Ox (Table S 4). However when including the requirement for a basic amino acid at position -4 or -3 (Sugden et al, 1999) the number of sequences identified in foreground and the background was too low for meaningful statistical analysis.
We generated target site motifs using the pLogo tool (O'Shea et al, 2013) (supplementary methods, Figure S2, Figure S3). The only significantly enriched position in the WT was D at the +1 (dataset I).
In the CCA1-Ox, a D was significantly enriched at the -5 position in dataset I, as well as Ms at the +1 or -6 positions and a V at -2 in dataset II. However, none of these over-representations were consistent between datasets (Figure S 2, Figure S 3).

Rhythmically phosphorylated kinases and phosphatases
Since kinases and phosphatases themselves can be regulated by phosphorylation, we were interested in rhythmic phosphopeptides of kinases and phosphatases. Identification of rhythmic kinase or phosphatase activities in the WT could therefore help to discover new clock output pathways via protein phosphorylation, as a consequence of normal clock function, or components of an NTO.
For two kinases, CRK8 and AT5G61560, we found rhythmic phosphopeptides in the WT with nonchanging protein abundance, indicating that rhythmicity is due to phosphorylation rather than changes in protein abundance (Figure 6a,b). CRK8 is a member of the CDPK-SnRK1 superfamily (Hrabak et al, 2003).
All rhythmically phosphorylated phosphatases in our data are members of the protein phosphatase 2 C (PP2C) family ( Figure 7) and they are rhythmic only in the WT. Some PP2Cs can dephosphorylate and inactivate SnRK1 (Rodrigues et al, 2013) and may therefore also play a role in phospho dawn.
For most of the kinases and phosphatases shown here, the CCA1-Ox vaguely follows the WT pattern, but it is not classified as rhythmic by JTK_CYCLE (with the exception of At5g14720, Figure 6d), therefore these phosphorylation events are likely to be an output of the TTFL.

Functions of rhythmic phosphoproteins in the WT include transport, regulation of translation, photosynthesis, carbon and nitrogen metabolism
Phosphoproteins that were rhythmic in the WT are associated with a variety of physiological functions (  (Schepetilnikov et al, 2013;Jackson et al, 2010) (Choudhary et al, 2015). CINV1 is required for carbon metabolism and growth (Barratt et al, 2009), SWEET12 facilitates sucrose transport (Chen et al, 2012).

Only few proteins are rhythmically phosphorylated in the CCA1-Ox in both datasets
The overlap of rhythmic phosphopeptides between datasets for the CCA1-Ox is small (Table 1, Table   2). With 22 to 26h periods, phosphopeptides of a Leucine rich repeat protein kinase as well as Dormany/auxin associated family protein are shared, though the latter also showed a significant increase in its protein abundance in dataset I (supplementary data S2) in the CCA1-Ox, therefore changes may be due to increasing protein expression in constant light. When allowing shorter periods, NIA1 as well as a remorin family protein were rhythmic.

Circadian rhythms in protein abundance in the global proteomics datasets
We applied the same rhythmicity, peak time distribution and GO analyses to the global proteomics data (Table 3c,d, supplementary data 3). In both datasets, enriched GO terms in the WT at the end of the subjective day (ZT12 and ZT36) were photosynthesis related terms, 'response to glucose', 'regulation of protein dephosphorylation' and an oxidoreductase activity with NAD or NADP as acceptor. The latter was also enriched in the CCA1-Ox in dataset I (Table 3d). In the CCA1-Ox, two terms related to cell wall metabolic processes were enriched in both datasets during the subjective night (ZT20 and ZT44) (Table 3d). Looking at individual rhythmic proteins shared between the datasets, we identified six for he WT, three for the CCA1-Ox (Table 2d,e). One protein, INOSITOL 3-PHOSPHATE SYNTHASE 1 (MIPS1, AT4G39800) was rhythmic in both datasets and both genotypes.
Among the non-shared rhythmic proteins in the WT were several ribosomal proteins and proteins of the photosynthetic machinery (supplementary data 2), in agreement with published data (Choudhary et al, 2015).

F2KP is rhythmically phosphorylated at the conserved residue Ser276
As an example of how our datasets can be used to investigate new clock output pathways, we analysed the molecular function of a rhythmic phosphosite in the bifunctional enzyme F2KP. Several phosphosites were detected in F2KP with only S276 showing a circadian rhythm in the WT in dataset II (Figure 5d). Only one F2KP peptide was detected in the global protein analysis and only in dataset II. Since its changes over time are not significant and not in phase with the phosphopeptide including S276 it is unlikely that the rhythm in S276 phosphorylation is caused by changes in F2KP protein abundance.

Point mutation of S276 affects F2KP kinase activity in vitro
To test if the S276 phosphorylation site is functionally relevant, we measured the maximum F6P,2K activity of phosphomimetic mutants S276D and S276A with the unmutated WT control enzyme in vitro. Two independent expressions and subsequent purifications using a GST tag were carried out to ensure reproducibility. S276A had an approximately 2.5 fold increased activity compared to the unmutated version, while S276D had only slightly increased activity, which was statistically significant in one out of two independent experiments (Figure 9a,b, Figure S5a,b,d). Equivalent amount of expressed F2KP protein in assays was verified by western blotting (Figure 9, Figure S5) with two different antibodies (Figure 9c, Figure S5c).

Discussion
In this study, we generated and analysed global and phospho-proteomics time courses for two mains reasons: Firstly, by studying WT plants in constant conditions, we aimed to find out how pervasive rhythmicity of protein phosphorylation and abundance is in constant conditions and a normally functioning clock; this also involved identification of rhythmic proteins / phosphorylation events and could allow discovery of novel clock output pathways through protein phosphorylation. Secondly, by applying the same approach to the CCA1-Ox line, we aimed to identify candidate processes, proteins and phosphorylation events that are rhythmic in spite of an impaired transcriptional oscillator, and may therefore be driven by an NTO.

Rhythmic protein phosphorylation is more wide-spread than rhythmic protein abundance
According to our data (Table 1), the proportion of rhythmic phosphopeptides is larger than for proteins, but smaller than for rhythmic transcripts (Covington et al, 2008). The fact that the proportion of rhythmic identifications is higher in dataset I than in dataset II is likely due to dampening of rhythms in the later sampled dataset II. In the CCA1-Ox, some rhythmic phosphopeptides were detected, however amplitudes were lower than in the WT and most changing phosphopeptides had a shorter, and therefore not circadian period.
We found rhythmic phosphorylation of proteins related to a variety of processes in the WT, many of which are in agreement with a previous study (Choudhary et al, 2015), such as proteins involved in nitrogen metabolism and translation. GO analysis of our phosphoproteomics data also suggests circadian rhythmicity of ion transport or homeostasis related processes in the WT at subjective dawn, in agreement with rhythms of calcium ion concentration in Arabidopsis (Xu et al, 2007) and magnesium ions in Ostreococcus and other organisms (Feeney et al, 2016). GO analysis for the global protein data confirmed enrichment of photosynthesis related terms (Choudhary et al, 2016). Of particular interest is also the enrichment of oxidoreductase activity with NAD or NADP as acceptors, and protein dephosphorylation processes at subjective dusk.

Identification of candidate protein kinases involved in phospho-dawn
Our kinase prediction and enrichment analysis revealed enrichment of some CMGC subgroups, such as MAPK, CK2, GSK, DYRK, CDK or DAPK. CK2 is involved in the circadian clock function in Arabidopsis by phosphorylating CCA1 (Sugano et al, 1999;Daniel et al, 2004). A previous study reported enrichment of predominantly CK2 predictions among significantly changing phosphopeptides, though usage of the whole genome as enrichment background may have introduced a bias for highly abundant kinases and targets (Choudhary et al, 2015). Roles for MAPK and GSK have been reported for the circadian clock function in other eukaryotes (Akashi & Nishida, 2000;Iitaka et al, 2005;Besing et al, 2015;Martinek et al, 2001). Our kinase prediction analysis focused on the phospho-dawn peptides where the most consistent enrichment -at both dawn time points, in both datasets, and both genotypes -was for CAMK groups. On the database EKPD, which GPS3.0 is based on, the CAMK group in plants comprises 89 members of the CDPK-SnRK superfamily of kinases.

SnRK1 as a candidate protein kinase for phospho-dawn
Even though all of these 89 CDPK-SnRK member are potential candidates, SnRK1 appeared the most plausible one for three main reasons. Firstly, relevance of SnRK1 in circadian timing has previously been reported: SnRK1 genetically interacts with the clock protein TIC (Sánchez-Villarreal et al, 2018), affects circadian period in the light (Shin et al, 2017), and is involved in entrainment to sugar (Frank et al, 2018). Secondly, SnRK1 is an important metabolic hub, and in normal light-dark conditions the morning is associated with profound metabolic changes in plants, such as transition from using starch to direct photoassimilates, but also transition from normal starch usage to a starvation response if light intensities are low in the morning but starch is alomost depleted. Thirdly, among the phospho dawn phosphosites are known, rather than just predicted SnRK1 targets and interactors. Targets include F2KP, sucrose phosphate synthase and nitrate reductases NIA1 and NIA2 ( Figure 5). In addition, PP2Cs which may regulate SnRK1 have phosphosites that peak at subjective dawn (Woods et al, 2003;Rodrigues et al, 2013). DUF581 interacts with SnRK1 and serves as platform for SnRK1 signalling (Nietzsche et al, 2014). DUF581 protein abundance was rhythmic in global dataset I (not detected in dataset II), one phosphopeptide was rhythmic in the WT in both datasets and in the CCA1-Ox in dataset I (Table 2 and supplementary data S2).
Interestingly, SnRK1 signaling is regarded as antagonistic to TOR signaling (Nukarinen et al, 2016) but both RPS6A and RPS6B, ribosomal proteins and TOR targets, are rhythmically phosphorylated in the WT, and in dataset I also in the CCA1-Ox (supplementary data S2), with their phosphorylation peaks at subjective dawn, confirming published data (Choudhary et al, 2015). This indicates that the interplay between SnRK1 and TOR may be more complex than simply antagonistic.

Could an NTO contribute to phospho-dawn?
Since the Kai proteins are not found in eurkaryotes, they cannot be the drivers of eukaryotic clocks.
However, the fact that KaiC autophosphorylation creates 24h rhythms in a test tubes shows that it is biochemically possible phosphorylation reactions to drive circadian rhytms. In addition, since kinases exist that are involved in clocks across kingdoms, phosphorylation in circadian timing may be more conserved than the actual clock genes. The conserved kinases are therefore candidates for NTO components or links between the NTO and the TTFL. For these reason, we hypothesized that protein phosphorylation could be involved in a plant NTO mechanism. We used the CCA1-Ox as an approximation of a plant without a functional canonical TTFL. We cannot fully exclude a remaining low amplitude canonical TTFL, or an alternative transcriptional oscillator operation in the CCA1-Ox since transcriptomics studies in this line have not been published. However, simulation of increased CCA1 levels with the Arabidopsis clock model predicts period lengthening (data not shown), while with the relaxed JTK_CYCLE p-value < 0.05 (rather than q-value) we observe either approximately 24h or shorter periods in this line. Therefore, there may be room for speculations that either the TTFL has an extremely robust metabolic or protein kinase activity output at subjective dawn even if otherwise severely impaired, or that an NTO gives a time cue that has a small contribution to phospho dawn. It did not escape our attention that indeed some properties of SnRK1, or related protein kinases, would make them good candidates for NTO components, given likely conservation of the NTO and parallels to the cyanobacterial NTO. For example, SnRK1 is conserved with animals and fungi (Polge & Thomas, 2007;Hrabak et al, 2003), with AMPK being involved in circadian timing in mammals (Lamia et al, 2009). In addition, the activating kinases of SnRK1, SnAK1 and SnAK2 (Crozet et al, 2010) and other members of the family auto-phosphorylate, some of them very slowly, like KaiC (Oh et al, 2012;Iwasaki et al, 2002;Nemoto et al, 2015). Like AMPK (Oakhill et al, 2011;Gowans et al, 2013) and SNF-1 (Mayer et al, 2011), SnRK1 is regulated by adenylates as it is protected by AMP from dephosphorylation by PP2C in vitro (Sugden et al, 1999) and interacts with adenosine kinase (Mohannath et al, 2014). Adenylates have previously been suggested as NTO components in an alga (Goto et al, 1985). In summary, even though the preservation of the phospho-dawn and its predicted kinases in the CCA1-Ox and WT are intriguing, we did not find overwhelming evidence for genuinely rhythmic phosphopeptides and the existence of an NTO in the CCA1-Ox.

Phosphomimetic studies illustrate how datasets can be used to discover novel clock output pathways
The specific roles of most of the rhythmic phosphorylation sites identified in our study have not been investigated. To link circadian phosphorylation of a protein to its function, we analysed the effect of a phosphosite of the enzyme F2KP on its activity using phosphomimetic mutations in vitro. F2KP is one of the regulators of carbon partitioning into starch and sucrose  and is necessary to maintain normal growth in fluctuating light conditions (McCormick & Kruger, 2015).
The phosphosite of interest, S276, is within the plant-specific regulatory N-terminal domain (Villadsen et al, 2000;Villadsen & Nielsen, 2001) but is not among the phosphosites in the 14-3-3 binding site ( Kulma et al, 2004). S276 is also a likely SnRK1 target (Cho et al, 2016), even though it lacks a typical SnRK1 target pattern. S276 is conserved in many plant species or replaced by a threonine residue ( Figure S 4). Our in-vitro phosphomimetic experiments suggest that a lack of negative charge at position 276 increases F2KP's kinase activity. F-2,6-BP levels in the plant increase slowly across the day in short day conditions (Kulma et al, 2004). Since F2KP is the only enzyme that converts between F-2,6-BP and F-6-P, F2KP activity must gradually increase during the day. Our results suggest that S276 phosphorylation could be responsible as S276 phosphorylation decreases during the day, as out measurements suggest, increases its activity. This would favour channeling of carbon into starch more at the end of the day than in the morning. Since starch builds up at a fairly constant rate across the day (Smith et al, 2004), this does not appear to translate directly into increasing carbon fluxes into starch.
Rather, we hypothesize that F2KP regulation may play together with other regulators of carbon partitioning and may be more relevant in fluctuating than constant light, where it is important for growth (McCormick & Kruger, 2015), for example to guarantee sufficient carbon supply for sucrose formation in the morning but prioritizing building up of starch resources toward the end of the day.
In conclusion, we have demonstrated how these phosphoproteomics datasets can be used to explore new circadian clock output pathways.          Tables   Table 1:

Key numbers for global and phosphoproteomics datasets. a) Numbers of quantified and changing identifications in each dataset, b) shared and added numbers of both datasets. In b)
protein IDs of peptides are used since peptide IDs are not comparable between experiments. Table 2: List of all proteins and phosphopeptides that are rhythmic in WT or CCA1-Ox in both datasets (i.e. shared rhythmic IDs). All p-values are based on JTK_CYCLE analysis. * in some phosphopeptides, the location of phosphorylated residues differs slightly (either shifted by 1-2 residues or one of several phosphates is missing) between datasets, in which case the phosphorylated residues of both datasets are noted (dataset I / dataset II). Bold: phase different of peak is up to 4h. fam. = family, pr. = protein Table 2 continued              Table S 1: GPS3 kinase prediction, using ANOVA p-values. ANOVA analysis on each phosphosite time course was carried out on the data after outlier exclusion, treating 0 values as missing values. GPS3 predictions from Phosphopeptides with ANOVA p-values less than 0.05 were used as foreground for testing enrichment of kinases predictions in all significantly changing phosphopeptides, while phosphopeptides with p<0.05 and a given peak time were used to analyse for enrichment at specific times of the day (24 or 48h). In the case of several predicted kinases within a family of kinases for the same phosphosite, only the one with the highest difference of score and cutoff was considered. Predictions for each kinase were counted in the foreground as well as the background and a Fisher's exact test was used to test for significant enrichment. Fisher's test p-values are shown in the table. Only p-values <0.05 are shown and only for overrepresentation of predictions rather than under-representation.

Outlier analysis
Two measures were used to identify statistical outliers on arcsinh transformed data: