Robust, coherent, and synchronized circadian clock-controlled oscillations along Anabaena filaments
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 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 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).
To characterize quantitatively the degree of synchronization between cellular clocks along a filament, we used the synchronization index (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 , 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 () (Table 1). The value of 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.
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 in different filaments oscillated nearly in phase mainly during the first oscillations. To evaluate quantitatively the degree of phase synchronization between filaments, we calculated by choosing one cell per filament in a number of filaments measured simultaneously (Materials and methods). We obtained , 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 -specific, but clock-independent, as it was also observed in filaments of the 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 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 as a function of distance along a filament (see Materials and methods). We found that 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 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 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 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 expression is due primarily to significant coupling between the clocks in neighboring cells due to cell-cell communication. Of note, the value of calculated for a background, in which 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 . To quantify these variations, we calculated the square of the coefficient of variation , where denotes the standard deviation of the fluorescence intensity of contiguous cells along a filament. A plot of 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 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 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 for the different pairs and obtained , , and . 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 evaluated for the three genes was 0.78 ± 0.08.
To expose an underlying oscillatory behavior in the kai genes data and support the notion that the oscillations are circadian, we applied persistent homology methods (Pereira and de Mello, 2015; Otter et al., 2017) to two-dimensional phase portraits of their respective time series (see Appendix 1 – supplemental 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 hr, hr, and () 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 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 for 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 (), T (), and D (). 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 µM matches the concentration reported previously ( µM; see Rust et al., 2007).
Theoretical model of arrays of coupled noisy clocks
Next, we generalized the single-clock model above to an array of coupled circadian clocks as in Anabaena, which is endowed with cell-cell communication via septal proteins (Herrero et al., 2016; 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.
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 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 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 indicated that the cell-cell variability along filaments as measured by 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 in 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 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
Strains
Strains bearing a chromosomally encoded , 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 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 protocolStrains and derived strains were grown photoautotrophically in BG11 medium containing NaNO3, 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 protocolStrains 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 NaNO3 and 10 mM NaHCO3, 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 protocolFor 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 (Pfaffl, 2001).
Image segmentation
Request a detailed protocolAll 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 protocolSynchronization was measured by the order parameter proposed by Garcia-Ojalvo et al., 2004:
where denotes a time average, indicates an average over all cells, and µ denotes the average of the fluorescence intensity of each cell fi. Hence, is defined as the ratio of the standard deviation of to the standard deviation of fi, 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. was then calculated, and this procedure was repeated for different choices of cells, at least three times for each experiment. All the evaluations of were carried out over a full period of oscillation in either one of the first two oscillations, except for the background, for which 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 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 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
where
where fi denotes the fluorescence intensity in cell in a filament of cells.
Robustness of Anabaena dataset fit to the Synechococcus parameters
Request a detailed protocolEach kinetic parameter () was randomly selected from a Gaussian distribution centered at the nominal value () estimated by Rust et al., 2007 using the normrnd MATLAB function. The standard deviation is assigned as . 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 and [KaiA] (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 (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 whose elements 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 .
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:
where the superscript indicates that the species belongs to the cell . The interconversion rates between phosphoforms are given by
with
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 stands for the total number of S molecules in cell at time , and, in the same way, we will use to denote the total number of molecules for . The total number of KaiC phosphoforms remains constant in each cell of the filament, and we set this value equal to . 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).
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 , as shown in Appendix 1—figure 1, counts the number of neighbors of any given cell. The chemical equations corresponding to the envisaged coupling read
where and label connected nodes.
Since the total number of molecules within each cell of the filament remains constant, the variable can be neglected because it can be recovered from the conservation law . For this reason, at each time the state of the system is specified by the -dimensional vector
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
while the interconversion transitions are
To account for cell-cell communication, we also add the transitions for long-range interaction between two connected cells and :
The dynamics is fully described in terms of the following master equation:
whose explicit version reads
The master equation can be written in a more compact form though the use of the step operators that act on a generic function as and :
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 into two contributions, the deterministic one, which is obtained taking the limit denoted by , and the fluctuations whose amplitude scales as , as
We introduce a new probability function , where the is a -dimensional vector whose components are
so that the left-hand side of the master Equation (15) can be recast in the form
with .
Operating an expansion of the step operators with respect to the small parameter , performing the change of variable (Equation 16) into the right-hand side of the master equation, and collecting terms proportional to , we get
where
and the matrix denotes the discrete Laplacian operator. The concentration reads and 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
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 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
where
and calculate its three eigenvalues. The nature of the equilibrium points and their stability is summarized in Appendix 1—figure 2, where we plot the equilibrium points to varying . Continuous lines denote stable equilibrium points, while dashed lines 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 . The value 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 , the system undergoes a saddle-node bifurcation with the sudden creation of two additional fixed points, one stable and the other unstable (region III).
To better investigate the region of the 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 and . 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 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.
Analysis of fluctuations
Considering the next-to-leading order of the van Kampen expansion and collecting together terms proportional to , we find the differential equation for the probability distribution , which results in the following Fokker–Planck equation:
where is
The two matrices and have dimension , and they can be expressed as the sum of two distinct terms, one due to the spatial contribution that we label with the superscript , and the other related to the interactions within each cell that we label as . Moreover, we assume that all the entries of the two matrices are evaluated at the equilibrium point . Thus, we have
with
where is the Jacobian matrix relative to system (Equation 19) evaluated at the steady state. For , we have
with
Spatial correlations
The Fokker–Planck Equation (25) describes the fluctuations about the equilibrium point , and it is equivalent to the following coupled Langevin equations:
where is the component of vector defined in Equation (17) while is the ith component of the noise with the following properties:
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) is defined by
where is the temporal Fourier transform of the th component of . The degree of correlation between species in cell and species in cell can therefore be quantified through the complex coherence function (CCF), a matrix , whose entries are the ratio of the elements of PSDM to a normalization coefficient:
To obtain an analytic expression of , we solve the Langevin Equation (31) performing a temporal Fourier transform. In this way, for each component of the fluctuations vector on cell we get
with , where denotes the identity matrix. Substituting Equation (36) into Equation (34), one eventually obtains
where is the complex conjugate transpose of and use has been made of the following properties:
Accounting for a smooth kernel of interactions
In this section, we will extend the analysis so as to account for a 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, . 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:
where is a constant and 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 as
and the corresponding Laplacian matrix now becomes
Accordingly, the two matrices (Equations 27 and 29) appearing in the Fokker–Planck Equation (25) must be replaced by
In the following, we will focus on two different choices for the function . Our first choice is to consider an exponentially decaying kernel
were . The second option that we set to explore assumes a power-law decay in the form
with . In the next section, we will challenge different scenarios against experimental data.
Comparison with experiments
Setting in Equation (37) and considering the diagonal terms and returns 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 when postulating a uniform range of interaction or by varying and the characteristic parameter of the decay function , 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 and . Notice that 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 (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.
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 transcriptional fusion in a WT background has been reported previously (Corrales-Guerrero et al., 2013).
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 NaNO3, 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).
-
Dryad Digital RepositoryDemographic noise seeds robust synchronized oscillations in the circadian clock of Anabaena.https://doi.org/10.5061/dryad.sxksn031n
-
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/EMBLID 771. Complete Genomic Sequence of the Filamentous Nitrogen-fixing Cyanobacterium Anabaena sp. Strain PCC 7120.
-
NCBI Gene Expression OmnibusID GSE50922. Circadian Control of Global Gene Expression by the Cyanobacterial Master Regulator RpaA.
-
NCBI GenBankID CP000100.1. Synechococcus elongatus PCC 7942, complete genome.
References
-
Daily expression pattern of protein-encoding genes and small noncoding RNAs in Synechocystis sp. strain PCC 6803Applied and Environmental Microbiology 80:5195–5206.https://doi.org/10.1128/AEM.01086-14
-
Circadian rhythms from multiple oscillators: lessons from diverse organismsNature Reviews Genetics 6:544–556.https://doi.org/10.1038/nrg1633
-
Synchronization of stochastic oscillators in biochemical systemsPhysical Review E 88:012107.https://doi.org/10.1103/PhysRevE.88.012107
-
Circadian rhythms in cyanobacteriaMicrobiology and Molecular Biology Reviews 79:373–385.https://doi.org/10.1128/MMBR.00036-15
-
Tissue-specific circadian clocks in plantsCurrent Opinion in Plant Biology 29:44–49.https://doi.org/10.1016/j.pbi.2015.11.003
-
Independent coordinates for strange attractors from mutual informationPhysical Review A 33:1134–1140.https://doi.org/10.1103/PhysRevA.33.1134
-
Exact stochastic simulation of coupled chemical reactionsThe Journal of Physical Chemistry 81:2340–2361.https://doi.org/10.1021/j100540a008
-
Principles of rhythmicity emerging from cyanobacteriaEuropean Journal of Neuroscience 51:13–18.https://doi.org/10.1111/ejn.14434
-
Optimal conditions for genetic transformation of the Cyanobacterium anacystis nidulans R2Journal of Bacteriology 158:36–42.https://doi.org/10.1128/JB.158.1.36-42.1984
-
FIMO: scanning for occurrences of a given motifBioinformatics 27:1017–1018.https://doi.org/10.1093/bioinformatics/btr064
-
Variation of the folding and dynamics of the Escherichia coli chromosome with growth conditionsMolecular Microbiology 86:1318–1333.https://doi.org/10.1111/mmi.12071
-
Systems biology of cellular rhythms: from cacophony to symphonyCurrent Opinion in Genetics & Development 20:571–573.https://doi.org/10.1016/j.gde.2010.10.003
-
The multicellular nature of filamentous heterocyst-forming cyanobacteriaFEMS Microbiology Reviews 40:831–854.https://doi.org/10.1093/femsre/fuw029
-
ConferenceA novel method for topological embedding of time-series data26th European Signal Processing Conference (EUSIPCO). IEEE. pp. 2350–2354.https://doi.org/10.23919/EUSIPCO.2018.8553502
-
Overexpression of the response regulator rpaA causes an impaired cell division in the Cyanobacterium synechocystis sp. PCC 6803The Journal of General and Applied Microbiology 66:121–128.https://doi.org/10.2323/jgam.2020.01.004
-
Global analysis of circadian expression in the Cyanobacterium synechocystis sp. strain PCC 6803Journal of Bacteriology 187:2190–2199.https://doi.org/10.1128/JB.187.6.2190-2199.2005
-
Costs of Clock-Environment misalignment in individual cyanobacterial cellsBiophysical Journal 111:883–891.https://doi.org/10.1016/j.bpj.2016.07.008
-
BookDigital Spectral Analysis with ApplicationsEnglewood Cliffs, NJ: Prentice Hall.
-
Predator-prey cycles from resonant amplification of demographic stochasticityPhysical Review Letters 94:218102.https://doi.org/10.1103/PhysRevLett.94.218102
-
BookCell division cycles and circadian rhythmsIn: 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
-
Independence of circadian timing from cell division in cyanobacteriaJournal of Bacteriology 183:2439–2444.https://doi.org/10.1128/JB.183.8.2439-2444.2001
-
Persistent homology for time series and spatial data clusteringExpert Systems with Applications 42:6026–6038.https://doi.org/10.1016/j.eswa.2015.04.010
-
A new mathematical model for relative quantification in real-time RT-PCRNucleic Acids Research 29:e45.https://doi.org/10.1093/nar/29.9.e45
-
Minimal tool set for a prokaryotic circadian clockBMC Evolutionary Biology 17:169.https://doi.org/10.1186/s12862-017-0999-7
-
BookJavaplex: A research software package for persistent (co) homologyIn: 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
-
Noise in biologyReports on Progress in Physics 77:026601.https://doi.org/10.1088/0034-4885/77/2/026601
-
Crystal structure of the C-terminal clock-oscillator domain of the cyanobacterial KaiA proteinNature Structural & Molecular Biology 11:623–631.https://doi.org/10.1038/nsmb781
Article and author information
Author details
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
© 2021, Arbel-Goren et al.
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.
Metrics
-
- 2,837
- views
-
- 347
- downloads
-
- 20
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Physics of Living Systems
Many proteins have been recently shown to undergo a process of phase separation that leads to the formation of biomolecular condensates. Intriguingly, it has been observed that some of these proteins form dense droplets of sizeable dimensions already below the critical concentration, which is the concentration at which phase separation occurs. To understand this phenomenon, which is not readily compatible with classical nucleation theory, we investigated the properties of the droplet size distributions as a function of protein concentration. We found that these distributions can be described by a scale-invariant log-normal function with an average that increases progressively as the concentration approaches the critical concentration from below. The results of this scaling analysis suggest the existence of a universal behaviour independent of the sequences and structures of the proteins undergoing phase separation. While we refrain from proposing a theoretical model here, we suggest that any model of protein phase separation should predict the scaling exponents that we reported here from the fitting of experimental measurements of droplet size distributions. Furthermore, based on these observations, we show that it is possible to use the scale invariance to estimate the critical concentration for protein phase separation.
-
- Computational and Systems Biology
- Physics of Living Systems
Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.