Robust, coherent, and synchronized circadian clock-controlled oscillations along Anabaena filaments

Circadian clocks display remarkable reliability despite significant stochasticity in biomolecular reactions. We study the dynamics of a circadian clock-controlled gene at the individual cell level in Anabaena sp. PCC 7120, a multicellular filamentous cyanobacterium. We found significant synchronization and spatial coherence along filaments, clock coupling due to cell-cell communication, and gating of the cell cycle. Furthermore, we observed low-amplitude circadian oscillatory transcription of kai genes encoding the post-transcriptional core oscillatory circuit and high-amplitude oscillations of rpaA coding for the master regulator transducing the core clock output. Transcriptional oscillations of rpaA suggest an additional level of regulation. A stochastic one-dimensional toy model of coupled clock cores and their phosphorylation states shows that demographic noise can seed stochastic oscillations outside the region where deterministic limit cycles with circadian periods occur. The model reproduces the observed spatio-temporal coherence along filaments and provides a robust description of coupled circadian clocks in a multicellular organism.


INTRODUCTION
Endogenous circadian clocks allow the alignment of cellular physiology with diurnal light/darkness cycles on Earth, endowing organisms, from unicellular cyanobacteria to multicellular plants and mammals, with a selective fitness advantage (Cohen and Golden, 2015).Significant progress has been achieved in understanding circadian clock architectures and function in cyanobacteria, which are arguably, the simplest organisms exhibiting self-sustained oscillations (Hasty et al., 2010;Rust et al., 2007;Lambert and Rust, 2016;Dong et al., 2010;Teng et al., 2013;Gan and O'Shea, 2017).The molecular mechanisms behind autonomous circadian clocks have been elucidated primarily in the unicellular Synechococcus elongatus.These investigations have shown that the core of the circadian clock consists of three proteins, KaiA, KaiB and KaiC (Ishiura et al., 1998), whose oscillating behavior can be reconstituted in vitro (Nakajima et al., 2005).The basic mechanism behind the clock is based on the stimulation of KaiC autophosphorylation by the binding of KaiA, and the autodephosphorylation that ensues when KaiB binds to KaiC, blocking KaiA's stimulatory function.A salient feature of the circadian clock in Synechococcus is the high temporal precision it can exhibit, despite the fact that biochemical reactions in a cell are stochastic events and that clock components may be subject to variations in molecular copy numbers between cells, variations known as demographic noise (Tsimring, 2014).Many studies have addressed the robustness of circadian rhythms to demographic noise in Kai proteins (Mihalcescu et al., 2004;Chabot et al., 2007;Teng et al., 2013;Pittayakanchit et al., 2018;Chew et al., 2018), but copy number variations in KaiC phosphoforms, which impact directly on clock function, has not been previously considered.
The multicellular character of higher organisms and of some cyanobacterial species (Herrero et al., 2016), naturally prompts the question of how do ensembles of noisy circadian clocks perform in a multicellular organismal setting.Theoretical considerations have motivated the notion that reliable collective oscillations may result from the coupling of "sloppy", noisy clocks (Enright, 1980).It has been suggested that coupling of circadian clocks in unicellular organisms by quorum sensing interactions may result in emergent synchronization (Garcia-Ojalvo et al., 2004), and experimental evidence in support of this notion has been reported in a synthetic system (Danino et al., 2010).
Anabaena sp.strain PCC 7120 (henceforth Anabaena) is a multicellular cyanobacterium in which cells are arranged in a one-dimensional configuration, with local, nearest-neighbor cell-cell coupling through septal junctions (Herrero et al., 2016).Evidence of coupling of metabolic pathways along a filament due to cellcell communication has been reported (Mullineaux et al., 2008).In contrast to Synechococcus, information about the mechanism of the circadian clock of Anabaena is scant.Sequence BLAST analysis has shown that Anabaena contains homologs of the kaiA, kaiB, and kaiC genes of Synechococcus, and structural studies indicate that the interactions between the respective proteins are similar (Garces et al., 2004).Furthermore, Anabaena also contains factors coupling the kai post-transcriptional oscillator to input signals, and to output factors such as RpaA, CikA and SasA that couple the clock to the genes it regulates (Schmelling et al., 2017).The roles of these genes in Anabaena remain to be elucidated.A first characterization of the circadian rhythms in bulk cultures of Anabaena has shown that the clock is autonomous, running freely under constant light conditions following exposure to two 12-hour light-dark cycles, similarly to Synechococcus (Kushige et al., 2013).However, this study also showed that, in contrast to Synechococcus, none of the kai genes is expressed with a large amplitude.Nonetheless, about eighty kai-controlled genes that exhibit high-amplitude circadian oscillations were identified using DNA microarray analysis, a behavior that was abolished in a kaiABC deletion mutant.
Here we present an experimental and theoretical study of circadian clocks in multicellular Anabaena.Its one-dimensional character allowed us to follow clocks in each and every cell along a filament, and shed light on the interplay between demographic noise and cell-cell communication in setting synchrony and spatial coherence along filaments.In our experiments we followed in vivo the output of clocks in individual cells by monitoring the expression from the promoter of pecB, a clock-controlled gene of high amplitude oscillations (Kushige et al., 2013).This gene is part of the pecBACEF operon, and codes for the beta subunits of phycoerythrocyanin, a structural component of the phycobilisome rod that plays a major role in light harvesting for photosynthesis (Swanson et al., 1992).On the theoretical side, we first incorporated the effects of demographic noise into a three-component model of a single clock describing the phosphorylation states of KaiC, as in Synechococcus.Next, we extended this single-cell stochastic model to describe a one-dimensional array of coupled clocks, as in multicellular Anabaena, allowing us also to analyze the spatio-temporal coherence properties of noisy oscillations in Anabaena filaments.

Circadian clocks of individual cells in growing filaments are highly synchronized
We followed circadian oscillations in Anabaena from a chromosomal gfp fusion to the N-terminus of the clock-controlled protein, PecB, here denoted as P pecB−gf p at the level of individual cells along filaments.
Prior and during the experiments, filaments were grown under constant light conditions.Typical images of wild-type filaments at different time points are shown in Fig. 1A (see also Figure 1 Video 1).One salient feature in the images is significant synchrony, i.e. cellular oscillations progressed in individual cells along filaments together and with a similar period.The images in Fig. 1A correspond to successive maxima and minima of the mean cell fluorescence intensity, µ, as a function of time, which we plot in Fig. 2A.Similar experiments with a strain in which the kaiABC genes were deleted (∆kaiABC) resulted in a low-level, nonoscillatory signal (Fig. 2A), confirming the regulation of pecB expression by the circadian clock genes.As a control, we monitored expression from the promoter of hetR, a gene coding for the master regulator of development in Anabaena.The expression from the hetR promoter P hetR−gf p did not exhibit an oscillatory behavior (Fig. 2A), a result consistent with previous microarray experiments (Kushige et al., 2013).This does not preclude a possible interaction between the circadian clock and differentiation.On the other hand, the autofluorescence from photosynthetic pigments did not display oscillatory behavior (see Fig. 1B, and Fig. 2A).To characterize quantitatively the degree of synchronization between cellular clocks along a filament, we used the synchronization index R (Garcia-Ojalvo et al., 2004) (Materials and Methods), which can be readily calculated from the measured fluorescence intensity in each cell, and which varies between 0 (no synchronization) and 1 (full synchronization).To compute R, several cells along a filament were selected, and their fluorescence intensity was followed over a full period of oscillation.Contiguous cells, which include sister cells from the same mother are highly synchronized, as shown by the large value of R (0.89 ± 0.04) (Table 1).The value of R for cells initially separated by intervals of ten other cells -chosen in an attempt to avoid initial correlations between their clocks-was significantly indistinguishable from that computed for contiguous cells, underscoring the large degree of synchronization of clocks along a filament.
To test the extent to which different filaments are synchronized under the same conditions, we measured the average fluorescence intensity per cell in different filaments as a function of time (Fig. 2B).The expression from P pecB−gf p in different filaments oscillated nearly in phase mainly during the first oscillations.To evaluate quantitatively the degree of phase synchronization between filaments, we calculated R by choosing one cell per filament in a number of filaments measured simultaneously (Materials and Methods).We obtained R = 0.75 ± 0.04, a value that is significantly smaller than that obtained for cells within the same filament (Fig. 2C and Table 1).We surmise that initial synchronization may be due to phase resetting following the change in conditions e.g.illumination, upon transfer of cells from bulk culture to the microscope for realtime measurements.This change also could account for the decay in fluorescence intensity observed during the first circa 24 hours of our experiments (Fig. 2A and B).Furthermore, this decay is P pecB−gf p -specific, but clock-independent, as it was also observed in filaments of the ∆kaiABC background (Fig. 2A).

Circadian clocks along filaments are coupled by cell-cell communication
Another salient feature in the snaphots of Fig. 1A is the high spatial coherence of the expression from P pecB−gf p along filaments, i.e. all cells have nearly the same phase along their circadian cycle.To quantify the extent to which clocks are actually correlated, we calculated the spatial autocorrelation function of fluorescence intensity g as a function of distance along a filament (see Methods).We found that g decays to zero for separations of five cells or more (Fig. 2D).
To evaluate the contribution of phase inheritance following cell division to the observed autocorrelation, a simulated filament was generated from each measured filament by dividing each cell into two for three generations, partitioning the fluorescence of a mother cell binomially between the two daughters, and then multiplying the results by two, in order to conserve the average fluorescence per cell.The autocorrelation functions of these simulated filaments were then computed and averaged.The resulting mean autocorrelation (Fig. 2D, magenta line) decreased to zero already at the second neighbor, suggesting that another factor, e.g.cell-cell communication, contributes significantly to the coupling of fluctuations of pecB expression along a filament.To support this notion further, we calculated the spatial autocorrelation function of P pecB−gf p expression in filaments of a strain ∆sepJ/∆fraCD in which genes coding for three septal proteins SepJ, FraC and FraD involved in cell-cell communication were deleted (Nürnberg et al., 2015).Circadian oscillations in expression from P pecB−gf p in filaments of this strain were observed (Figure 2 -figure supplement 1A).
The resulting spatial autocorrelation function decreased over significantly shorter lengthscales (about 2 cells) than the wild type (Fig. 2D).Therefore, we calculated the synchronization index between contiguous cells in this genetic background and found that R was significantly smaller than that measured for contiguous cells within filaments of the wild-type background, but comparable to that obtained for cells in different filaments (Table 1).Consistently with this smaller value of R, traces of individual cells and their respective lineages were considerably more noisy (see Figure 2 -figure supplement 1B) than lineages in wild-type filaments (Fig. 2C).These findings indicate that the correlation in P pecB−gf p expression is due primarily to significant coupling between the clocks in neighboring cells, due to cell-cell communication.Of note, the value of R calculated for a ∆kaiABC background, in which P pecB−gf p expression is clock-independent, was similar to that for ∆sepJ/∆fraCD, in which cell-cell communication is impaired (Table 1).

Cell-cell variability oscillates out of phase with the circadian rhythm
Variations in gene expression between cells along a filament may limit both synchrony and spatial coherence.
These variations are evident in Fig. 1A even though their amplitude is small relative to the clock-modulated activity of P pecB−gf p .To quantify these variations, we calculated the square of the coefficient of variation CV = σ/µ, were σ denotes the standard deviation of the fluorescence intensity of contiguous cells along a filament.A plot of CV 2 as a function of time showed that the cell-cell variability itself displays oscillatory behavior, attaining maxima approximately in the middle of periods during which the mean fluorescence intensity increases (Fig. 2A).

Coupling between cell-cycle and clocks
In a cellular setting, circadian and cell cycle oscillations take place concurrently.In a variety of organisms, from prokaryotes to mammals, the circadian clock has been observed to gate cell division, enabling cell division to take place at some phases of the circadian cycle but inhibiting at others (Mori, 2009;Yang et al., 2010).The synchronization index R for strains with the indicated genotypes (Materials and Methods) was measured from the fluorescence intensities of P pecB−gf p expression in the same cells followed over a full circadian period in a filament, either in clusters of contiguous cells (contiguous), or for cells separated by intervals of ten cells (separate).To measure synchronization between filaments, R was computed from about ten cells, each belonging to a different filament.For each genetic background, the mean and standard error of the mean (SEM) of R was determined from a number of independent repeats (3-6), carried out in n different experimental runs.Significance (P-value) in interstrain comparisons was established by the Mann-Whitney U-test, and an asterix (*) represents rejection of the null hypothesis that samples come from distributions with equal medians.the time at which cell division takes place along the circadian cycle in individual cell traces, under conditions of constant illumination.The fluorescence intensity traces of two contiguous ancestor cells bearing the P pecB−gf p fusion and their respective lineage are shown in Fig. 2C.In Fig. 2E we show a histogram of the timing of cell division events as function of the phase at which they occur along the circadian clock, obtained from traces similar to those in Fig. 2C.Far from being equiprobable along the circadian cycle, cell division events showed a marked tendency to occur near the beginning or the end of a circadian cycle (minima in µ) as for Synechococcus (Yang et al., 2010).Fig. 2C also shows that the clock phase was faithfully inherited by any two daughters following cell division.
Transcriptional oscillations of kai genes and the master transducer regulator rpaA Previous Northern blot measurements and microarray analysis of kai genes showed no reliable, high-amplitude oscillatory expression of any of the kai genes in Anabaena.To study with higher sensitivity the expression of kai genes, we carried out RT-qPCR measurements of wild-type and ∆kaiABC strains in bulk cultures.Since the value of the R index was smaller between filaments than within (Table 1), the experiments were carried out under constant light conditions following two 12 h/12 h light-dark cycles, to enhance synchronization (see Materials and Methods).The relative expression of the three kai genes indeed showed noisy, low amplitude temporal modulations (Fig. 3A).The oscillations were largely in phase, and transcription occured mainly during subjective day.To assess quantitatively the extent of coordination, we calculated the synchronization index R for the different pairs and obtained R A,B =0.85 ± 0.05, R A,C =0.79 ± 0.09 and R B,C = 0.87 ± 0.06.Error bars were obtained by bootstrap methods (Methods).Since the differences between these values are not significant, we conclude that transcription of the three genes is coordinated.In fact, the value of R evaluated for the three genes was 0.78 ± 0.08.
To expose an underlying oscillatory behavior in the kai genes data and support the notion that the oscillations are circadian, we applied persistent homology methods (Pereira and de Mello, 2015;Otter et al., 2017) to two-dimensional phase portraits of their respective time series (see Appendix 1 -supplemental meth-ods).The delay τ for each phase portrait corresponds to the first minimum of the auto-mutual information of the time series (Fraser and Swinney, 1986), and for a periodic, nearly sinusoidal signal, corresponds to a quarter of the signals period (Kennedy et al., 2018).We obtained τ = 7.1 ± 1.2 hr, τ = 6.7 ± 1.2 hr and τ = 7.3 ± 1.1 (n ≥ 3) hr for kaiA, kaiB and kaiC respectively, all consistent with a circadian period of oscillation (Appendix 1).A Vietoris-Rips filtration of the clouds of points in the phase portraits (Appendix 1), indeed showed evidence for a persistent cycle in the transcription of each of the kai genes (Figure 3 figure supplement 1).
To shed light on how the state of the clock is relayed to the genes it controls, e.g.pecB, we measured the relative expression of rpaA in wild-type and ∆kaiABC strains by RT-qPCR.RpaA has been shown to transduce the phosphorylation state of KaiC to clock-controlled genes in other cyanobacteria (Markson et al., 2013;Iijima et al., 2015), and is highly conserved (Schmelling et al., 2017).The role RpaA plays in the circadian oscillations of Anabaena is unknown.We found that the relative expression of rpaA displays high amplitude oscillatory behavior (Fig. 3B).Furthermore, pecB displays oscillatory behavior, with similar amplitude and phase.The oscillatory behavior of both genes was abolished in the ∆kaiABC mutant (Fig. 3B     and C).Note that the expression of both rpaA and pecB peaks during subjective night.

Incorporating demographic noise into a model of a single clock
In order to develop a model of an array of coupled noisy clocks as in Anabaena, we first characterized the spectral properties of uncoupled clocks in individual cells of Synechococcus.We carried out experiments following the expression of YFP from the kaiBC promoter in single cells of Synechococcus growing within patterned gels (Fig. 4A).Circadian oscillations in the lineages of two sister cells are shown in Fig. 4B (see also Fig. 5B for the associated power spectrum).Expression from the kaiBC promoter exhibited circadian oscillations with a period of about 25 hrs, similar to those observed previously (Teng et al., 2013).Our next goal was to generalize a deterministic model for individual clocks in Synechococcus (Rust et al., 2007), to include the effects of demographic stochasticity.We followed the interaction network depicted in Fig. 5A (adapted from (Rust et al., 2007)).Our mathematical model takes into account the dynamics of three forms of KaiC, namely the single-phosphorylated forms S-KaiC (phosphorylated at Serine 431), T-KaiC (phosphorylated at Threonine 432), and the double-phosphorylated form D-KaiC, while the unphosphorylated U-KaiC can be deduced from the conservation of the total number of molecules of KaiC.Crucial for the appearance of oscillations is the negative feedback mediated by S-KaiC, through inactivation of KaiA via KaiB function.
Furthermore, the condition on the active KaiA monomers (see Eq. ( S4) in (Rust et al., 2007)) is modeled here by a continuous function (Appendix 1 -Eq.( 6)) that makes analytical progress possible.The parameter γ in our non-linear phosphorylation and dephosphorylation rates k XY for X, Y = {U, T, D, S} sets the steepness of the inverted sigmoidal dependence of the rates on KaiA (Figure 5 -figure supplement 1).
The dynamics of the stochastic model is described by a master equation accounting for discrete variations in molecular copy numbers of phosphoforms, instead of deterministic, ordinary differential equations (Rust et al., 2007;Brettschneider et al., 2010).We then use the van Kampen system-size expansion to carry out a linear noise approximation that yields, to leading order, an extended set of ordinary differential equations for ).The analysis of these equations allows us to determine the region in parameter space within which the model exhibits sustained deterministic oscillations, Fig. 5C.The parameters of the model have been set to the values determined in vitro in (Rust et al., 2007), and reported in Appendix 1 -table 1, except for γ and [KaiA] which were allowed to change freely.Note that deterministic oscillations with a circadian period were limited only to a small strip near the stability boundary.At the next-to-leading order, the expansion allows to evaluate the effects of demographic noise and to calculate the power spectrum of fluctuations for each species' abundance, due to finite size effects (see Appendix 1).A fit of the theoretical power spectrum to the experimentally measured one provides an adequate interpolation, upon adjusting the two fitting parameters, γ and [KaiA], see Fig. 5B.
In drawing the comparison between theory and experiments we assumed that the fluorescence intensity is an immediate proxy of the phosphoform T-KaiC (Teng et al., 2013).The fitted parameters position the systems outside the region of deterministic oscillations (circle in Fig. 5C), suggesting that circadian rhythms can be a manifestation of noise-driven oscillations (Fig. 5D).Remarkably, the fitted value for [KaiA] = 1.308µM matches the concentration reported previously ([KaiA] = 1.3µM , see (Rust et al., 2007)) .

