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

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Article and author information
Metrics

Abstract

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 et al., 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, have 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 cell-cell 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 hr 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 are expressed with a large amplitude. Nonetheless, about 80 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.

Results

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_{p e c B - g f p}$ at the level of individual cells along filaments. Prior to and during the experiments, filaments were grown under constant light conditions. Typical images of wild-type (WT) filaments at different time points are shown in Figure 1A (see also Figure 1, Video 1). One salient feature in the images is significant synchrony, that is, cellular oscillations progressed in individual cells along filaments together and with a similar period. The images in Figure 1A correspond to successive maxima and minima of the mean cell fluorescence intensity, µ, as a function of time, which we plot in Figure 2A. Similar experiments with a strain in which the kaiABC genes were deleted ( $Δ$ kaiABC) resulted in a low-level, non-oscillatory signal (Figure 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_{h e t R - g f p}$ did not exhibit an oscillatory behavior (Figure 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 Figure 1B, Figure 2A).

Figure 1

Download asset Open asset

Circadian oscillation in *Anabaena*.

(A) GFP fluorescence in a filament of an *Anabaena* strain bearing a $P_{p e c B - g f 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 (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 with 1 supplement see all

Download asset Open asset

Characterization of a clock-controlled gene in *Anabaena*.

(A) Average cell fluorescence intensity from $P_{p e c B - g f 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_{h e t R - g f p}$ (cyan line); and temporal dependence of the cell-cell variability $C V^{2}$ in expression of $P_{p e c B - g f 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_{p e c B - g f p}$ fusion in the lineages of two contiguous cells as a function of time. (D) Average spatial autocorrelation function of $P_{p e c B - g f 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 (A) and (C), corresponding to filaments of the *sepJ/fraCD* genetic background, see Figure 2—figure supplement 1.

Video 1

Download asset

posterframe for video — Real-time expression of a clock-controlled gene and filament autofluorescence during circadian oscillations in WT *Anabaena*.

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 \pm 0.04$ ) (Table 1). The value of $R$ for cells initially separated by intervals of 10 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.

Table 1

Synchronization of expression of a clock-controlled gene in cells within and between Anabaena filaments.

The synchronization index $R$ for strains with the indicated genotypes (Materials and methods) was measured from the fluorescence intensities of $P_{p e c B - g f 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 10 cells (separate). To measure synchronization between filaments, $R$ was computed from about 10 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 (Rust et al., 2007; Lambert et al., 2016; Dong et al., 2010; Teng et al., 2013), carried out in $n$ different experimental runs. Significance (p-value) in interstrain comparisons was established by the Mann–Whitney U-test, and * represents rejection of the null hypothesis that samples come from distributions with equal medians. WT: wild type.

Genotype	Cell cluster	R (mean ± SEM)	n	Comparison with strain	p-Value
WT	Contiguous	0.89 ± 0.04	3	WT (separate)	0.117
WT	Separate	0.85 ± 0.01	2
WT	Different filaments	0.75 ± 0.04	2	WT (separate)	0.026*
$Δ s e p J Δ f r a C D$	Contiguous	0.73 ± 0.05	4	WT (contiguous)	0.001*
$Δ k a i A B C$	Contiguous	0.71 ± 0.03	3	WT	0.001*

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 (Figure 2B). The expression from $P_{p e c B - g f 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 \pm 0.04$ , a value that is significantly smaller than that obtained for cells within the same filament (Figure 2C, Table 1). We surmise that initial synchronization may be due to phase resetting following the change in conditions, for example, illumination, upon transfer of cells from bulk culture to the microscope for real-time measurements. This change also could account for the decay in fluorescence intensity observed during the first circa 24 hr of our experiments (Figure 2A, B). Furthermore, this decay is $P_{p e c B - g f p}$ -specific, but clock-independent, as it was also observed in filaments of the $Δ k a i A B C$ background (Figure 2A).

Circadian clocks along filaments are coupled by cell-cell communication

Another salient feature in the snapshots of Figure 1A is the high spatial coherence of the expression from $P_{p e c B - g f p}$ along filaments, that is, 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 Materials and methods). We found that $g$ decays to zero for separations of five cells or more (Figure 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 (Figure 2D, magenta line) decreased to zero already at the second neighbor, suggesting that another factor, for example, 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_{p e c B - g f p}$ expression in filaments of a $Δ$ sepJ/ $Δ$ fraCD strain 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_{p e c B - g f p}$ in filaments of this strain were observed (Figure 2—figure supplement 1A). The resulting spatial autocorrelation function decreased over significantly shorter lengthscales (about two cells) than the WT (Figure 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 WT 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 WT filaments (Figure 2C). These findings indicate that the correlation in $P_{p e c B - g f 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 $Δ k a i A B C$ background, in which $P_{p e c B - g f 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 Figure 1A even though their amplitude is small relative to the clock-modulated activity of $P_{p e c B - g f p}$ . To quantify these variations, we calculated the square of the coefficient of variation $C V = σ / μ$ , where $σ$ denotes the standard deviation of the fluorescence intensity of contiguous cells along a filament. A plot of $C V^{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 (Figure 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). To test whether the cell cycle and circadian clocks are coupled in Anabaena cells, we monitored 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_{p e c B - g f p}$ fusion and their respective lineage are shown in Figure 2C. In Figure 2E, we show a histogram of the timing of cell division events as a function of the phase at which they occur along the circadian clock, obtained from traces similar to those in Figure 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). Figure 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 real-time quantitative polymerase chain reaction (RT-qPCR) measurements of WT 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 hr/12 hr 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 (Figure 3A). The oscillations were largely in phase, and transcription occurred 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 \pm 0.05$ , $R_{A, C} = 0.79 \pm 0.09$ , and $R_{B, C} = 0.87 \pm 0.06$ . Error bars were obtained by bootstrap methods (Materials and 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.

Figure 3 with 2 supplements see all

Download asset Open asset

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 (Materials and methods). A persistence homology analysis of these data is presented in Figure 3—figure supplement 1. (B, 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). Gray 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 (C) (Figure 3—figure supplement 2).

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 methods). 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 signal's period (Kennedy et al., 2018). We obtained $τ = 7.1 \pm 1.2$ hr, $τ = 6.7 \pm 1.2$ hr, and $τ = 7.3 \pm 1.1$ ( $n \geq 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, for example, pecB, we measured the relative expression of rpaA in WT and $Δ k a i A B C$ 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 (Figure 3B). Furthermore, pecB displays oscillatory behavior, with similar amplitude and phase. The oscillatory behavior of both genes was abolished in the $Δ$ kaiABC mutant (Figure 3B, 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 (Figure 4A). Circadian oscillations in the lineages of two sister cells are shown in Figure 4B (see also Figure 5B for the associated power spectrum). Expression from the kaiBC promoter exhibited circadian oscillations with a period of about 25 hr, 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 Figure 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 Equation S4 in Rust et al., 2007) is modeled here by a continuous function (Appendix 1 – Equation 6) that makes analytical progress possible. The parameter $γ$ in our nonlinear phosphorylation and dephosphorylation rates $k_{X Y}$ 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 concentrations of S ( $ϕ_{S} = [S-KaiC]$ ), T ( $ϕ_{T} = [T-KaiC]$ ), and D ( $ϕ_{D} = [D-KaiC]$ ). The analysis of these equations allows us to determine the region in parameter space within which the model exhibits sustained deterministic oscillations (Figure 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 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 Figure 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 system outside the region of deterministic oscillations (circle in Figure 5C), suggesting that circadian rhythms can be a manifestation of noise-driven oscillations (Figure 5D). Remarkably, the fitted value for $[KaiA] = 1.308$ µM matches the concentration reported previously ( $[KaiA] = 1.3$ µM; see Rust et al., 2007).

Figure 4

Download asset Open asset

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 with 1 supplement see all

Download asset Open asset

Stochastic model for circadian oscillations in *Synechococcus*.

(A) Schematic representation of interconversion between KaiC phosphoforms modulated by the activity of KaiA in an individual clock. The different phosphoform 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 phosphoforms $X, Y$ with the indicated rates $k_{x y}$ . 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 (B). The diamond stands for best fit parameters obtained for *Anabaena* (Figure 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 (C).

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; Figure 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 that 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 Figure 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 (Figure 2C, Figure 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 similarly able to 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 Figure 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 (Figure 6C) upon adjusting the two fitting parameters, $γ$ and [KaiA]. Again, the fitted values position the system (diamond in Figure 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 Figure 5C. The bifurcation line that 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.

Figure 6

Download asset Open asset

Stochastic model for circadian oscillations in *Anabaena*.

(A) Schematic representation of the *Anabaena* filament showing coupling of circadian clocks via cell-cell communication (red arrows). (B) Gillespie simulations of quasi-cycles of T-KaiC in a continuous stretch of 10 cells along a filament. The reaction parameters correspond to the diamond plotted in Figure 5C. The total amount of KaiC phosphoforms was set to 5000, and the number of steps of the Gillespie algorithm was $1.6 \times 10^{7}$ . (C) Average power spectrum of single-cell fluorescence intensity along filaments (red symbols) fit to the data with the prediction from the stochastic model (blue line). The best fit values correspond to the diamond shown in Figure 5C. (D) Complex coherence function measuring the correlation expression from $P_{p e c B - g f p}$ in 35 cell segments at the frequency of temporal oscillations. Red full circles correspond to experimental data, blue squares represent the fit to the experimental data with the prediction of the stochastic model, and empty red circles represent the coherence function of the experimental data in which cells have been reshuffled. 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 Figure 2D. The intercell coupling was obtained from the fit in (C).

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_{p e c B - g f 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 Figure 6D, the coherence at the characteristic frequency of oscillations was a monotonically decaying function of distance, an observation that points to the existence of non-local interactions by cell-cell communication along the filament. As complementary evidence, we notice that the magnitude of the coherence function took a constant value when spatial correlations were broken by randomly reshuffling the position of the cells along the filament. Remarkably, the prediction of the stochastic model captures the experimental behavior, as shown in Figure 6C. The range of interaction of the exponential kernel is self-consistently quantified by the fitting strategy and yields a result that agrees with that obtained from the spatial autocorrelation analysis (Figure 2D). More specifically, this conclusion is consistent with the reduction of spatial autocorrelation in the expression of a clock-regulated gene when either cell-cell communication is perturbed or clock genes are deleted, and with the notion that inheritance following cell division alone cannot account for all the spatial variation of expression along filaments (Figure 2D). In Figure 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 protein 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 (Figure 2D) and the complex coherence function (Figure 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 state, 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_{p e c B - g f 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 decentralized 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_{p e c B - g f p}$ indicated that the cell-cell variability along filaments as measured by $C V^{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 $C V^{2}$ in $P_{p e c B - g f 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, that is, 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 Figure 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) suggests 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 are small, but that of pecB and many other targets are 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 are 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 mechanism behind 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 (Figure 5B, Figure 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 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 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 and 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 suggest 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.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Strain, strain background (Anabaena)	$P_{p e c B - g f p}$ , WT	This paper		Anabaena PCC 7120 WT, bearing a pecB promoter fusion to gfp
Strain, strain background (Anabaena)	$P_{p e c B - g f p}$ , $Δ k a i A B C$	This paper		Anabaena PCC 7120 deletion mutant of the $k a i A B C$ genes, bearing a pecB promoter fusion to gfp
Strain, strain background (Anabaena)	$P_{h e t R - g f p}$	doi: 10.1371/journal.pgen.1005031	CSL64	Anabaena PCC 7120 WT, bearing a hetR promoter fusion to gfp
Strain, strain background (Anabaena)	$P_{p e c B - g f p}$ , $Δ s e p J / Δ f r a C / Δ f r a D$	This paper		Anabaena PCC 7120 deletion mutant of the $s e p J$ , $f r a C$ , $f r a D$ genes (CSVM141), bearing a pecB promoter fusion to gfp
Strain, strain background (Synechococcus elongatus)	YFP-SsrA	This paper	PCC 7942	Synechococcus elongatus PCC 7942 (wild-type) expressing YFP- SsrA
Recombinant DNA reagent	EB2316 (plasmid)	Addgene plasmid	87753	http://n2t.net/addgene: 87753
Recombinant DNA reagent	pSpark (plasmid)	Canvax	C0001	https://lifescience.canvaxbiotech.com/wpcontent/uploads/sites/2/2015/08pSpark-DNA-Cloning.pdf
Recombinant DNA reagent	pCSRO $s a c B -$ containing cloning vector	doi: 10.1128/JB.00181-13
Commercial assay or kit	Fast SYBR Green Master Mix	Applied Biosystems	4385612
Commercial assay or kit	QuantiTect Reverse Transcription kit	QIAGEN	205311

Strains

Strains bearing a chromosomally encoded $P_{p e c B - g f p}$ , were obtained by conjugation with the following backgrounds: WT Anabaena sp. (also known as Nostoc sp.) strain PCC 7120; $Δ$ sepJ/ $Δ$ fraCD, 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. Strain CSL64 bearing a chromosomally encoded $P_{h e t R - g f p}$ transcriptional fusion in a WT background has been reported previously (Corrales-Guerrero et al., 2015). The S. elongatus 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 (Addgene plasmid # 87753; http://n2t.net/addgene:87753; RRID: Addgene 87753).

Culture conditions

Request a detailed protocol

Strains and derived strains were grown photoautotrophically in BG11 medium containing NaNO₃, 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 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 hr 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

Request a detailed protocol

Strains were grown as described previously (Di Patti et al., 2018). When required, antibiotics, Sm and 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 for about 5 days), were adjusted so as to have a chlorophyll a content of 2–4 µg/mL, 24 hr 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₃ and 10 mM NaHCO₃, 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) and tungsten halogen light (10 µmol, 3000K color). 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 10 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 a 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 OD750 of 0.3–0.4. Samples were then entrained by two light-dark cycles (12 hr–12 hr) before measurements commenced. The cultures were then diluted to a OD750 of 0.1 using BG11 medium, and 2 µL of the culture was pipetted onto a glass-bottom culture dish. A patterned pad prepared as above but solidified on an optical grating (Hadizadeh Yazdi 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

Request a detailed protocol

For 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 hr–12 hr) and then grown under constant light. Then, measurements were started 24 hr after, and total RNA was extracted from the PCC 7120 and $Δ$ kaiABC cultures in two biological replicates, every 4 hr. 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 min 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 dye and specific primers (Appendix 1—table 2). The amplification protocol was one cycle at 95°C for 10 min, 40 cycles of 95°C for 15 s, and 60°C for 60 s. 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 3–4 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^{- Δ Δ c t}$ (Pfaffl, 2001).

Image segmentation

Request a detailed protocol

All image processing and data analysis were 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, was obtained as output for further statistical analysis.

Analysis of synchronization and spatial correlations along filaments

Request a detailed protocol

Synchronization was measured by the order parameter proposed by Garcia-Ojalvo et al., 2004:

R = \frac{⟨ μ^{2} ⟩ - ⟨ μ ⟩^{2}}{\bar{⟨ f_{i}^{2} ⟩ - ⟨ f_{i} ⟩^{2}}}

where $⟨ \cdot ⟩$ denotes a time average, $\bar{.}$ 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 $Δ k a i A B C$ background, for which $R$ was calculated for an interval of 24 hr, 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 U-test.

The spatial autocorrelation function of fluorescence intensities along filaments was 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 $Δ k a i A B C$ background, autocorrelation functions were calculated at about 30 and 50 hr, corresponding to the first minimum and maximum of oscillations observed in filaments of WT background. The spatial autocorrelation function was calculated from the formula

g_{k} = \frac{c_{k}}{c_{0}}

where

c_{k} = \frac{1}{N} \sum_{i = 1}^{N - k} (f_{i} - \bar{f}) (f_{i + k} - \bar{f})

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

Request a detailed protocol

Each kinetic parameter ( $k$ ) was randomly selected from a Gaussian distribution centered at the nominal value ( $\bar{k}$ ) estimated by Rust et al., 2007 using the normrnd MATLAB function. The standard deviation is assigned as $σ = \bar{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 \pm 1.4$ and [KaiA] $= 1.3 \pm 0.24$ (mean ± SD).

Appendix 1

Persistent homology analysis of periodic behavior

To apply persistent homology methods for the detection of oscillatory behavior in the noisy low-amplitude time series of kaiA, kaiB, and kaiC, a two-dimensional representation of each time series was constructed by a time-delayed embedding, that is, a plot of the signal versus itself at two different time points separated by a delay $τ$ (see Figure 3—figure supplement 1). To overcome the sparsity of data (12 time points) when plotting the time-delayed embedding, each time series was first padded by intermediate points determined by linear interpolation. The delay $τ$ was then evaluated as the first minimum of the mutual information of the series, as determined by a fifth degree polynomial fit, following previously published procedures (Fraser and Swinney, 1986). The existence of an underlying oscillation in the data was assessed by using a MATLAB algorithm based on the JavaPlex package (Tausz et al., 2014), which yields as output a barcode diagram in which each bar represents the range of scales over which a cycle persists as observed in the phase portrait, regarded as a simplicial complex. Note that for a periodic signal the delay $τ$ corresponds to a quarter of the signal's period (Kennedy et al., 2018).

Stochastic model

We model the Anabaena filament as a one-dimensional chain, where each node represents a cell with its individual circadian clock. We use the letter $α$ to label the node and denote by $Ω$ the length of the chain, so that $α = 1, \dots, Ω$ (see Appendix 1—figure 1). As explained in the main text, the coupling between clocks that result from metabolic factors yields a long-range interaction between clocks along the filament. To account for this crucial phenomenon, we schematize the chain as a network and, therefore, introduce the adjacency matrix $W$ whose elements $W_{α β}$ are 1 if nodes $α$ and $β$ are connected, 0 otherwise. The connectivity of each node is assumed to be constant for each node and is denoted by $k_{α}$ .

Appendix 1—figure 1

Download asset Open asset

Schematic representation of an *Anabaena* filament.

The parameter $k_{α}$ measures the connectivity of each node.

Referring back to Figure 3A, we consider, in each cell, the cyclic interconversion between phosphoforms of KaiC, T, S, D (ST), and U encapsulated in the following chemical equations:

\begin{array}{ll} U^{α} \overset{k_{U T}}{\to} T^{α} T^{α} \overset{k_{T U}}{\to} U^{α} D^{α} \overset{k_{D T}}{\to} T^{α} S^{α} \overset{k_{S U}}{\to} U^{α} \\ U^{α} \overset{k_{U S}}{\to} S^{α} T^{α} \overset{k_{T D}}{\to} D^{α} D^{α} \overset{k_{D S}}{\to} S^{α} S^{α} \overset{k_{S D}}{\to} D^{α} \end{array}

where the superscript $α$ indicates that the species belongs to the cell $α$ . The interconversion rates between phosphoforms are given by

k_{X Y} (\frac{n_{S}^{α}}{N}) = k_{X Y}^{0} + k_{X Y}^{A} f (\frac{n_{S}^{α}}{N}) for (X, Y) \in {U, T, D, S}

with

f (\frac{n_{S}^{α}}{N}) = \frac{Γ}{1 + \exp [\frac{γ}{[KaiA]} (\frac{n_{S}^{α} [KaiC]}{N} - \frac{[KaiA]}{2})]}

where [KaiA] (resp. [KaiC]) is the concentration of KaiA (resp. [KaiC]). The values of the rates are set as those of Lambert et al., 2016 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 it interacts with S-KaiC instantaneously (Lambert et al., 2016). The typical sigmoidal shape of the nonlinear function (Teng et al., 2013) is plotted in Figure 5—figure supplement 1 (blue dashed line), where we also plotted the corresponding function used in the paper by Lambert et al., 2016 (red line).

Appendix 1—table 1

Parameters used in the simulations.

$k_{U T}^{0}$	0 h⁻¹	$k_{U T}^{A}$	0.479077 h⁻¹
$k_{T D}^{0}$	0 h⁻¹	$k_{T D}^{A}$	0.212923 h⁻¹
$k_{S D}^{0}$	0 h⁻¹	$k_{S D}^{A}$	0.505692 h⁻¹
$k_{U S}^{0}$	0 h⁻¹	$k_{U S}^{A}$	0.0532308 h⁻¹
$k_{T U}^{0}$	0.21 h⁻¹	$k_{T U}^{A}$	0.0798462 h⁻¹
$k_{D T}^{0}$	0 h⁻¹	$k_{D T}^{A}$	0.1730000 h⁻¹
$k_{D S}^{0}$	0.31 h⁻¹	$k_{D S}^{A}$	−0.319885 h⁻¹
$k_{S U}^{0}$	0.11 h⁻¹	$k_{S U}^{A}$	−0.133077 h⁻¹
$k_{1 / 2}$	0.43 µM	[KaiC]	3.4 µM

All the values are from Lambert et al., 2016.

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 (see 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 describe this effect through long-range 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

T^{α} \overset{ν}{\to} T^{β} D^{α} \overset{ν}{\to} D^{β} S^{α} \overset{ν}{\to} S^{β}

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 $N = n_{U}^{α} + n_{S}^{α} + n_{T}^{α} + n_{D}^{α}$ . For this reason, at each time $t$ the state of the system is specified by the $3 Ω$ -dimensional vector $𝒏 (t)$

n (t) = (n_{T}^{1} (t), n_{D}^{1} (t), n_{S}^{1} (t); \dots; n_{T}^{Ω}, n_{D}^{Ω}, n_{S}^{Ω}) = (n^{1} (t), \dots, n^{Ω} (t))

The microscopic dynamics described by rules (Equations 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 $𝕋 (𝒏^{'} | 𝒏)$ from state $𝒏$ to a new state $𝒏^{'}$ for all the chemical equations that constitute the stochastic model. To denote the new state $𝒏^{'}$ , we only explicitly write the components of the vector that change. Accordingly, the phosphorylation and dephosphorylation transition rates in each cell $α$ read

T (n_{i}^{α} + 1 | n) = \frac{n_{U}^{α}}{N} k_{U i}^{(α)} i \in {T, S}

T (n_{i}^{α} - 1 | n) = \frac{n_{i}^{α}}{N} k_{i U}^{(α)} i \in {T, S}

while the interconversion transitions are

T (n_{i}^{α} - 1, n_{j}^{α} + 1 | n) = \frac{n_{i}^{α}}{N} k_{i j}^{(α)} (i, j) \in {(T,D), (D,T), (D,S), (S,D)}

To account for cell-cell communication, we also add the transitions for long-range interaction between two connected cells $α$ and $β$ :

T (n_{i}^{α} - 1, n_{i}^{β} + 1 | n) = \frac{n_{i}^{α}}{N} ν i \in {T, D, S}

The dynamics is fully described in terms of the following master equation:

\frac{\partial P (𝒏, t)}{\partial t} = \sum_{𝒏^{'} \neq 𝒏} [𝕋 (𝒏 | 𝒏^{'}) P (𝒏^{'}, t) - 𝕋 (𝒏^{'} | 𝒏) P (𝒏, t)]

whose explicit version reads

\begin{array}{ll} \frac{\partial P (n, t)}{\partial t} & = \sum_{α = 1}^{Ω} {T (n | n_{T}^{α} - 1) P (n_{T}^{α} - 1; t) - T (n_{T}^{α} - 1 | n) P (n, t) + \\ + T (n | n_{T}^{α} - 1, n_{D}^{α} + 1) P (n_{T}^{α} - 1, n_{D}^{α} + 1; t) - T (n_{T}^{α} - 1, n_{D}^{α} + 1 | n) P (n, t) + \\ + T (n | n_{T}^{α} + 1, n_{D}^{α} - 1) P (n_{T}^{α} + 1, n_{D}^{α} - 1; t) - T (n_{T}^{α} + 1, n_{D}^{α} - 1 | n) P (n, t) + \\ + T (n | n_{D}^{α} - 1, n_{S}^{α} + 1) P (n_{D} - 1^{α}, n_{S}^{α} + 1; t) - T (n_{D}^{α} - 1, n_{S}^{α} + 1 | n) P (n, t) + \\ + T (n | n_{S}^{α} - 1) P (n_{S}^{α} - 1; t) - T (n_{S}^{α} - 1 | n) P (n, t) + \\ + T (n | n_{D}^{α} + 1, n_{S}^{α} - 1) P (n_{D}^{α} + 1, n_{S}^{α} - 1; t) - T (n_{D}^{α} + 1, n_{S}^{α} - 1 | n) P (n, t) + \\ + T (n | n_{T}^{α} + 1) P (n_{T}^{α} + 1; t) - T (n_{T}^{α} + 1 | n) P (n, t) + \\ + T (n | n_{S}^{α} + 1) P (n_{S}^{α} + 1; t) - T (n_{S}^{α} + 1 | n) P (n, t) + \end{array}

\begin{array}{ll} + \sum_{β = 1}^{Ω} W_{α β} [T (n | n_{T}^{α} + 1, n_{T}^{β} - 1) P (n_{T}^{α} + 1, n_{T}^{β} - 1; t) - T (n_{T}^{α} + 1, n_{T}^{β} - 1 | n) P (n, t) + \\ + T (n | n_{D}^{α} + 1, n_{D}^{β} - 1) P (n_{D}^{α} + 1, n_{D}^{β} - 1; t) - T (n_{D}^{α} + 1, n_{D}^{β} - 1 | n) P (n, t) + \\ + T (n | n_{S}^{α} + 1, n_{S}^{β} - 1) P (n_{S}^{α} + 1, n_{S}^{β} - 1; t) - T (n_{T}^{α} + 1, n_{T}^{β} - 1 | n) P (n, t) + \\ + T (n | n_{T}^{α} - 1, n_{T}^{β} + 1) P (n_{T}^{α} - 1, n_{T}^{β} + 1; t) - T (n_{T}^{α} - 1, n_{T}^{β} + 1 | n) P (n, t) + \\ + T (n | n_{D}^{α} - 1, n_{D}^{β} + 1) P (n_{D}^{α} - 1, n_{D}^{β} + 1; t) - T (n_{D}^{α} - 1, n_{D}^{β} + 1 | n) P (n, t) + \\ + T (n | n_{S}^{α} - 1, n_{S}^{β} + 1) P (n_{S}^{α} - 1, n_{S}^{β} + 1; t) - T (n_{T}^{α} - 1, n_{T}^{β} + 1 | n) P (n, t)]} \end{array}

The master equation can be written in a more compact form though the use of the step operators $ϵ_{i}^{\pm (α)}$ that act on a generic function $g (𝒏)$ as $ϵ_{i}^{\pm (α)} g (n_{i}^{α}) = g (n_{i}^{α} \pm 1) \forall i \in {T, D, S}$ and $α = 1, \dots, Ω$ :

\begin{array}{ll} \frac{\partial P (n, t)}{\partial t} = & \sum_{α = 1}^{Ω} {[(ϵ_{T}^{- (α)} - 1) T (n_{T}^{α} + 1 | n) + (ϵ_{S}^{- (α)} - 1) T (n_{S}^{α} + 1 | n) + \\ + (ϵ_{T}^{+ (α)} - 1) T (n_{T}^{α} - 1 | n) + (ϵ_{T}^{+ (α)} ϵ_{D}^{- (α)} - 1) T (n_{T}^{α} - 1, n_{D}^{α} + 1 | n) + \\ + (ϵ_{T}^{- (α)} ϵ_{D}^{+ (α)} - 1) T (n_{T}^{α} + 1, n_{D}^{α} - 1 | n) + (ϵ_{D}^{+ (α)} ϵ_{S}^{- (α)} - 1) T (n_{D}^{α} - 1, n_{S}^{α} + 1 | n) + \\ + (ϵ_{S}^{+ (α)} - 1) T (n_{S}^{α} - 1 | n) + (ϵ_{D}^{- (α)} ϵ_{S}^{+ (α)} - 1) T (n_{D}^{α} + 1, n_{S}^{α} - 1 | n)] P (n, t) + \\ + \sum_{β = 1}^{Ω} W_{α β} [(ϵ_{T}^{+ (α)} ϵ_{T}^{- (β)} - 1) T (n_{T}^{α} - 1, n_{T}^{β} + 1 | n) + \\ + (ϵ_{D}^{+ (α)} ϵ_{D}^{- (β)} - 1) T (n_{D}^{α} - 1, n_{D}^{β} + 1 | n) + \\ + (ϵ_{S}^{+ (α)} ϵ_{S}^{- (β)} - 1) T (n_{S}^{α} - 1, n_{S}^{β} + 1 | n)] P (n, t)} \end{array}

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, which is obtained taking the limit $N \to \infty$ denoted by $ϕ_{X}^{(α)}$ , and the fluctuations $ξ_{X}^{(α)}$ whose amplitude scales as $1 / \sqrt{N}$ , as

\begin{array}{ll} \frac{n_{T}^{α} (t)}{N} = ϕ_{T}^{α} (t) + \frac{ξ_{T}^{α}}{\sqrt{N}} \\ \frac{n_{D}^{α} (t)}{N} = ϕ_{D}^{α} (t) + \frac{ξ_{D}^{α}}{\sqrt{N}} α = 1, \dots, Ω \\ \frac{n_{S}^{α} (t)}{N} = ϕ_{S}^{α} (t) + \frac{ξ_{S}^{α}}{\sqrt{N}} \end{array}

We introduce a new probability function $Π (𝝃, t) \equiv P (𝒏 (ϕ (t), 𝝃); t)$ , where the $𝝃$ is a $3 Ω$ -dimensional vector whose components are

ξ = (ξ_{T}^{1}, ξ_{D}^{1}, ξ_{S}^{1}, \dots, ξ_{T}^{Ω}, ξ_{D}^{Ω}, ξ_{S}^{Ω}) \equiv (ξ^{1}, \dots, ξ^{Ω})

so that the left-hand side of the master Equation (15) can be recast in the form

\begin{array}{ll} \frac{\partial P}{\partial t} & = \sum_{α = 1}^{Ω} [\frac{1}{N} \frac{\partial Π}{\partial τ} - \frac{1}{\sqrt{N}} \sum_{i \in {T, D, S}} \frac{\partial Π}{\partial ξ_{i}^{α}} \frac{d ϕ_{i}^{α}}{d τ}] \end{array}

with $τ = t / N$ .

Operating an expansion of the step operators with respect to the small parameter $1 / \sqrt{N}$ , performing the change of variable (Equation 16) into the right-hand side of the master equation, and collecting terms proportional to $1 / \sqrt{N}$ , we get

\frac{d ϕ_{i}^{α}}{d τ} = h_{i}^{(α)} (ϕ^{α}) + ν \sum_{β = 1}^{Ω} Δ_{α β} ϕ_{i}^{β}

where

h_{T}^{(α)} (ϕ^{α}) = ϕ_{U}^{α} k_{U T} + ϕ_{D}^{α} k_{D T} - ϕ_{T}^{α} k_{T U} - ϕ_{T}^{α} k_{T D}

h_{D}^{(α)} (ϕ^{α}) = ϕ_{T}^{α} k_{T D} + ϕ_{S}^{α} k_{S D} - ϕ_{D}^{α} k_{D T} - ϕ_{D}^{α} k_{D S}

h_{S}^{(α)} (ϕ^{α}) = ϕ_{U}^{α} k_{U S} + ϕ_{D}^{α} k_{D S} - ϕ_{S}^{α} k_{S U} - ϕ_{S}^{α} k_{S D}

and the matrix $𝚫$ denotes the discrete Laplacian operator. The concentration $ϕ_{U}$ reads $[KaiC] - ϕ_{T} - ϕ_{D} - ϕ_{S}$ and $[KaiC] ϕ_{X} \to ϕ_{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 (Equation 19) simplifies to

\frac{d ϕ_{i}}{d τ} = h_{i} (ϕ^{α})

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 (Equation 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 $J$

J = (\begin{matrix} - k_{U T} - k_{T D} - k_{T U} & - k_{U T} + k_{D T} & J 1 \\ k_{T D} & - k_{D T} - k_{D S} & J 2 \\ - k_{U S} & - k_{U S} + k_{D S} & J 3 \end{matrix})

where

\begin{array}{ll} J 1 & = - k_{U T} + f^{'} (ϕ_{S}^{*}) [ϕ_{U}^{*} k_{U T}^{A} + ϕ_{T}^{*} (- k_{T U}^{A} - k_{T D}^{A}) + ϕ_{D}^{*} k_{D T}^{A}] \\ J 2 & = k_{S D} + f^{'} (ϕ_{S}^{*}) [ϕ_{T}^{*} k_{T D}^{A} + ϕ_{D}^{*} (- k_{D T}^{A} - k_{D S}^{A}) + ϕ_{S}^{*} k_{S D}^{A}] \\ J 3 & = - k_{U S} - k_{S D} - k_{S U} + f^{'} (ϕ_{S}^{*}) [ϕ_{U}^{*} k_{U S}^{A} + ϕ_{D}^{*} k_{D S}^{A} + ϕ_{S}^{*} (- k_{S U}^{A} - k_{S D}^{A})] \end{array}

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 denote 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 of 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).

Appendix 1—table 2

Oligodeoxynucleotide primers used in this work.

Name	Sequence 5′–3′
alr0523-EcoRI-Fw	TTTTGAATTCGCTTATAAACAGCAGTTAACAGGCT
alr0523-Rev	TGCTACCTCCACCGCCTGCCTGTTCAACTACTTTGGA
4G-GFP-Fw	GCGGTGGAGGTAGCAAAGGAGAAGAACTTTTCAC
GFP-Rev	GCCTGAATTCTTATTTGTATAGTTCATCCATGCC
alr0523-1-Fw	GAATTCGCTTATAAACAGCAGTTAACAGGCT
alr0523-1-Rev	CTAGCACCTCCACCGCCTGCCTGTTCAACTACTTTGGA
4G-GFP-Fw in plasmid PCSV3	ATTTGAAACTGCGCCACGGATC
Rev plasmid PCSV3	GACCATGACGGATTAGCTCAGTAG
alr0523(7120)–1	CGT GAG TCT CCA ACG GAG GC
Kai-1	GAAACTGCAGGCAGAATAGGAAATCTCTAC
Kai-2	CCAAATGATATCGTGCTGACAAACCTACAGTGC
Kai-3	CAGCACGATATCATTTGGTATCGTACTATATTC
Kai-4	CTTTCTGCAGGTTGTCCAGCCAGCAGGGTAG
kaiA-2	CAGGGTGAGGCGATAATCCAT
kaiA-1	GCCAGAGTACTTGTTTCTAAGCAAC
CK1-R	CGATTCCGAAGCCCAACCT
kaiC-4	CGAGCTACCAACCGAAAG
kaiB-1	CGGCAATACTCCAAACTCAG
PkaiBC-1	GGTCTATCCCACGAGAAACC
YFP-2	GGTAGCTCAGGTAGTGGTTGTC
all5167-1q (forward)	GCTCAAGCAATTCGTCACTGTTCC
all5167-2q (reverse)	AAAGATTGCGTCGGTCTGGTGT
rnpB-1q (forward)	CTCTTGGTAAGGGTGCAAAGGTG
rnpB-4q (reverse)	GGCTCTCTGATAGCGGAACTGG
kaiC-3q (forward)	ATGAAGCAGTGGGAGTGGTG
kaiC-8q (reverse)	ACGTTACGGGCTATGACCAC
kaiB-1q (forward)	ACCAAATTCAGTCAGGGCGT
kaiB-4q (reverse)	GCCAATCAGAACTCTTTCCCG
kaiA-3q (forward)	CAACTCAAATCAGATTATCGCCA
kaiA-4q (reverse)	CTGCCCTCTAGTCGTAGCTG
pecB-1q (forward)	ATATTTAATGCTGGTGGTGCTTGTT
pecB-4q (reverse)	GCAGCGATCGTCCATGACACTAC
rpaA-1q (forward)	TTTAACGCCGGAGCAGATGA
rpaA-4q (reverse)	TGTCCGTGACGTTGTAGCAA

Appendix 1—figure 2

Download asset Open asset

Stability diagram.

Stability of the equilibrium points for the three species as a function of $γ$ and $[K a i A]$ . The values of all the other parameters are specified in Appendix 1—table 1. Continuous lines denote stable equilibrium points, and dashed lines denote unstable equilibrium points.

To better investigate the region of the parameter's 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 Figure 3B. From Figure 3B, it appears that the region with period $T ≃ 24$ hr 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 (Equation 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.

Appendix 1—figure 3

Download asset Open asset

Deterministic limit cycle region.

The portion of the plane delimited by the blue line marks the region of the parameter's space where deterministic regular oscillations occur.

Appendix 1—figure 4

Download asset Open asset

Deterministic simulation.

Results of the numerical integration of system (Equation 23) for $γ = 8$ and $[K a i A] = 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].

Appendix 1—figure 5

Download asset Open asset

Comparison between deterministic and stochastic simulations.

The blue lines denote the result of the numerical integration of Equation (23) while the noisy red lines represent the stochastic simulation of the system through the Gillespie, 1977 algorithm. For all the plots for $γ = 8$ and $[K a i A] = 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].

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 $Π (𝝃, τ)$ , which results in the following Fokker–Planck equation:

\frac{\partial Π (𝝃, t)}{\partial τ} = \sum_{α = 1}^{Ω} {- \sum_{i = 1}^{3} \frac{\partial}{\partial ξ_{i}^{α}} [𝒜_{i}^{α} (𝝃) Π (𝝃, τ)] + \frac{1}{2} \sum_{i, j = 1}^{3} \sum_{β = 1}^{Ω} \frac{\partial^{2}}{\partial ξ_{i}^{α} \partial ξ_{j}^{β}} [ℬ_{i j, α β} Π (𝝃, τ)]}

where $𝒜$ is

𝒜_{i}^{α} = \sum_{j = 1}^{3} \sum_{β = 1}^{Ω} ℳ_{i j, α β} ξ_{j}^{β}

The two matrices $ℳ$ and $ℬ$ have dimension $3 Ω \times 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 $(S P)$ , and the other 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

ℳ_{i j, α β} = ℳ_{i j}^{(N S)} δ_{α β} + ℳ_{i j}^{(S P)} Δ_{α β}

with

ℳ^{(N S)} = J

ℳ^{(S P)} = (\begin{matrix} ν & 0 & 0 \\ 0 & ν & 0 \\ 0 & 0 & ν \end{matrix})

where $J$ is the Jacobian matrix relative to system (Equation 19) evaluated at the steady state. For $ℬ$ , we have

ℬ_{i j, α β} = ℬ_{i j}^{(N S)} δ_{α β} + ℬ_{i j}^{(S P)} Δ_{α β}

with

ℬ^{(N S)} \equiv (\begin{matrix} ϕ_{U}^{*} k_{U T} + ϕ_{D}^{*} k_{D T} + & - ϕ_{T}^{*} k_{T D} - ϕ_{D}^{*} k_{D T} & 0 \\ + ϕ_{T}^{*} k_{T U} + ϕ_{T}^{*} k_{T D} \\ - ϕ_{T}^{*} k_{T D} - ϕ_{D}^{*} k_{D T} & ϕ_{T}^{*} k_{T D} + ϕ_{S}^{*} k_{S D} + & - ϕ_{D}^{*} k_{D S} - ϕ_{S}^{*} k_{S D} \\ + ϕ_{D}^{*} k_{D T} + ϕ_{D}^{*} k_{D S} \\ 0 & - ϕ_{D}^{*} k_{D S} - ϕ_{S}^{*} k_{S D} & ϕ_{U}^{*} k_{U S} + ϕ_{D}^{*} k_{D S} + \\ + ϕ_{S}^{*} k_{S U} + ϕ_{S}^{*} k_{S D} \end{matrix})

ℬ^{(S P)} = (\begin{matrix} - 2 ν ϕ_{T}^{*} & 0 & 0 \\ 0 & - 2 ν ϕ_{D}^{*} & 0 \\ 0 & 0 & - 2 ν ϕ_{S}^{*} \end{matrix})

Spatial correlations

The Fokker–Planck Equation (25) describes the fluctuations about the equilibrium point $ϕ^{*}$ , and it is equivalent to the following coupled Langevin equations:

\frac{d ξ_{i}^{α}}{d τ} = \sum_{j = 1}^{3} \sum_{β = 1}^{Ω} ℳ_{i j, α β} ξ_{j}^{β} + η_{i}^{α}

where $ξ_{i}^{α}$ is the component of vector $𝝃^{α}$ defined in Equation (17) while $η_{i}^{α}$ is the ith component of the noise $𝜼^{α}$ with the following properties:

⟨ η_{i}^{α} (τ) ⟩ = 0

⟨ η_{i}^{α} (τ) η_{j}^{β} (τ^{'}) ⟩ = ℬ_{i j, α β} δ (τ - τ^{'})

Starting from 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., 2019; Zankoc et al., 2017). The (PSDM) $P_{i j, α β} (ω)$ is defined by

P_{i j, α β} (ω) = ⟨ {\tilde{ξ_{i}}}^{α} (ω) {\tilde{ξ_{j}}}^{β *} (ω) ⟩

where ${\tilde{ξ}}_{i}^{α} (ω)$ is the temporal Fourier transform of the $i$ th 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 Ω \times 3 Ω$ matrix $C$ , whose entries are the ratio of the elements of PSDM to a normalization coefficient:

C_{i j, α β} (ω) = \frac{P_{i j, α β}}{\sqrt{P_{i i, α α} P_{j j, β β}}}

To obtain an analytic expression of $P_{i j, α β}$ , 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

{\tilde{ξ}}_{i}^{α} (ω) = \sum_{j = 1}^{3} \sum_{β = 1}^{Ω} {(- i ω δ_{i j, α β} - ℳ_{i j, α β})}^{- 1} η_{j}^{β} (ω) = \sum_{j = 1}^{3} \sum_{β = 1}^{Ω} Φ_{i j, α β}^{- 1} (ω) η_{j}^{β} (ω)

with $Φ (ω) = - i ω 𝟙 - ℳ$ , where $𝟙$ denotes the $3 Ω \times 3 Ω$ identity matrix. Substituting Equation (36) into Equation (34), one eventually obtains

\begin{array}{ll} P_{i j, α β} (ω) & = ⟨ {\tilde{ξ_{i}}}^{α} (ω) {\tilde{ξ_{j}}}^{β *} (ω) ⟩ = \\ = \sum_{l, m = 1}^{3} \sum_{γ, δ = 1}^{Ω} Φ_{i l, α γ}^{- 1} (ω) ⟨ {\tilde{η}}_{l}^{γ} (ω) {\tilde{η}}_{m}^{* δ} (ω) ⟩ (Φ_{m j, δ β}^{†})^{- 1} (ω) = \\ = \sum_{l, m = 1}^{3} \sum_{γ, δ = 1}^{Ω} Φ_{i l, α γ}^{- 1} (ω) (ℬ_{l m, γ δ}) (Φ_{m j, δ β}^{†})^{- 1} (ω) = \\ = (Φ^{- 1} ℬ (Φ^{†})^{- 1})_{i j, α β} \end{array}

where $Φ^{†} = {(Φ^{*})}^{T}$ is the complex conjugate transpose of $Φ$ and use has been made of the following properties:

⟨ {\tilde{η}}_{i}^{α} (ω) ⟩ = 0

⟨ {\tilde{η}}_{i}^{α} (ω) {\tilde{η}}_{j}^{β} (ω^{'}) ⟩ = ℬ_{i j, α β}

Accounting for a smooth kernel of interactions

In this section, we will extend the analysis so as to account for 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 that are far apart from each other. The newly imposed modulation reflects in the transition rates (12) as follows:

T (n_{i}^{α} - 1, n_{i}^{β} + 1 | 𝒏) = \frac{n_{i}}{N} ν f (| α - β |)

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

W_{α β}^{'} = W_{α β} f (| α - β |)

and the corresponding Laplacian matrix $Δ^{'}$ now becomes

Δ_{α β}^{'} = W_{α β}^{'} - k_{α} δ_{α β}

Accordingly, the two matrices (Equations 27 and 29) appearing in the Fokker–Planck Equation (25) must be replaced by

ℳ_{i j, α β} = ℳ_{i j}^{(N S)} δ_{α β} + ℳ_{i j}^{(S P)} Δ_{α β}^{'}

ℬ_{i j, α β} = ℬ_{i j}^{(N S)} δ_{α β} + ℬ_{i j}^{(S P)} Δ_{α β}^{'}

In the following, we will focus on two different choices for the function $f$ . Our first choice is to consider an exponentially decaying kernel

f (| α - β |) = e^{- σ | α - β |}

were $σ > 0$ . The second option that we set to explore assumes a power-law decay in the form

f (| α - β |) = \frac{1}{| α - β |^{ρ}}

with $ρ > 0$ . In the next section, we will challenge different scenarios against experimental data.

Comparison with experiments

Setting $ν = 0$ in Equation (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 Figure 4C) and Synechococcus (see Figure 3C). 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 Equation (34).

The results of the comparison between theory and experiments are plotted in Figure 5 and Appendix 1—figure 6. The plot in the main text corresponds to the exponential long-range cell-cell communication (Equation 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 squares) and a decaying function (Equation 45) (blue squares). In all cases, the experimental points are nicely interpolated by the theoretical curves. Accounting for a smooth kernel of interaction yields more accurate fittings.

Appendix 1—figure 6

Download asset Open asset

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, and the squares represent fits to the experimental data with the prediction of the stochastic model using the constant kernel (black squares) and the power-law kernel (blue squares). Lines between symbols are a guide to the eye. 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.

Supplemental methods

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 pCSAL33 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 CSV3 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 CSV3, and Rev plasmid CSV3). 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 WT chromosomes, in BG11 medium supplemented with antibiotics. The Anabaena strains used as recipients in the conjugations were the WT 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_{h e t R - g f p}$ transcriptional fusion in a WT background has been reported previously (Corrales-Guerrero et al., 2013).

Appendix 1—figure 7

Download asset Open asset

Construction of an *Anabaena* $Δ$ *kaiABC* deletion mutant in the PpecB-gfp genetic background.

PCR, DNA restriction/ligation, and transformation into *Escherichia coli* were performed by standard techniques. Conjugation from *E. coli* to *Anabaena* was performed as described by Elhai et al., 1997, and sucrose (sacB)-based positive selection for double recombinants was performed 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.

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 hr 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 oligodeoxynucleotide sequences used for construction are listed in Appendix 1—table 2.

Culture conditions

Strains and derived strains were grown photoautotrophically in BG11 medium containing NaNO₃, 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 of photons [spectrum centered at 450 nm]) from a cool-white LED array. When required, Sm and 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 hr 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 Cm was added.

Sequence

Region alr0523

ttttGAAttcGCTTATAAACAGCAGTTAACAGGCTATTAATAATTATGATGAATATTAGTCTTGAAAATTATTATCATTTTTAGCCACAAAATAAAAACCGAAGATTGTTATCTAAAACACAGATTTTTTTTGATAAAAATCCTGCTTTAGCTTTGTTTATCAATTTCTACTGTCAGCATGAAAGTATATATAAGCTAGATTTGCTACCGCCGTAGAAATTATACTCAAGCATTTATTAAGTAAAGTAAACAATATTATTTAATTGTTAAACGGTTTTTCACCTAAAGGCATTACAACCTTTATAAAACAATAATTTTTAGAGAAGTCTGTATGGATTAACACTATCTTACAAAACTTAATAATAAATTCACTCAAGCCATATATAAAATACGAACTAAGGTTTTCAAAGCAATTAAAAAACAGAAAACAAATTGCTGTTTTCCTAATAAGAAATCGCACATGATTGCAGTCCCCCAAGATAATATATGGGAGTTAATTCTGCTTGAATATTTGTTTTCTAAATAGTAAGAATAATTGCAATCGACCTTATAAAAAGCTGCAATGACCTTTAGGAGGAAAGAAAGATGCTCGATGCTTTTTCCAAAGTAGTTGAACAGGCAggcggtggaggtagca

Search of putative RpaA regulatory sequences of the kaiABC, rpaA, pecB, and ftsZ promoter regions in Anabaena

The positions for all transcription start sites were extracted from Mitschke et al., 2011. For in silico analysis of the Anabaena promoter regions, we use a position-specific probability matrix that defines the RpaA binding motif in Synechococcus (Markson et al., 2013) in FIMO tool (Grant et al., 2011).

Data availability

Source data files, Video 1 and Key resources table have been deposited in Dryad (https://doi.org/10.5061/dryad.sxksn031n).

The following data sets were generated

1. Arbel-Goren R
2. Buonfiglio V
3. Patti FD
4. Camargo S
5. Zhitnitsky A
6. Valladares A
7. Flores E
8. Herrero A
9. Fanelli D
10. Stavans J
(2020) Dryad Digital Repository
Demographic noise seeds robust synchronized oscillations in the circadian clock of Anabaena.

https://doi.org/10.5061/dryad.sxksn031n

The following previously published data sets were used

1. Kaneko T
2. Nakamura Y
3. Wolk PC
4. Kuritz T
5. Sasamoto S
6. Watanabe A
7. Iriguchi M
8. Ishikawa A
9. Kawashima K
10. Kimura T
11. Kishida Y
12. Kohara M
13. Matsumoto M
14. Matsuno A
15. Muraki A
16. Nakazaki N
17. Shimpo S
18. Sugimoto M
19. Takazawa M
20. Yamada M
21. Yasuda M
22. Tabata S
(2001) Accession numbers, AP003581 (nucleotide positions 1-348,050), AP003582 (348,001- 690,650), AP003583 (690,601-1,030,250), AP003584 (1,030,251-1,378,550), AP003585 (1,378,501-1,720,550), AP003586 (1,720,501-2,069,550), AP003587 (2,069,501- 2,413,050), AP003588 (2,413,001-2,747,520), AP003589 (2,747,471-3,089,350), AP003590 (3,089,301-3,422,800), AP003591 (3,422,751-3,770,150), AP003592 (3,770,101- 4,118,350), AP003593 (4,118,301-4,451,850), AP003594 (4,451,801-4,795,050), AP003595 (4,795,001-5,142,550), AP003596 (5,142,501-5,491,050), AP003597 (5,491,001- 5,833,850), AP003598 (5,833,801-6,176,600), and AP003599 (6,176,551-6,413) The sequence data analyzed in this study have been registered in DDBJ/GenBank/EMBL
ID 771. Complete Genomic Sequence of the Filamentous Nitrogen-fixing Cyanobacterium Anabaena sp. Strain PCC 7120.

http://www.kazusa.or.jp/cyanobase/
(2013) NCBI Gene Expression Omnibus
ID GSE50922. Circadian Control of Global Gene Expression by the Cyanobacterial Master Regulator RpaA.

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32317
1. Copeland A
2. Lucas S
3. Lapidus A
4. Barry K
5. Detter JC
6. Glavina T
7. Hammon N
8. Israni S
9. Pitluck S
10. Schmutz J
11. Larimer F
12. Land M
13. Kyrpides N
14. Lykidis A
15. Golden S
16. Richardson P
(2014) NCBI GenBank
ID CP000100.1. Synechococcus elongatus PCC 7942, complete genome.

https://www.ncbi.nlm.nih.gov/nuccore/CP000100

References

(2007) Cyanobacterial clock, a stable phase oscillator with negligible intercellular coupling
PNAS 104:7051–7056.

https://doi.org/10.1073/pnas.0609315104
- PubMed
- Google Scholar
1. Beck C
2. Hertel S
3. Rediger A
4. Lehmann R
5. Wiegard A
6. Kölsch A
7. Heilmann B
8. Georg J
9. Hess WR
10. Axmann IM
(2014) Daily expression pattern of protein-encoding genes and small noncoding RNAs in Synechocystis sp. strain PCC 6803
Applied and Environmental Microbiology 80:5195–5206.

https://doi.org/10.1128/AEM.01086-14
- PubMed
- Google Scholar
(2005) Circadian rhythms from multiple oscillators: lessons from diverse organisms
Nature Reviews Genetics 6:544–556.

https://doi.org/10.1038/nrg1633
- PubMed
- Google Scholar
(2010) A sequestration feedback determines dynamics and temperature entrainment of the KaiABC circadian clock
Molecular Systems Biology 6:389.

https://doi.org/10.1038/msb.2010.44
- PubMed
- Google Scholar
1. Cai YP
2. Wolk CP
(1990) Use of a conditionally lethal gene in Anabaena sp. strain PCC 7120 to select for double recombinants and to entrap insertion sequences
Journal of Bacteriology 172:3138–3145.

https://doi.org/10.1128/JB.172.6.3138-3145.1990
- PubMed
- Google Scholar
(2007) Stochastic gene expression out-of-steady-state in the cyanobacterial circadian clock
Nature 450:1249–1252.

https://doi.org/10.1038/nature06395
- PubMed
- Google Scholar
1. Challenger JD
2. McKane AJ
(2013) Synchronization of stochastic oscillators in biochemical systems
Physical Review E 88:012107.

https://doi.org/10.1103/PhysRevE.88.012107
- PubMed
- Google Scholar
1. Chew J
2. Leypunskiy E
3. Lin J
4. Murugan A
5. Rust MJ
(2018) High protein copy number is required to suppress stochasticity in the cyanobacterial circadian clock
Nature Communications 9:3004.

https://doi.org/10.1038/s41467-018-05109-4
- PubMed
- Google Scholar
1. Cohen SE
2. Golden SS
(2015) Circadian rhythms in cyanobacteria
Microbiology and Molecular Biology Reviews 79:373–385.

https://doi.org/10.1128/MMBR.00036-15
- PubMed
- Google Scholar
(2013) Functional dissection and evidence for intercellular transfer of the heterocyst-differentiation PatS morphogen
Molecular Microbiology 88:1093–1105.

https://doi.org/10.1111/mmi.12244
- PubMed
- Google Scholar
(2015) Spatial fluctuations in expression of the heterocyst differentiation regulatory gene hetR in Anabaena filaments
PLOS Genetics 11:e1005031.

https://doi.org/10.1371/journal.pgen.1005031
- PubMed
- Google Scholar
(2018) FtsZ of filamentous, Heterocyst-Forming cyanobacteria has a conserved N-Terminal peptide required for normal FtsZ polymerization and cell division
Frontiers in Microbiology 9:2260.

https://doi.org/10.3389/fmicb.2018.02260
- PubMed
- Google Scholar
(2010) A synchronized quorum of genetic clocks
Nature 463:326–330.

https://doi.org/10.1038/nature08753
- PubMed
- Google Scholar
(2018) Robust stochastic turing patterns in the development of a one-dimensional cyanobacterial organism
PLOS Biology 16:e2004877.

https://doi.org/10.1371/journal.pbio.2004877
- PubMed
- Google Scholar
1. Dong G
2. Yang Q
3. Wang Q
4. Kim YI
5. Wood TL
6. Osteryoung KW
7. van Oudenaarden A
8. Golden SS
(2010) Elevated ATPase activity of KaiC applies a circadian checkpoint on cell division in Synechococcus elongatus
Cell 140:529–539.

https://doi.org/10.1016/j.cell.2009.12.042
- PubMed
- Google Scholar
(1997) Reduction of conjugal transfer efficiency by three restriction activities of Anabaena sp. strain PCC 7120
Journal of Bacteriology 179:1998–2005.

https://doi.org/10.1128/JB.179.6.1998-2005.1997
- PubMed
- Google Scholar
1. Endo M
(2016) Tissue-specific circadian clocks in plants
Current Opinion in Plant Biology 29:44–49.

https://doi.org/10.1016/j.pbi.2015.11.003
- PubMed
- Google Scholar
1. Enright JT
(1980) Temporal precision in circadian systems: a reliable neuronal clock from unreliable components?
Science 209:1542–1545.

https://doi.org/10.1126/science.7433976
- PubMed
- Google Scholar
1. Fraser AM
2. Swinney HL
(1986) Independent coordinates for strange attractors from mutual information
Physical Review A 33:1134–1140.

https://doi.org/10.1103/PhysRevA.33.1134
- PubMed
- Google Scholar
1. Gan S
2. O'Shea EK
(2017) An unstable singularity underlies stochastic phasing of the circadian clock in individual cyanobacterial cells
Molecular Cell 67:659–672.

https://doi.org/10.1016/j.molcel.2017.07.015
- PubMed
- Google Scholar
1. Garces RG
2. Wu N
3. Gillon W
4. Pai EF
(2004) Anabaena circadian clock proteins KaiA and KaiB reveal a potential common binding site to their partner KaiC
The EMBO Journal 23:1688–1698.

https://doi.org/10.1038/sj.emboj.7600190
- PubMed
- Google Scholar
(2004) Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing
PNAS 101:10955–10960.

https://doi.org/10.1073/pnas.0307095101
- PubMed
- Google Scholar
1. Gillespie DT
(1977) Exact stochastic simulation of coupled chemical reactions
The Journal of Physical Chemistry 81:2340–2361.

https://doi.org/10.1021/j100540a008
- Google Scholar
1. Golden SS
(2020) Principles of rhythmicity emerging from cyanobacteria
European Journal of Neuroscience 51:13–18.

https://doi.org/10.1111/ejn.14434
- PubMed
- Google Scholar
1. Golden SS
2. Sherman LA
(1984) Optimal conditions for genetic transformation of the Cyanobacterium anacystis nidulans R2
Journal of Bacteriology 158:36–42.

https://doi.org/10.1128/JB.158.1.36-42.1984
- PubMed
- Google Scholar
1. Gould PD
2. Domijan M
3. Greenwood M
4. Tokuda IT
5. Rees H
6. Kozma-Bognar L
7. Hall AJ
8. Locke JC
(2018) Coordination of robust single cell rhythms in the Arabidopsis circadian clock via spatial waves of gene expression
eLife 7:e31700.

https://doi.org/10.7554/eLife.31700
- PubMed
- Google Scholar
(2011) FIMO: scanning for occurrences of a given motif
Bioinformatics 27:1017–1018.

https://doi.org/10.1093/bioinformatics/btr064
- PubMed
- Google Scholar
(2012) Variation of the folding and dynamics of the Escherichia coli chromosome with growth conditions
Molecular Microbiology 86:1318–1333.

https://doi.org/10.1111/mmi.12071
- PubMed
- Google Scholar
(2010) Systems biology of cellular rhythms: from cacophony to symphony
Current Opinion in Genetics & Development 20:571–573.

https://doi.org/10.1016/j.gde.2010.10.003
- PubMed
- Google Scholar
(2019) On chaotic dynamics in transcription factors and the associated effects in differential gene regulation
Nature Communications 10:71.

https://doi.org/10.1038/s41467-018-07932-1
- PubMed
- Google Scholar
(2016) The multicellular nature of filamentous heterocyst-forming cyanobacteria
FEMS Microbiology Reviews 40:831–854.

https://doi.org/10.1093/femsre/fuw029
- PubMed
- Google Scholar
1. Iijima H
2. Shirai T
3. Okamoto M
4. Kondo A
5. Hirai MY
6. Osanai T
(2015) Changes in primary metabolism under light and dark conditions in response to overproduction of a response regulator RpaA in the unicellular Cyanobacterium synechocystis sp. PCC 6803
Frontiers in Microbiology 6:888.

https://doi.org/10.3389/fmicb.2015.00888
- PubMed
- Google Scholar
1. Ishiura M
2. Kutsuna S
3. Aoki S
4. Iwasaki H
5. Andersson CR
6. Tanabe A
7. Golden SS
8. Johnson CH
9. Kondo T
(1998) Expression of a gene cluster kaiABC as a circadian feedback process in cyanobacteria
Science 281:1519–1523.

https://doi.org/10.1126/science.281.5382.1519
- PubMed
- Google Scholar
Conference
(2018) A novel method for topological embedding of time-series data
26th European Signal Processing Conference (EUSIPCO). IEEE. pp. 2350–2354.

https://doi.org/10.23919/EUSIPCO.2018.8553502
- Google Scholar
1. Kizawa A
2. Osanai T
(2020) Overexpression of the response regulator rpaA causes an impaired cell division in the Cyanobacterium synechocystis sp. PCC 6803
The Journal of General and Applied Microbiology 66:121–128.

https://doi.org/10.2323/jgam.2020.01.004
- PubMed
- Google Scholar
1. Kucho K
2. Okamoto K
3. Tsuchiya Y
4. Nomura S
5. Nango M
6. Kanehisa M
7. Ishiura M
(2005) Global analysis of circadian expression in the Cyanobacterium synechocystis sp. strain PCC 6803
Journal of Bacteriology 187:2190–2199.

https://doi.org/10.1128/JB.187.6.2190-2199.2005
- PubMed
- Google Scholar
1. Kushige H
2. Kugenuma H
3. Matsuoka M
4. Ehira S
5. Ohmori M
6. Iwasaki H
(2013) Genome-wide and heterocyst-specific circadian gene expression in the filamentous Cyanobacterium anabaena sp. strain PCC 7120
Journal of Bacteriology 195:1276–1284.

https://doi.org/10.1128/JB.02067-12
- PubMed
- Google Scholar
1. Lambert G
2. Chew J
3. Rust MJ
(2016) Costs of Clock-Environment misalignment in individual cyanobacterial cells
Biophysical Journal 111:883–891.

https://doi.org/10.1016/j.bpj.2016.07.008
- PubMed
- Google Scholar
(2013) Circadian control of global gene expression by the cyanobacterial master regulator RpaA
Cell 155:1396–1408.

https://doi.org/10.1016/j.cell.2013.11.005
- PubMed
- Google Scholar
Book
1. Marple SL
(1987)
Digital Spectral Analysis with Applications

Englewood Cliffs, NJ: Prentice Hall.
- Google Scholar
1. Masuda K
2. Kitaoka R
3. Ukai K
4. Tokuda IT
5. Fukuda H
(2017) Multicellularity enriches the entrainment of Arabidopsis circadian clock
Science Advances 3:e1700808.

https://doi.org/10.1126/sciadv.1700808
- PubMed
- Google Scholar
1. McKane AJ
2. Newman TJ
(2005) Predator-prey cycles from resonant amplification of demographic stochasticity
Physical Review Letters 94:218102.

https://doi.org/10.1103/PhysRevLett.94.218102
- PubMed
- Google Scholar
(2004) Resilient circadian oscillator revealed in individual cyanobacteria
Nature 430:81–85.

https://doi.org/10.1038/nature02533
- PubMed
- Google Scholar
1. Mitschke J
2. Vioque A
3. Haas F
4. Hess WR
5. Muro-Pastor AM
(2011) Dynamics of transcriptional start site selection during nitrogen stress-induced cell differentiation in Anabaena sp. PCC7120
PNAS 108:20130–20135.

https://doi.org/10.1073/pnas.1112724108
- PubMed
- Google Scholar
1. Mohamed A
2. Jansson C
(1989) Influence of light on accumulation of photosynthesis-specific transcripts in the Cyanobacterium synechocystis 6803
Plant Molecular Biology 13:693–700.

https://doi.org/10.1007/BF00016024
- PubMed
- Google Scholar
(1996) Circadian gating of cell division in cyanobacteria growing with average doubling times of less than 24 hours
PNAS 93:10183–10188.

https://doi.org/10.1073/pnas.93.19.10183
- PubMed
- Google Scholar
Book
1. Mori T
(2009) Cell division cycles and circadian rhythms
In: Ditty J. L, Mackey S. R, Johnson C. H, editors. Bacterial Circadian Programs. Springer. pp. 183–204.

https://doi.org/10.1016/j.ceb.2013.07.013
- Google Scholar
1. Mori T
2. Johnson CH
(2001) Independence of circadian timing from cell division in cyanobacteria
Journal of Bacteriology 183:2439–2444.

https://doi.org/10.1128/JB.183.8.2439-2444.2001
- PubMed
- Google Scholar
1. Mullineaux CW
2. Mariscal V
3. Nenninger A
4. Khanum H
5. Herrero A
6. Flores E
7. Adams DG
(2008) Mechanism of intercellular molecular exchange in heterocyst-forming cyanobacteria
The EMBO Journal 27:1299–1308.

https://doi.org/10.1038/emboj.2008.66
- PubMed
- Google Scholar
1. Nakajima M
2. Imai K
3. Ito H
4. Nishiwaki T
5. Murayama Y
6. Iwasaki H
7. Oyama T
8. Kondo T
(2005) Reconstitution of circadian oscillation of cyanobacterial KaiC phosphorylation in vitro
Science 308:414–415.

https://doi.org/10.1126/science.1108451
- PubMed
- Google Scholar
(2015) Intercellular diffusion of a fluorescent sucrose analog via the septal junctions in a filamentous Cyanobacterium
mBio 6:e02109.

https://doi.org/10.1128/mBio.02109-14
- PubMed
- Google Scholar
(2017) A roadmap for the computation of persistent homology
EPJ Data Science 6:17.

https://doi.org/10.1140/epjds/s13688-017-0109-5
- PubMed
- Google Scholar
(2016) Discrete gene replication events drive coupling between the cell cycle and circadian clocks
PNAS 113:4063–4068.

https://doi.org/10.1073/pnas.1507291113
- PubMed
- Google Scholar
1. Pereira CMM
2. de Mello RF
(2015) Persistent homology for time series and spatial data clustering
Expert Systems with Applications 42:6026–6038.

https://doi.org/10.1016/j.eswa.2015.04.010
- Google Scholar
1. Pfaffl MW
(2001) A new mathematical model for relative quantification in real-time RT-PCR
Nucleic Acids Research 29:e45.

https://doi.org/10.1093/nar/29.9.e45
- PubMed
- Google Scholar
1. Pittayakanchit W
2. Lu Z
3. Chew J
4. Rust MJ
5. Murugan A
(2018) Biophysical clocks face a trade-off between internal and external noise resistance
eLife 7:e37624.

https://doi.org/10.7554/eLife.37624
- PubMed
- Google Scholar
1. Reppert SM
2. Weaver DR
(2002) Coordination of circadian timing in mammals
Nature 418:935–941.

https://doi.org/10.1038/nature00965
- PubMed
- Google Scholar
1. Rust MJ
2. Markson JS
3. Lane WS
4. Fisher DS
5. O'Shea EK
(2007) Ordered phosphorylation governs oscillation of a three-protein circadian clock
Science 318:809–812.

https://doi.org/10.1126/science.1148596
- PubMed
- Google Scholar
1. Schmelling NM
2. Lehmann R
3. Chaudhury P
4. Beck C
5. Albers SV
6. Axmann IM
7. Wiegard A
(2017) Minimal tool set for a prokaryotic circadian clock
BMC Evolutionary Biology 17:169.

https://doi.org/10.1186/s12862-017-0999-7
- PubMed
- Google Scholar
(1992) Genes encoding the phycobilisome rod substructure are clustered on the Anabaena chromosome: characterization of the phycoerythrocyanin operon
Journal of Bacteriology 174:2640–2647.

https://doi.org/10.1128/JB.174.8.2640-2647.1992
- PubMed
- Google Scholar
Book
(2014) Javaplex: A research software package for persistent (co) homology
In: Hong H, Yap C, editors. International Congress on Mathematical Software. Springer. pp. 129–136.

https://doi.org/10.1007/978-3-662-44199-2_23
- Google Scholar
1. Teng SW
2. Mukherji S
3. Moffitt JR
4. de Buyl S
5. O'Shea EK
(2013) Robust circadian oscillations in growing cyanobacteria require transcriptional feedback
Science 340:737–740.

https://doi.org/10.1126/science.1230996
- PubMed
- Google Scholar
1. Tseng R
2. Goularte NF
3. Chavan A
4. Luu J
5. Cohen SE
6. Chang YG
7. Heisler J
8. Li S
9. Michael AK
10. Tripathi S
11. Golden SS
12. LiWang A
13. Partch CL
(2017) Structural basis of the day-night transition in a bacterial circadian clock
Science 355:1174–1180.

https://doi.org/10.1126/science.aag2516
- PubMed
- Google Scholar
1. Tsimring LS
(2014) Noise in biology
Reports on Progress in Physics 77:026601.

https://doi.org/10.1088/0034-4885/77/2/026601
- PubMed
- Google Scholar
1. Uzumaki T
2. Fujita M
3. Nakatsu T
4. Hayashi F
5. Shibata H
6. Itoh N
7. Kato H
8. Ishiura M
(2004) Crystal structure of the C-terminal clock-oscillator domain of the cyanobacterial KaiA protein
Nature Structural & Molecular Biology 11:623–631.

https://doi.org/10.1038/nsmb781
- PubMed
- Google Scholar
Book
1. van Kampen NG
(2014) Stochastic Processes in Physics and Chemistry
Springer.

https://doi.org/10.1007/3-540-45396-2
- Google Scholar
(2010) Circadian gating of the cell cycle revealed in single cyanobacterial cells
Science 327:1522–1526.

https://doi.org/10.1126/science.1181759
- PubMed
- Google Scholar
1. Zankoc C
2. Fanelli D
3. Ginelli F
4. Livi R
(2017) Intertangled stochastic motifs in networks of excitatory-inhibitory units
Physical Review E 96:022308.

https://doi.org/10.1103/PhysRevE.96.022308
- PubMed
- Google Scholar
1. Zankoc C
2. Fanelli D
3. Ginelli F
4. Livi R
(2019) Desynchronization and pattern formation in a noisy feed-forward oscillator network
Physical Review E 99:012303.

https://doi.org/10.1103/PhysRevE.99.012303
- PubMed
- Google Scholar
1. Zavřel T
2. Faizi M
3. Loureiro C
4. Poschmann G
5. Stühler K
6. Sinetova M
7. Zorina A
8. Steuer R
9. Červený J
(2019) Quantitative insights into the cyanobacterial cell economy
eLife 8:e42508.

https://doi.org/10.7554/eLife.42508
- PubMed
- Google Scholar
1. Zhao W
2. Guo Q
3. Zhao J
(2007) A membrane-associated Mn-superoxide dismutase protects the photosynthetic apparatus and nitrogenase from oxidative damage in the Cyanobacterium anabaena sp. PCC 7120
Plant and Cell Physiology 48:563–572.

https://doi.org/10.1093/pcp/pcm025
- PubMed
- Google Scholar

Article and author information

Author details

Rinat Arbel-Goren

Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel

Contribution
Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-7253-2036
Valentina Buonfiglio

Dipartimento di Fisica e Astronomia, Università di Firenze, INFN and CSDC, Sesto Fiorentino, Italy

Contribution
Software, Formal analysis

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1528-2955
Francesca Di Patti

Consiglio Nazionale delle Ricerche, Istituto dei Sistemi Complessi, Sesto Fiorentino, Italy

Contribution
Conceptualization, Resources, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1288-0079
Sergio Camargo

Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel

Contribution
Resources, Formal analysis, Investigation, Visualization

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5688-3215
Anna Zhitnitsky

Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel

Contribution
Formal analysis

"This ORCID iD identifies the author of this article:" 0000-0001-9598-4812
Ana Valladares

Instituto de Bioquímica Vegetal y Fotosíntesis, CSIC and Universidad de Sevilla, Sevilla, Spain

Contribution
Resources, Validation, Investigation

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6683-1355
Enrique Flores

Instituto de Bioquímica Vegetal y Fotosíntesis, CSIC and Universidad de Sevilla, Sevilla, Spain

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7605-7343
Antonia Herrero

Instituto de Bioquímica Vegetal y Fotosíntesis, CSIC and Universidad de Sevilla, Sevilla, Spain

Contribution
Conceptualization, Resources, Supervision, Funding acquisition, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1071-6590
Duccio Fanelli

Dipartimento di Fisica e Astronomia, Università di Firenze, INFN and CSDC, Sesto Fiorentino, Italy

Contribution
Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Visualization, Methodology, Project administration, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8545-9424
Joel Stavans

Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel

Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing

For correspondence
joel.stavans@weizmann.ac.il

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0396-7797

Funding

Minerva Foundation

Joel Stavans

Fondazione Ente Cassa di Risparmio di Firenze

Duccio Fanelli

European Regional Development Fund (BUF2016-77097-P)

Antonia Herrero

European Regional Development Fund (BFU2017-88202-P)

Enrique Flores

Ministry of Foreign Affairs (EXPLICS)

Francesca Di Patti

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank M Feingold, M Camarena, T Unger, JE Frías, N Rosenthal, and O Vardi for help, and M Rust for strains. This study was supported by the Minerva Foundation (JS); grant numbers BFU2017-88202-P (EF) and BUF2016-77097-P (AH) from Plan Nacional de Investigación, Spain, co-financed by the European Regional Development Fund; Fondazione Ente Cassa di Risparmio, Florence, Italy (DF, FDP, VB); Project EXPLICS granted by the Italian Ministry of Foreign Affairs and International Cooperation (FDP).

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.