Theoretical model of arrays of coupled noisy clocks
Next, we generalized the single-clock model above to an array of coupled circadian clocks as in Anabaena, which is endowed with cell-cell communication via septal proteins (Herrero et al., 2016) (Fig. 6A).We postulate that the intercellular transfer of factors such as sugars (Mullineaux et al., 2008;Nürnberg et al., 2015) may affect the behavior of neighboring clocks, yielding an effective long-ranged coupling between clock inputs across the filament.The interaction is here modeled by an exponential kernel which extends over a few cells.We assumed that the core clock mechanisms of Synechococcus and Anabaena are similar (Kushige et al., 2013;Schmelling et al., 2017;Garces et al., 2004), and followed the interaction network depicted in Fig. 5A.We further assumed that the rates of interconversion between KaiC phosphoforms for Anabaena are similar to those measured previously (Rust et al., 2007), in a reconstituted in vitro system consisting of KaiA, KaiB and KaiC of Synechococcus (Appendix 1 -table 2).This is supported by the fact that individual cells in both organisms exhibit oscillations with circadian periods (Figs.2C and 4B) that in the case of Synechococcus, are rather insensitive to changes in KaiC concentrations (Chew et al., 2018), constraining the values of these rates.Furthermore, KaiA of Anabaena, which is about two thirds shorter than KaiA of Synechococcus, is able to similarly dimerize and enhance the phosphorylation of KaiC of other cyanobacteria in vitro, as well as elicit oscillations when its gene is transferred to Synechococcus cells, despite the evolutionary divergence between both organisms (Uzumaki et al., 2004).Individual cells can therefore display self-sustained oscillations only if the parameters γ and [KaiA] take values inside the colored region of Fig. 5C.These oscillations are characterized by a circadian period only within a narrow strip near the stability boundary.The fluorescence intensity is assumed to reflect the output of the clock, with a phase difference as shown previously (Kushige et al., 2013).We hence calculated the power spectrum of the fluorescence signal on every cell along the filament and averaged together the results.The experimentally computed power spectrum shows a clear peak, which is nicely interpolated by the theoretically-predicted curve (Fig. 6C) upon adjusting the two fitting parameters, γ and [KaiA].Again, the fitted values position the system (diamond in Fig. 5C) outside the region of deterministic order, suggesting that demographic noise could trigger the observed oscillations.To test the robustness of the fit we have implemented a bootstrap procedure (see Materials and Methods) by perturbing each kinetic parameter with respect to the values reported in (Rust et al., 2007).The imposed perturbation was about 10%, for each individual kinetic parameter.The obtained best fit estimates for [KaiA] and γ were found to display a degree of variability of approximately 20%, with reference to averaged values reported in Fig. 5C.The bifurcation line which sets the separation between the deterministic limit cycle and the stationary stable fixed point became modulated depending on the set of assigned kinetic constants.Each pair of fitted γ and [KaiA] falls by definition in the domain where the deterministic oscillations are impeded and the stochastic finite size corrections drive the emergence of the resonant quasi cycles.
In Anabaena, the oscillations displayed by different cells along a filament are synchronized, most likely by cell-cell communication.To support further the notion of clock coupling, we studied the spatial coherence of noisy oscillations along an Anabaena filament.We assume that P pecB−gf p expression reports faithfully the output of the clock, albeit with a delay due to the transduction of KaiC's phosphorylation state up to activation of the promoter, probably by phosphorylated RpaA.To study spatial coherence, we calculated the complex coherence function (Marple, 1987), whose magnitude measures the degree of correlation between cells within a filament at different distances.As shown in Fig. 6D, the coherence at the characteristic frequency of .Lines between symbols are a guide to the eye.The fit was carried out by adjusting two parameters, the strength of the imposed spatial coupling and the characteristic scale of the exponential kernel, see Appendix 1. Remarkably, the range of the interaction as obtained from the fit is compatible with that estimated from the spatial autocorrelation depicted in Fig. 2D.The inter-cell coupling was obtained from the fit in panel C.
division alone cannot account for all the spatial variation of expression along filaments (Fig. 2D).In Fig. 6C quasi-cycles recorded on different cells of the stochastic Anabaena model are depicted showing a remarkable degree of synchronization.In Appendix 1 the analysis is extended so as to account for both a constant and a decaying power-law kernel of interaction.

DISCUSSION
Anabaena, a model organism in which each cell has a well-defined number of neighbors with which it communicates, has enabled us to study one-dimensional arrays of circadian clocks in a multicellular organism in space and time, by interrogating each cell individually.Remarkably, our experiments using a fluorescence reporter fused to the N-terminus of a clock-controlled gene show that filaments display large-amplitude oscillations, consistently with previous studies carried out in bulk cultures (Kushige et al., 2013).These oscillations, which run freely under constant light conditions and are therefore autonomous, are characterized by high spatial coherence and synchrony.In addition to inheritance following cell division as in unicellular Synechococcus (Amdaoud et al., 2007), the behaviors of both the spatial autocorrelation (Fig. 2D), and the complex coherence function (Fig. 6D) represent strong evidence that these two characteristics result from local coupling between clocks of neighboring cells due to cell-cell communication via septal junctions.
It is unlikely that this coupling results from the direct cell-to-cell transfer of KaiA, KaiB or KaiC clock components, as septal junctions in Anabaena are known to serve as conduits of only small molecules including metabolites, nutrients and small peptides (e.g.PatS-derived peptides involved in lateral inhibition during developmental pattern formation, mediated by SepJ (Corrales-Guerrero et al., 2015)).We postulate that metabolic factors such as sugars, which are transferred between cells via FraC and FraD proteins (Mullineaux et al., 2008;Nürnberg et al., 2015), may affect the behavior of neighboring clocks, yielding an effective long-ranged coupling between clock inputs (e.g.redox, glycogen abundance and ATP/ADP ratio) across the filament (Cohen and Golden, 2015;Golden, 2020).This is consistent with the shorter decay of the spatial autocorrelation function of P pecB−gf p fluorescence intensity, and with the reduced value of the synchronization index, measured in filaments of a strain in which the transport of sugars and/or other coupling factors is impaired (∆sepJ/∆fraCD).Of note, the value of the synchronization index in this strain is similar to that obtained for cells from different filaments.Thus, we have found that cells within a filament behave coordinately, whereas different filaments appear to be in different oscillatory phases.Nonetheless, coherent oscillations at the whole culture level were observed after light/darkness training of the cultures, showing that input signals can coordinate the circadian rhythm in cells from different filaments.
The high synchrony and spatial coherence observed in our experiments further suggest that circadian rhythms in Anabaena are not centrally coordinated as in higher multicellular organisms (Bell-Pedersen et al., 2005).In mammals for instance, clocks in peripheral tissues are centrally coordinated from a central pacemaker in the brain, the suprachiasmatic nucleus (Reppert and Weaver, 2002).In plants such as Arabidopsis, recent studies revealed the existence of waves of gene expression across the whole plant and the response of cells to positional information (Gould et al., 2018), consistent with weak coupling between cells and supporting a more decentralised model of clock coordination in plants (Endo, 2016).High synchrony and spatial coherence may be viewed as a necessary adaptation following the transition from a unicellular to a multicellular lifestyle (Masuda et al., 2017), while preserving key architectural features, gating of the cell cycle by the circadian clock and robust response to stresses, supporting the notion that a filament represents the organismic unit in Anabaena.
Our single cell measurements of gene expression from P pecB−gf p indicated that the cell-cell variability along filaments as measured by CV 2 is oscillating, achieving maxima approximately half-way during periods of increase of the fluorescence intensity.The phase difference between the oscillation in fluorescence intensity and its cell-cell variability is consistent with the phase of rpaA transcriptional oscillations observed in our RT-qPCR measurements.The oscillatory nature of CV 2 in P pecB−gf p expression may be due to the relay of the signal of the core clock by the transcriptional oscillations of a transcription factor, here RpaA (Heltberg et al., 2019).
Our experiments, carried out under constant light conditions i.e. non-varying external cues, provide clear evidence of gating of the cell cycle by the circadian clock, as for Synechococcus and many other organisms (Mori, 2009).In Synechococcus, cell doubling times vary considerably with light intensity (Teng et al., 2013), and the cell cycle has no effect on the circadian clock, irrespective of the cell cycle rate (Mori et al., 1996;Mori and Johnson, 2001).Furthermore, as for Synechococcus (Mihalcescu et al., 2004;Mori et al., 1996;Yang et al., 2010), clock phase is faithfully inherited by any two daughters following cell division as illustrated in Fig. 2C.In fact, the post-transcriptional design of cyanobacterial circadian circuits has been suggested to provide insulation from effects due to variable cell division rates (Paijmans et al., 2016).Cell doubling times in Anabaena can also vary significantly with light intensity (Zhao et al., 2007).The mechanism of how the circadian clock controls the timing of cell division has been studied in Synechococcus as well as in Synechocystis, and it has been suggested that phosphorylated RpaA regulates the bacterial tubulin-analog FtsZ, inhibiting the formation of the cytokinetic ring (Dong et al., 2010;Kizawa and Osanai, 2020).While the precise genetic and biochemical differences between the circadian clocks of Anabaena and Synechococcus remain to be elucidated, the presence of homologs of core clock components, output coupling factors such as RpaA (Schmelling et al., 2017) and cell division promoters such as FtsZ in Anabaena (Corrales-Guerrero et al., 2018), suggest that the mechanisms for gating in these two organisms may be similar.Note that ftsZ has an upstream putative RpaA binding site motif (see Figure 3 -figure supplement 2).Interestingly, nitrogen-fixing cells -heterocysts-that are formed under nitrogen deprivation in Anabaena, lose the ability to divide, and yet, the heterocyst-enriched fraction in bulk experiments has been shown to exhibit clear circadian oscillations (Kushige et al., 2013).This supports the notion that cell division does not affect the circadian clock in Anabaena as in Synechococcus.
The RT-qPCR results and their analysis, provide conclusive and quantitative evidence for circadian oscillations in the expression of kai genes in Anabaena.The calculation of the synchronization index R for different pairs of kai genes yields large and similar values, suggesting that the expression of the three genes is highly coordinated, consistent with their possible expression as an operon.Furthermore, a persistent homology analysis of the time series of the three genes yields clear evidence of oscillatory behavior, consistent with a circadian period.In agreement with microarray measurements (Kushige et al., 2013), the oscillations we detect are of small amplitude, similar to those of Synechocystis (Kucho et al., 2005;Beck et al., 2014), but unlike those of Synechococcus.Given that the oscillations in the expression of kai genes is small, but that of pecB and many other targets is large (Kushige et al., 2013), we reasoned that a possible way of transducing the core clock signal to clock-controlled genes may be furnished by the downstream transcription factor RpaA.The kai-dependent transcriptional oscillations of rpaA we observed in Anabaena suggest an additional level of regulation of clock-controlled genes, which differs from the post-translational mechanism observed in Synechococcus.In particular, the oscillations in kai genes is consistent with the existence of a transcription-translation feedback loop via RpaA (Markson et al., 2013).Indeed, a well defined candidate binding site motif of RpaA is located upstream of the kai genes (Figure 3 -figure supplement 2).Additional candidate binding site motifs of RpaA are found upstream of the rpaA locus itself as well as the pecBACEF operon.The large amplitude oscillatory behavior of rpaA transcription we observed remains to be elucidated.
In order to describe theoretically circadian oscillations in individual Anabaena filaments, we started by incorporating demographic noise into a deterministic model of a single clock that includes only the kai genes as a core, such as the one describing oscillations in Synechococcus, and then extended the model to one describing a one-dimensional array of coupled noisy clocks as in Anabaena.The stochastic models capture well the peaks of finite width in the power spectra of fluctuations observed in our experiments (Figs.5B and     6C).Importantly, the parameters we obtain from the fits to the experimental power spectra lie well outside the region in parameter space where deterministic oscillations occur, and in particular, the smaller region in response to combined nitrogen deprivation (Herrero et al., 2016;Di Patti et al., 2018).In this context, note that Anabaena filaments display circadian oscillations under nitrogen-deprived conditions, as shown previously (Kushige et al., 2013).Lastly, the model also reproduces the spatial coherence of oscillations in filaments, providing independent evidence of clock coupling through cell-cell communication.
In spite of the fact that Kai protein copy numbers run in the thousands in Synechococcus, phase fluctuations between cells are readily visible both in our experiments as well as in previous ones (Chew et al., 2018;Chabot et al., 2007).Stochastic modeling indicates that these numbers may be needed to compensate for noise amplification introduced by the post-translational feedback loop provided by KaiA, and reducing the numbers of Kai proteins leads to lower clock precision (Chew et al., 2018).We point out that demographic noise may be more significant than the copy numbers of Kai proteins may suggest: the functional form of KaiC is hexameric, which further partitions into its different phosphoforms (Chew et al., 2018).In addition, it has been shown that KaiB spends most of its time in an inactive tetrameric state, preventing its interaction with other Kai proteins, and that KaiA also oscillates between active and inactive states by binding to KaiB (Tseng et al., 2017).Together, these regulatory contributions might reduce active KaiC effective numbers, thereby increasing fluctuations.
While the copy numbers of Kai proteins in Anabaena have not been measured, it is likely that they are significantly smaller than in Synechococcus, given the low transcriptional levels observed in our RT-qPCR experiments, as well as in DNA microarrays (Kushige et al., 2013).The expression of kai genes in Synechocystis is weak (Kucho et al., 2005), and protein copy numbers are correspondingly small (Zavřel et al., 2019).Thus, demographic noise effects may be indeed at least as important in Anabaena as they are for Synechococcus.One may hypothesize that cell-cell communication and the resulting coupling of clocks, compensate for the smaller number of Kai proteins, setting the high synchrony and spatial coherence we observe in Anabaena.The picture that emerges may be more general and applicable to other multicellular organisms (Gould et al., 2018).
The transition from unicellular to multicellular organisms demanded coordination between physiological processes in different cells, in order to enhance organismal fitness and adaptation to stresses.This is true in particular of circadian clocks and the genes they regulate.By analyzing filaments at the individual cell level, we found concrete evidence supporting coordination mediated by cell-cell communication, allowing clocks in different cells to be coupled.Furthermore, the high synchrony and spatial coherence of cell clocks along filaments observed in the experiments suggests that cell-cell communication contributes significantly to these two properties.Our theoretical model of single core clocks and arrays thereof shows that far from being detrimental, demographic noise may seed oscillations that can be synchronized by clock coupling.This provides a robust description of circadian oscillations in a multicellular organism such as Anabaena.

Strains
Strains, bearing a chromosomally encoded P pecB−gf p , were obtained by conjugation with the following backgrounds: wild type (WT) Anabaena sp.(also known as Nostoc sp.) strain PCC 7120; ∆sepJ/∆fraCD, with the strain CSVM141 in which sepJ, fraC and fraD were deleted (Nürnberg et al., 2015); ∆kaiABC in which the kaiABC genes were deleted (for details see Appendix 1 -Supplemental methods), as recipients.

Culture conditions
Strains and derived strains were grown photoautotrophically in BG11 medium containing NaNO 3 , supplemented with 20 mM HEPES (pH 7.5) with shaking at 180 rpm, at 30 • C, as described previously (Corrales-Guerrero et al., 2013;Di Patti et al., 2018).Growth took place under constant illumination (10 µmol m −2 s −1 ) of photons (spectrum centered at 450 nm) from a cool-white LED array.When required, streptomycin sulfate (Sm), and spectinomycin dihydrochloride pentahydrate (Sp) were added to the media at final concentrations of 2 µg/mL for liquid and 5 µg/mL for solid media (1% Difco agar).The densities of the cultures were adjusted so as to have a chlorophyll content of 2-4 µg/mL 24 h prior to the experiment, following published procedures (Di Patti et al., 2018).For time lapse measurements, filaments in cultures were harvested and concentrated 50 fold.Synechococcus cultures were grown as above, and when required, 7.5 µg chloramphenicol (Cm) was added.

Samples for time-lapse microscopy
Strains were grown as described previously (Di Patti et al., 2018).When required, antibiotics, streptomycin sulfate (Sm) and spectinomycin dihydrochloride pentahydrate (Sp), were added to the media, at final concentrations of 2 µg/mL for liquid and 5 µg/mL for solid media.The densities of the cultures, grown under an external LED array (15 µmol m −2 s −1 ) for about five days, were adjusted so as to have a chlorophyll a content of 2-4 µg/mL, 24 h prior to the experiment following published procedures (Di Patti et al., 2018).
For time-lapse, single-cell measurements of Anabaena, 5 µL of culture concentrated 100-fold were pipetted onto an agarose low-melting gel pad (1.5%) in BG11 medium containing NaNO 3 and 10 mM NaHCO 3 , which was placed on a microscope slide.The pad with the cells was then covered with a #0 mm coverslip and then placed on the microscope at 30 • C. The cells grew under light from both an external LED array (15 µmol m −2 s −1 ) and tungsten halogen light (10 µmol m −2 s −1 ), 3000K colour).Under these illumination conditions, the doubling time of cells is similar to that in bulk cultures (Di Patti et al., 2018).The change in illumination conditions when transferring cells from bulk cultures to the microscope results in high synchronization within filaments.Images of about ten different fields of view were taken every 30 min on a Nikon Eclipse Ti-E microscope controlled by the NIS-Elements software using a 60 N.A 1.40 oil immersion phase contrast objective lens (Nikon plan-apochromat 60 1.40) and an Andor iXon X3 EMCCD camera.Focus was maintained throughout the experiment using an Perfect Focus System (Nikon).All the filters used are from Chroma.The filters used were ET480/40X for excitation, T510 as dichroic mirror, ET535/50M for emission (GFP set), ET500/20x for excitation, T515lp as dichroic mirror, and ET535/30m for emission (EYFP set), and ET430/24x for excitation, 505dcxt as dichroic mirror, and HQ600lp for emission (chlorophyll set).
Samples were excited with a pE-2 fluorescence LED illumination system (CoolLED).
For measurements of Synechococcus, the cultures were grown as above to a OD 750 of 0.3-0.4.Samples were then entrained by two light-dark cycles (12 h-12 h) before measurements commenced.The cultures were then diluted to a OD 750 of 0.1 using BG11 medium, and 2 µL of the culture was pipetted onto a glassbottom culture dish.A patterned pad prepared as above but solidified on an optical grating (Hadizadeh Y. et al., 2012) was placed atop the cell suspension.Then, 1.5% agarose melted in BG11, cooled to 37 • C was poured on top of the well to cover the pad.After solidifying this last layer, 5 ml of BG11 was added to the culture dishes to maintain the moisture level of the agarose pad during the experiment.This device was placed in the microscope and time-lapse measurements were carried out as for Anabaena.

RT-qPCR measurements
For real time quantitative polymerase chain reaction (RT-qPCR) measurements, strains were grown under the same conditions as described above.When the cultures reached the beginning of their exponential phase, about 1.8 µg/mL of chlorophyll a content, they were entrained by two light-dark cycles (12 h-12 h) and then grown under constant light.Then, measurements were started 24 h after, and total RNA was extracted from the PCC 7120 and ∆kaiABC cultures in two biological replicates, every 4 h.For each sample, 20 mL of cells were collected by filtration and washed with buffer TE50 (50 mM Tris-HCl at pH 8.0, 100 mM EDTA).
Cells were resuspended in 2 mL TE50 buffer, and then they were centrifuged at 11,500 g for 2 minutes at 4 • C. Supernatant was then removed and cell pellets were flash-frozen in nitrogen liquid for storage at 80 • C.
Total RNA was isolated by using hot phenol as described (Mohamed and Jansson, 1989) with modifications.
All samples were treated with DNase I, and RNA quantity and purity were assessed with a Nanodrop One spectrophotometer as well as by agarose gel electrophoresis.
For cDNA synthesis, 600 ng of total RNA of each sample were used for reverse transcription with the QuantiTect Reverse Transcription kit (Qiagen).Then PCR amplification of 17.4 ng of each cDNA was carried out using Fast SYBR Green Master Mix (Applied Biosystems) that includes an internal reference based on the ROX TM dye, and specific primers (Appendix -table 2).The amplification protocol was one cycle at 95 • C for 10 min, 40 cycles of 95 • C for 15 sec, and 60 • C for 60 sec.RT-qPCR was carried out using an Applied Biosystems StepOnePlus instrument equipped with the StepOne Software v2.3.After the amplification was completed, a melting point calculation protocol was carried out in order to check that only the correct product was amplified in each reaction.Reactions were run in triplicate in three to four independent experiments from two biological replicates.The transcript levels of kaiA (alr2884), kaiB (alr2885), kaiC (alr2886), pecB (alr0523) and rpaA (all0129) were normalized to the transcript levels of the housekeeping genes rnpB and all5167 (ispD) to obtain the ∆∆ct value.Relative gene expression was calculated as 2 −∆∆ct (Pfaffl, 2001).

Image segmentation
All image processing and data analysis was carried out using Matlab (MathWorks).Filament and individual cell recognition was performed on phase contrast images using an algorithm developed in our laboratory.The program's segmentation was checked in all experiments and corrected manually for errors in recognition.The total fluorescence from GFP (for Anabaena) and chlorophyll (autofluorescence) channels of each cell, as well as the cell area, were obtained as output for further statistical analysis.

Analysis of synchronization and spatial correlations along filaments
Synchronization was measured by the order parameter proposed by Garcia-Ojalvo et al (Garcia-Ojalvo et al., 2004): where ⟨•⟩ denotes a time average, .indicates an average over all cells, and µ denotes the average of the fluorescence intensity of each cell f i .Hence R is defined as the ratio of the standard deviation of µ(t) to the standard deviation of f i , averaged over all cells.For the measurement of synchronization within the same filament, groups of 8-11 cells were chosen, whether separated or contiguous (sharing a common ancestor as determined from a lineage analysis).For evaluation of inter-filament synchronization, one cell per filament was chosen randomly in different fields of view.R was then calculated and this procedure was repeated for different choices of cells, at least three times for each experiment.All the evaluations of R were carried out over a full period of oscillation, in either one of the first two oscillations, except for the ∆kaiABC background, for which R was calculated for an interval of 24 hours, during which other strains display the first full oscillation.The final result comprises the mean of at least three independent repeats, in at least two independent experiments.Errors in the quoted values of R therefore represent standard errors (SEM).
Statistical analyses were performed in Matlab using Mann-Whitney ′ s U-test.
The spatial autocorrelation function of fluorescence intensities along filaments were calculated using the autocorr Matlab command with 30 lags, from at least 25 filaments of 50 cells each, and at least two independent experiments.For WT and ∆sepJ/∆fraCD backgrounds, autocorrelation functions of individual filaments were calculated at maxima and minima of circadian oscillations and the results were averaged.
For filaments of the ∆kaiABC background, autocorrelation functions were calculated at about 30 and 50 hours, corresponding to the first minimum and maximum of oscillations observed in filaments of wild-type background.The spatial autocorrelation function was calculated from the formula: where where f i denotes the fluorescence intensity in cell i in a filament of N cells.

Robustness of Anabaena dataset fit to the Synechococcus parameters
Each kinetic parameter (k) was randomly selected from a Gaussian distribution centered at the nominal value ( k) estimated by (Rust et al., 2007) using the normrnd Matlab function.The standard deviation is assigned as σ = k/10.For every complete set of (randomly selected) kinetic constants, we proceeded with the fit of γ and [KaiA] to interpolate the experimental power spectrum.Each pair of fitted values was stored to eventually compute averaged estimates, together with the error, as quantified by the associated standard deviation.By averaging over 200 independent realizations of the implemented procedure yields γ = 7.2 ± 1.4 and [KaiA] 3.4 µM Appendix 1 -table 1: Parameters used in the simulations.All the values are from the paper (Rust et al., 2007).
where the superscript α indicates that the species belongs to the cell α.The interconversion rates between phosphoforms are given by where [KaiA] (resp. [KaiC]) is the concentration of KaiA (resp.[KaiC]).The values of the rates are set as those of (Rust et al., 2007) and are listed in Appendix 1 -table 1.The variable n α S stands for the total number of S molecules in cell α at time t and, in the same way, we will use n α X to denote the total number of X molecules for X = T, D, U .The total number of KaiC phosphoforms remains constant in each cell of the filament, and we set this value equal to N .Note that KaiB does not appear explicitly in the equations as it is assumed that it is in excess at all times and it is further assumed that interacts with S-KaiC instantaneously (Rust et al., 2007).The typical sigmoidal shape of the nonlinear function ( 6) is plotted in Figure 5 -figure supplement 1 (blue dashed line), where we also plotted the corresponding function used in the paper by Rust et al. (Rust et al., 2007) (red line).
To model cell-cell communication we postulate a constant kernel that extends over few cells, as depicted in Appendix 1 -figure 1.This hypothesis will be relaxed later (Accounting for a smooth kernel of interactions).More specifically, and despite the fact that there is no experimental evidence for intercellular transport of KaiC, we assume that coupling between clocks is metabolic and we describe this effect through longrange interactions between molecules at different nodes.The parameter k α , as shown in Appendix 1 -figure 1, counts the number of neighbors of any given cell.The chemical equations corresponding to the envisaged coupling read where α and β label connected nodes.
Since the total number of molecules within each cell of the filament remains constant, the variable n α U can be neglected because it can be recovered from the conservation law For this reason, at each time t the state of the system is specified by the 3Ω-dimensional vector n(t) The microscopic dynamics described by rules ( 4) and ( 7) is Markovian, and thus the configuration of the system only depends on that at the previous time.We introduce the transition rates T(n ′ |n) from state n to a new state n ′ for all the chemical equations that constitute the stochastic model.To denote the new state n ′ we only explicitly write the components of the vector that change.Accordingly, the phosphorylation and dephosphorylation transition rates in each cell α read while the interconversion transitions are To account for cell-cell communication, we also add the transitions for long range interaction between two connected cells α and β: The dynamics is fully described in terms of the following master equation whose explicit version reads + T(n|n α S + 1)P (n α S + 1; t) − T(n α S + 1|n)P (n, t)+ The master equation can be written in a more compact form though the use of the step operators ϵ ±(α) i that act on a generic function g(n) as ϵ ∀i ∈ {T, D, S} and α = 1, . . ., Ω: The master equation cannot be solved analytically and therefore we apply the van Kampen system-size expansion (van Kampen, 2014) to analyze both the deterministic behavior and the role of demographic noise.For this purpose, we split each discrete variable n α X (t) into two contributions, the deterministic one, that is obtained taking the limit N → ∞ denoted by φ We introduce a new probability function Π(ξ, t) ≡ P n φ(t), ξ ; t where the ξ is a 3Ω-dimensional vector whose components are so that the left hand side of the master equation ( 15) can be recast in the form with τ = t/N .Operating an expansion of the step operators with respect to the small parameter 1/ √ N , performing the change of variable ( 16) into the right hand side of the master equation, and collecting terms proportional to 1/ √ N we get where and the matrix ∆ denotes the discrete Laplacian operator.The concentration φ U reads [KaiC]−φ T −φ D −φ S and [KaiC]φ X → φ X for the ease of notation.

Stability analysis
We set to zero the coupling parameter ν, which amounts to neglecting the spatial terms.In this way system (19) simplifies to We fix all the parameters to the values specified in Appendix 1 table -1, while [KaiA] and γ are used as control parameters and denote by φ * = (φ * T , φ * D , φ * S ) the homogeneous stationary solution of system (23).To characterize the nature of the fixed points we perform a linear stability analysis.For this purpose, for each positive φ * we evaluate the Jacobian matrix where and calculate its three eigenvalues.The nature of the equilibrium points and their stability is summarized in Appendix 1 figure 2 where we plot the equilibrium points to varying γ.Continuous lines denote stable equilibrium points, while dashed lines stand for unstable solutions.The diagram shows three different regions: in region I the system admits only one stable complex equilibrium point that becomes unstable as long as γ > γ 1 .The value γ 1 marks a supercritical Hopf bifurcation though which we enter in region II.Here we have a couple complex eigenvalues with positive real part and the trajectories approach a stable limit cycle.For γ > γ 2 the system undergoes a saddle-node bifurcation with the sudden creation of two additional fixed points, one stable and the other unstable (region III).
To better investigate the region of the parameters space where the deterministic system exhibits regular oscillations similar to those observed in experiments, we focused on an interval of γ including the two critical points γ 1 and γ 2 .In Appendix 1 -figure 3 we delimit with blue line the portion of the parameter plane (γ, [KaiA]) where the concentrations of three phosphoforms of KaiC exhibit regular oscillations.The period of oscillations corresponding to the smaller selected region delimited by the black rectangle is reported in Fig. 3B of the main text.From Fig. 3B of the main text, it appears that the region with period T ≃ 24 h is close to the separatrix between the stable fixed point dynamics and the limit cycle regions.A typical result of a numerical integration of system (23) for a set of parameters outside the limit cycle region, but close to the border, is shown in Appendix 1 -figure 4. Concentrations exhibit damped oscillations that eventually reach the stationary state.Performing stochastic simulations for the same parameters as those in Appendix 1 -figure 4 yields the results displayed in Appendix 1 -figure 5. Quasi-regular oscillations develop when the deterministic solutions attain the corresponding equilibrium value (McKane and Newman, 2005).To shed light onto this phenomenon is entirely devoted the forthcoming discussion.

Analysis of fluctuations
Considering the next-to-leading order of the van Kampen expansion and collecting together terms proportional to 1/N , we find the differential equation for the probability distribution Π(ξ, τ), that results in the following Fokker Planck equation where A is The two matrices M and B have dimension 3Ω × 3Ω, and they can be expressed as the sum of two distinct terms, one due to the spatial contribution that we label with the superscript (SP ), and the other one related to the interactions within each cell that we label as (N S).Moreover we assume that all the entries of the two matrices are evaluated at the equilibrium point φ * .Thus we have with Appendix1-figure 4: Deterministic simulation.Results of the numerical integration of system ( 23) for γ = 8 and [KaiA] = 1.2.The other parameters are assigned as specified in Appendix 1 -table (1).φ X for X = T, D, S stands here for the relative abundance of the phosphoforms in units of [KaiC].
where J is the Jacobian matrix relative to system (19) evaluated at the steady state.For B we have with

Spatial correlations
The Fokker-Planck equation ( 25) describes the fluctuations about the equilibrium point φ * and it is equivalent to the following coupled Langevin equations Appendix1-figure 5: Comparison between deterministic and stochastic simulations.The blue lines denote the result of the numerical integration of (23) while the noisy red lines represent the stochastic simulation of the system through the Gillespie algorithm (Gillespie, 1977).For all the plots for γ = 8 and [KaiA] = 1.2 while the other parameters are set as those in Appendix 1 -table 1. φ X for X = T, D, S stands here for the relative abundance of the phosphoforms in units of [KaiC].
where ξ α i is the component of vector ξ α defined in Eq. ( 17) while η α i is the ith component of the noise η α with the following properties Starting form this equation, it is possible to derive an analytic expression for the so called Power spectrum density matrix (PSDM), a diagnostic tool that allows to highlight the presence of correlations between species at different distances along the chain (Challenger and McKane, 2013;Zankoc et al., 2017Zankoc et al., , 2019)).The (PSDM) P ij,αβ (ω) is defined by where ξα i (ω) is the temporal Fourier transform of the ith component of ξ α (τ ).The degree of correlation between species i in cell α and species j in cell β can therefore be quantified through the Complex Coherence Function (CCF), a 3Ω × 3Ω matrix C whose entries are the ratio of the elements of PSDM to a normalization coefficient To obtain an analytic expression of P ij,αβ we solve the Langevin equation ( 31) performing a temporal Fourier Transform.In this way, for each component of the fluctuations vector on cell α we get with Φ(ω) = −iω − M, where denotes the 3Ω × 3Ω identity matrix.Substituting Eq. ( 36) into Eq.( 34) one eventually obtains where Φ † = (Φ * ) T is the complex conjugate transpose of Φ and use has been made of the following properties:

Accounting for a smooth kernel of interactions
In this section we will extend the analysis so as to account for a a generalized kernel of interaction, which decays with the inter-particle distance.More concretely, each node is coupled with all the other nodes along the chain, but the strength of the coupling is modulated as a monotonous decreasing function of the distance, namely f (|α − β|).Adjacent nodes are hence prone to mutual influences, as compared to nodes which are far apart from each other.The newly imposed modulation of reflects in the transition rates (12) as follows: where ν is a constant and f (|α − β|) is a decreasing function of the distance between nodes α and β.This results in a modification of the diffusive-like term in the master equation which can be reabsorbed in a new definition of the adjacency matrix W αβ as: and the corresponding Laplacian matrix ∆ ′ now becomes Accordingly, the two matrices ( 27) and ( 29) appearing in the Fokker-Planck equation ( 25) must be replaced by In the following we will focus on two different choices for the function f .Our first choice is to consider an exponentially decaying kernel: were σ > 0. The second option that we set to explore assumes a power law decay in the form: with ρ > 0. In the next section we will challenge different scenarios against experimental data.

Comparison with experiments
Setting ν = 0 in Eq. ( 37) and considering the diagonal terms α = β and i = j, returns N identical temporal power spectra of the fluctuations.The analytical prediction can be used to fit experimental data for both Anabaena (see Fig. 4C of the main text) and Synechococcus (see Fig. 3C of the main text).In both cases the values of γ and [KaiA] predicted by the fit lay outside the region where deterministic oscillations occur.Further, the CCF can be invoked to quantitatively explain the correlation of the fluorescence signal recorded in different cells of the Anabaena filaments.More specifically, we will adjust the theoretically predicted CCF to the analogous quantity estimated experimentally, by varying the coupling parameter ν and the connectivity k α when postulating a uniform range of interaction, or by varying ν and the characteristic parameter of the decay function f (|α − β|) a smoothly decaying kernel of interaction.
To calculate the experimental spatial correlation we proceeded as follows.We focused on a given portion made of Ω cells of the full Anabaena filament.We recall that this latter is growing in time.Hence, as time progresses, the cells within the examined portion of filament get replaced by their offspring as follows cellular duplication.However, for the analysis of spatial correlations, we are only interested in the relative positioning of the cells.In other words, we analyze the fluorescence content of each cell in relation to the position occupied within the filament irrespectively of the fact that the monitored population of Ω cells is constantly renewed in time by new offspring generation.By following this procedure we obtained Ω time series for the fluorescence on each spot along the chain.Then these temporal signals are transformed via a standard FFT in time.We further select the signal corresponding to ω * , the frequency at which the temporal power spectrum exhibits a peak and insert it in Eq. ( 34).
The results of the comparison between theory and experiments are plotted in Fig. 5 of the main text and in Appendix 1 -figure 6.The plot in the main text corresponds to the exponential long-range cell-cell communication (Eq.( 44)), where the best fit parameters are ν = 0.61 and σ = 0.35.Notice that 1/σ ≃ 3 set the characteristic scale of the spatial interaction (in units of cells).In Appendix 1 -figure 6 we compare the fit of experimental data (red circles) obtained by using a fixed interaction kernel with connectivity k α = 8 (black dashed line) and a decaying function (Eq.( 45)) (blue solid line).In all cases, the experimental points are nicely interpolated by the theoretical curves.Accounting for a smooth kernel of interaction yields more accurate fittings.
Appendix1-figure 6: The complex coherence function: comparison between theory and experiments.Complex coherence function measuring the correlation of 35-cell segments at the frequency of temporal oscillations.Red circles correspond to experimental data, the two lines represent the fit to the experimental data with the prediction of the stochastic model using the constant kernel (dashed black line) and the power law kernel (blue solid line).The fitting procedure gives ν = 0.49 for the fixed connectivity k α = 8, while it returns ν = 1.11 and ρ = 1.96 for the power law function.

Cloning and strain construction
For construction of a fusion of the promoter of pecB to the gfp-mut2 gene in Anabaena sp.(also known as Nostoc sp.) strain PCC 7120, a custom DNA sequence containing an EcorRI restriction site preceded by the pecB sequence of Anabaena was synthesized and inserted into the pUC57 plasmid (Bio basic) (Region alr0523).The pecB construct was amplified with the primer pair alr0523-EcoRI-Fw/alr0523-Rev.A pC-SAL33 plasmid, which bears the gfp-mut2 sequence was amplified with the primer pair 4G-GFP-Fw/alr0523-GFP-Rev. The PCR amplified sequence products gfp-mut2 and pecB were mixed and amplified with the primer pair (alr0523 PCR-Fw/GFP-PCR2-Rev).The merged PCR product and a pCSV3 plasmid were digested with EcoRI, and ligated to produce pSVRG4.The transferred sequence was confirmed by sequencing with the primers (alr0523-1-Fw, alr0523-1-Rev, 4G-gfp-Fw in plasmid PCSV3, and Rev plasmid PCSV3).The final construct bears a sequence of 1337 bp that comprises the pecB promoter and the 5' end of the pecB gene linked to the gfp-mut2 gene by a sequence encoding four glycines.The pSVRG4 plasmid was transferred to Anabaena by conjugation, which was performed as described (Elhai et al., 1997).Clones resistant to streptomycin (Sm) and spectinomycin dihydrochloride pentahydrate (Sp) were selected and their genomic structure was tested by PCR (alr0523(7120)-1/GFP-Rev). Periodically-performed PCR analysis indicated that the strains were stable, not showing wild-type chromosomes, in BG11 medium supplemented with antibiotics.The Anabaena strains used as recipients in the conjugations were: the wild-type PCC 7120, ∆sepJ/∆fraCD, in which sepJ, fraC and fraD were deleted (Nürnberg et al., 2015), and ∆kaiABC in which the kaiABC genes were deleted (see Appendix 1 -figure 7).Strain CSL64 bearing a chromosomally encoded P hetR−gf p transcriptional fusion in a wild-type background has been reported previously (Corrales-Guerrero et al., 2015).
The Synechococcus strain containing the gene encoding a YFP-SsrA reporter whose expression is driven by the kaiBC promoter at neutral site II was constructed by transforming Synechococcus sp.PCC 7942 with plasmid EB2316, which was a gift from Erin O'Shea (Addgene plasmid # 87753; http://n2t.net/addgene:87753;RRID: Addgene 87753).The transformation of Synechococcus was carried out as previously described (Golden and Sherman, 1984).The cells were incubated for 48 h at 30 • C under illumination on nitrocellulose filters, and transformants were selected on Cm-containing BG11 plates.Confirmation of genomic structure was carried out by PCR with the primers kaiC-4, kaiB-1, PkaiBC-1 and YFP-2.All oligo sequences used for construction are listed in Appendix 1 -table 2.  (Elhai et al., 1997), and sucrose (sacB)-based positive selection for double recombinants as described by (Cai and Wolk, 1990).Anabaena ∆kaiABC homozygous mutants containing the kai deletion with insertion of the C.K1 gene cassette in direct orientation were obtained and confirmed by PCR analysis.Kai-1 to Kai-4 are oligodeoxynucleotide primers.

Culture conditions
Strains and derived strains were grown photoautotrophically in BG11 medium containing NaNO 3 , supplemented with 20 mM HEPES (pH 7.5) with shaking at 180 rpm, at 30 • C, as described previously (Corrales-Guerrero et al., 2013;Di Patti et al., 2018).Growth took place under constant illumination (10 µmol m −2 s −1 ) of photons (spectrum centered at 450 nm) from a cool-white LED array.When required, streptomycin sulfate (Sm), and spectinomycin dihydrochloride pentahydrate (Sp) were added to the media at final concentrations of 2 µg/mL for liquid and 5 µg/mL for solid media (1% Difco agar).The densities of the cultures were adjusted so as to have a chlorophyll content of 2-4 µg/mL 24 h prior to the experiment, following published procedures (Di Patti et al., 2018).For time lapse measurements, filaments in cultures were harvested and concentrated 50 fold.Synechococcus cultures were grown as above, and when required, 7.5 µg

Figure 1 :
Figure 1: Circadian oscillation in Anabaena.(A) GFP fluorescence in a filament of an Anabaena strain bearing a P pecB−gf p promoter fusion, growing under nitrogen-replete conditions.The snapshots were chosen near maxima and minima of the circadian oscillations.(B) Autofluorescence as a function of time in Anabaena.Snapshots correspond to those in part A and time 0 corresponds to the time at which filaments were placed in a device for microscope observation (for details see Materials and Methods).For a time-lapse movie see Video 1 (taken over 6 days).

Figure 2 :
Figure 2: Characterization of a clock-controlled gene in Anabaena.(A) Average cell fluorescence intensity from P pecB−gf p in a filament as a function of time for a wild-type genetic background (full black circles) and for a ∆kaiABC background (empty black circles); Intensity of autofluorescence as a function of time (blue circles); Average cell fluorescence intensity from P hetR−gf p (cyan line); Temporal dependence of the cell-cell variability CV 2 in expression of P pecB−gf p (red line).Data were taken from at least 50 contiguous cells along a filament.(B) Average fluorescence intensity as a function of time of filaments in different fields of view from the same experimental run.Each trace was obtained from at least 50 contiguous cells along each filament.(C) Expression from a P pecB−gf p fusion in the lineages of two contiguous cells as a function of time.(D) Average spatial autocorrelation function of P pecB−gf p expression along filaments of wild-type (blue), ∆sepJ/∆fraCD (green) and ∆kaiABC (orange) genetic backgrounds.Error bars represent standard errors.Magenta line: contribution to the spatial autocorrelation function of fluctuations from the wild-type data set, induced by binomial partitioning of molecules between daughter cells, following each of three consecutive cell divisions.Prior to divisions, the cell order in each filament was reshuffled.(E) Histogram of the phase of cell-division events along the circadian cycles, with 0 and 2π denoting two consecutive minima, from two independent experiments.For additional data similar to panels A and C, corresponding to filaments of the sepJ/fraCD genetic background, see Figure 2 -figure supplement 1.

Figure 3 :
Figure 3: Transcriptional oscillations in the core clock genes, rpaA and pecB.(A) Relative expression of kaiA (green), kaiB (red) and kaiC (blue) as a function of time measured by RT-qPCR (Methods).A persistence homology analysis of these data is presented in Figure 3 -figure supplement 1. (B) and (C) Relative expression levels of rpaA and pecB respectively, in wild-type (full circles) and ∆kaiABC strains (empty circles).Curves have been normalized by their temporal mean.Error bars represent the standard error of the mean of three independent experiments (see Materials and Methods).Grey shades represent periods of subjective night.For additional information about regulatory sequences of the kaiABC, rpaA, pecB promoter regions and RpaA binding sites in Anabaena, see Figure 3c -figure supplement 2.

Figure 4 :
Figure 4: Circadian oscillations in Synechococcus.(A) Growth and lineage of a cell in patterned agarose, expressing YFP from the kaiBC promoter.The snapshots were chosen near maxima and minima of the circadian oscillations.(B) Fluorescence intensity of YFP of individual cells obtained from two independent cell lineages (red and blue).

Figure 5 :Figure 6 :
Figure 5: Stochastic model for circadian oscillations in S. elongatus.(A) Schematic representation of interconversion between KaiC phosphoforms modulated by the activity of KaiA in an individual clock.The different phosphophorm states of KaiC are denoted by: U (unphosphorylated, U-KaiC); T (phosphorylated at Threonine, T-KaiC); S (phosphorylated at Serine, S-KaiC) and D (phosphorylated at both sites).Arrows denote transitions between the different phosphophorms X, Y with the indicated rates k xy .KaiB mediates the inactivation of KaiA by S, as described by the continuous function f (see Figure 5 -figure supplement 1).(B) Average power spectrum of single cell fluorescence (red symbols), fit to the data with the prediction from the stochastic model (blue line).(C) γ-[KaiA] plane where deterministic limit cycle oscillations in individual clocks occur.The color corresponds to the period of oscillations (in hours).The boundary of the colored region corresponds to a Hopf bifurcation.Note that deterministic oscillations with a circadian period are limited only to a small strip near the stability boundary at the bottom right.The circle identifies the values of γ and KaiA that we obtain by fitting experimental power spectra in part (B).The diamond stands for best fit parameters obtained for Anabaena (Fig. 6C).(D) Comparison between damped deterministic oscillations (blue line) and quasi-cycles, both at the circle point outside the region of the deterministic oscillations in panel (C).
corresponding to oscillations characterized by circadian periods.This strongly suggests that the oscillatory behavior may correspond to noise-seeded oscillations as observed in other systems (McKane and Newman, 2005), without the need of training by light-dark cycles.The noise-seeded enlargement of the range of biological parameters over which circadian oscillatory behavior can be observed may confer robustness to clocks, a decisive biological advantage.Robustness may enable proper circadian clock function of Anabaena under a variety of stresses, including those involved in metabolic rewiring and the ensuing differentiation

Appendix1-figure 2 :
Stability diagram.Stability of the equilibrium points for the three species as a function of γ and [KaiA].The values of all the other paramters are specified in Appendix 1 -table 1. Continuous lines denote stable equilibrium points, dashed line unstable equilibrium points.Appendix1-figure3: Deterministic limit cycle region.The portion of the plane delimited by the blue line marks the region of the parameters space where deterministic regular oscillations occur.
Conjugation into Anabaena P pecB -gfp and selection for Nm R Anabaena P pecB -gfp 'kaiABC chromosome C.K1 (Nm R ) Selection for sucrose R and PCR test for gene replacement in independent exconjugants Gene replacement by double cross-over Appendix1-figure 7: Construction of an Anabaena ∆kaiABC deletion mutant in the PpecB-gfp genetic background.PCR, DNA restriction/ligation and transformation into E. coli by standard techniques.Conjugation from E. coli to Anabaena as described by

Table 1 :
To test whether the cell-cycle and circadian clocks are coupled in Anabaena cells, we monitored Synchronization of expression of a clock-controlled gene in cells within and between Anabaena filaments.