Developmental function and state transitions of a gene expression oscillator in C. elegans

Gene expression oscillators can structure biological events temporally and spatially. Different biological functions benefit from distinct oscillator properties. Thus, finite developmental processes rely on oscillators that start and stop at specific times; a poorly understood behavior. 15 Here, we have characterized a massive gene expression oscillator comprising >3,700 genes in C. elegans larvae. We report that oscillations initiate in embryos, arrest transiently after hatching and in response to perturbation, and cease in adults. Experimental observation of the transitions between oscillatory and non-oscillatory states at a resolution where we can identify bifurcation points reveals an oscillator operating near a Saddle Node on Invariant Cycle (SNIC) bifurcation. 20 These findings constrain the architecture and mathematical models that can represent this oscillator. They also reveal that oscillator arrests occur reproducibly in a specific phase. Since we find oscillations to be coupled to developmental processes, including molting, this characteristic of SNIC bifurcations thus endows the oscillator with the potential to halt larval development at defined intervals, and thereby execute a developmental checkpoint function. 25


Introduction
Gene expression oscillations occur in many biological systems as exemplified by circadian rhythms in metabolism and behavior (Panda et al., 2002), vertebrate somitogenesis (Oates et al., 2012), plant lateral root branching (Moreno-Risueno et al., 2010), and C. elegans larval development (Hendriks et al., 2014). They are well-suited for timekeeping, acting as molecular 5 clocks that can provide a temporal, and thereby also spatial, structure for biological events (Uriu, 2016). This structure may represent external time, as illustrated by circadian clocks, or provide temporal organization of internal processes without direct reference to external time, as illustrated by somitogenesis clocks (Rensing et al., 2001).
Depending on these distinct functions, oscillators require different properties. Thus, 10 robust representation of external time requires a stable period, i.e., the oscillator has to be compensated for variations in temperature and other environmental factors. It also benefits from a phase-resetting mechanism to permit moderate realignments, if needed, to external time.
Intuitively, either feature seems unlikely to benefit developmental oscillators. By contrast, because developmental processes are finite, e.g., an organism has a characteristic number of 15 somites, developmental oscillators need a start and an end. How such changes in oscillator activity occur in vivo, and which oscillator features enable them, is largely unknown (Riedel-Kruse et al., 2007;Shih et al., 2015).
Here, we characterize the recently discovered 'C. elegans oscillator' (Hendriks et al., 2014;Kim et al., 2013) at high temporal resolution and across the entire period of C. elegans 20 development, from embryo to adult. The system is marked by a massive scale where ~3,700 genes exhibit transcript level oscillations that are detectable, with large, stable amplitudes and widely dispersed expression peak times (i.e., peak phases), in lysates of whole animals. For the purpose of this study, and because insufficient information exists on the identities of core oscillator versus output genes, we define the entire system of oscillating genes as 'the oscillator'.
We demonstrate that the oscillations are coupled to molting, i.e., the cyclical process of new cuticle synthesis and old cuticle shedding that occurs at the end of each larval stage. We observe and characterize onset and offset of oscillations both during continuous development and upon 5 perturbation, and find that these changes in the state of the oscillator system (or bifurcations), occur with a sudden change in amplitude. They also occur in a characteristic oscillator phase and thus at specific, recurring intervals. Because of the phase-locking of the oscillator and molting, arrests thus always occur at the same time during larval stages, around molt exit. This time coincides with the recurring windows of activity of a checkpoint that can halt larval development 10 in response to nutritionally poor conditions. Hence, our results indicate that the C. elegans oscillator functions as a developmental clock whose architecture supports a developmental checkpoint function. Indeed, the features of the bifurcations constrain oscillator architecture, excluding a simple negative-loop design, and possible parameters of mathematical models.

Thousands of genes with oscillatory expression during the four larval stages
Although previous reports agreed on the wide-spread occurrence of oscillatory gene expression in C. elegans larvae (Grün et al., 2014;Hendriks et al., 2014;Kim et al., 2013), the published data sets were either insufficiently temporally resolved or too short to characterize oscillations 5 across C. elegans larval development. Hence, to understand the extent and features of these oscillations better, including their continuity throughout development, we performed two extended time course experiments to cover the entire period of post-embryonic development plus early adulthood at hourly resolution. We extracted total RNA from populations of animals synchronized by hatching in the absence of food. light-gray off-diagonals), presumably reflecting progression through the four larval stages.
The larger dataset enabled us to improve on the previous identification of genes with oscillatory expression (Hendriks et al., 2014). Using cosine wave fitting, and an amplitude cutoff of 20.5, we classified 3,739 genes (24 % of total expressed genes) as 'oscillating' (i.e., rhythmically expressed) from TC2 (Fig. 1B, S1C and Table S1; Methods). Relative to the 20 previous result of 2,718 oscillating genes (18.9% of total expressed genes) in mRNA expression data of L3 and L4 animals (Hendriks et al., 2014), this adds 1,240 new genes and excludes 219 of the previously annotated oscillating genes. We consider this latter group to be most likely false positives from the earlier analysis, resulting from the fact that some genes behave substantially different during L4 compared to the preceding stages as shown below.
Visual inspection of a gene expression heatmap of the fused time course (TC3; Fig. 1C) revealed four cycles of gene expression for the oscillating genes. Oscillations were absent during the first few hours of larval development as well as in adulthood, from ~37 hours on, and both 5 their onset and offset appeared to occur abruptly. We will analyze these and additional features of the system and their implications in more detail in the following sections.

'Oscillating' genes are expressed in several tissues with dispersed peak phases
An examination of the calculated peak phases confirmed the visual impression that individual 10 transcripts peaked at a wide variety of time points, irrespective of expression amplitude (Fig.   1D). In circadian rhythms, peak phase distributions are typically clustered into three or fewer groups when examined in a specific tissue (Koike et al., 2012;Korenčič et al., 2014). However, the identity of oscillating genes differs across cell types and tissues, and for those genes that oscillate in multiple tissues, phases can differ among tissues (Zhang et al., 2014). Hence, we 15 wondered whether the broad peak phase distribution was a consequence of our analysis of RNA from whole animals, whereas individual tissues might exhibit a more defined phase distribution.
To understand in which tissues oscillations occur, we utilized a previous annotation of tissue-specifically expressed genes (Cao et al., 2017). 1,298, and thus a substantial minority (~35%) of oscillating genes, fell in this category for seven different tissues. They were strongly 20 (~2.5-fold) enriched in the hypodermis (epidermis) and pharynx, and more modestly (≤1.5-fold) in glia and intestine (Fig. 1E). By contrast, oscillating genes were greatly depleted from body wall muscle, neurons, and gonad. Hence, oscillatory gene expression occurs indeed in multiple tissues. However, although peak phase distributions deviated for each tissue to some degree from that seen for all oscillating genes, they were still widely distributed for each individual tissue ( Fig. 1F).
We conclude that a wide dispersion of peak phases appears to be an inherent oscillator feature rather than the result of a convoluted output of multiple, tissue-specific oscillators with 5 distinct phase preferences.

Oscillations initiate with a time lag in L1
The observation that oscillations were undetectable during the first few hours of larval development and started only after > 5 hours into L1 (Fig. 1A, C) surprised us. Hence, we 10 performed a separate experiment that covered the first 24 hours of larval development (TC4).
This confirmed our initial finding of a lack of oscillations during the first few hours of larval development ( Fig. 2A, B).
To understand how oscillations initiate after the initial quiescence, we looked at individual genes and observed that the start of detectable oscillations differed for individual 15 genes ( Fig. 2A, C). Nonetheless, the occurrence of first peaks was globally well correlated to the peak phases calculated from data in Fig. 1 (Fig. 2D, E). Moreover, the transcript levels of many genes with a late-occurring (11 -13 hours) first peak proceeded through a trough before reaching their first peak as exemplified in Fig. 2C for F11E6.3. We conclude that initiation of larval development after hatching is accompanied by a time lag prior to transition into an 20 oscillatory state, and that oscillations exhibit a structure of phase-locked gene expressing patterns as soon as they become detectable.

L1 larvae undergo an extended intermolt
Although the gene expression oscillations occur in the context of larval development, functional connections have been lacking. However, genes encoding cuticular components were reported to be enriched among previously identified oscillating genes (Hendriks et al., 2014;Kim et al., 2013), and Gene Ontology (GO-) term analysis of the new extended set of oscillating genes 5 confirms that the top 12 enriched terms all linked to cuticle formation and molting, or protease activity (Fig. 3A). These findings, and the fact that molting is itself a rhythmic process, repeated at the end of each larval stage, suggest the possibility of a functional link between molting and gene expression oscillations.
If such a link were true, we would predict that the initial period of quiescence in the early 10 L1 stage be accompanied by a lengthened stage, and, specifically, intermolt duration. Indeed, using a luciferase-based assay that reveals the period of behavioral quiescence, or lethargus, that is associated with the molt (Fig. S2A-C), others had previously reported an extended L1 relative to other larval stages (Olmedo et al., 2015). However, they reported an extension of both molt and intermolt. 15 As the previously used luciferase-expressing transgenic strains developed relatively slowly and with limited synchrony across animals, presumably due to their specific genetic make-up, we repeated the experiment with a newly generated strain that expressed luciferase from a single copy integrated transgene and that developed with improved synchrony and speed ( Fig. S2D-J, Methods). Our results confirmed that L1 was greatly extended relative to the other 20 larval stages (Fig. S2I). However, in contrast to the previous findings (Olmedo et al., 2015), but consistent with our hypothesis, the differences appeared largely attributable to an extended intermolt (Fig. S2G). The duration of the first molt (M1) was instead comparable to that of M2 and M3 (Fig. S2H).
Thus, an extended first intermolt coincides with the fact that no oscillator activity can be detected by RNA sequencing during the first five hours of this larval stage. Moreover, because we performed the experiment by hatching embryos directly into food, we can conclude that the 5 extended L1 stage is an inherent feature of C. elegans larval development, rather than a consequence of starvation-induced synchronization.

Development is coupled to oscillatory gene expression
The luciferase assay revealed that also the L4 stage took significantly longer than the two 10 preceding stages, though not as long as L1 (Fig. S2I). In this case, both the fourth intermolt and the fourth molt were extended (Fig. S2G, H). As apparent from the gene expression heatmap, and quantified below, the oscillation period during L4 was also extended. Hence, grossly similar trends appeared to occur in larval stage durations and oscillation periods, determined by the luciferase assay and RNA sequencing, respectively. We considered this as further evidence for a 15 coupling of the two processes.
To test this hypothesis explicitly, we sought to quantify the synchrony of oscillatory gene expression and developmental progression in individual animals at the same time. To this end, we established a microchamber-based time-lapse microscopy assay by adapting a previous protocol (Turek et al., 2015). In this assay, animals are hatched and grown individually in small 20 chambers where they can be tracked and imaged while moving freely, enabling their progression through molts. Using Mos1-mediated single copy transgene integration (MosSCI) (Frøkjaer-Jensen et al., 2012), we generated transgenic animals that expressed destabilized gfp from the promoter of qua-1, a highly expressed gene with a large mRNA level amplitude.
Consistent with the RNA sequencing data, we detected oscillations of the reporter with four expression peaks (Fig. 3B). Moreover, we observed similar rates of development as in the luciferase assays when we curated the molts (Fig. 3C, Table S2, Methods). The stage durations 5 were in good agreement with the averaged oscillation period times for each cycle, obtained through a Hilbert transform of GFP intensities, for all three larval stages, L2 through L4, for which oscillations period lengths could be reliably determined (Fig. 3D).
Single animal imaging enabled us to ask when molts occurred relative to oscillatory gene expression, and we observed a very uniform behavior across animals (Fig. 3B). To quantify this 10 relationship, we determined the gene expression phases at molt entries and exits. We obtained highly similar values across worms within one larval stage (Fig. 3E), and only a minor drift when comparing phases across larval stages. Two additional reporter transgenes, based on the promoters of dpy-9 and F11E6.3, which differ in peak expression phases from qua-1 and one another, yielded similar results (Fig. S3). 15 We considered two possible interpretations of the narrow distributions of oscillation phases at molt entry and exit: first, both oscillations and development could be under independent, but precise temporal control. In this model, certain developmental events would merely coincide with specific phases of oscillations rather than being coupled to them. Therefore, variations in the periods of oscillation and development would add up, non-linearly, 20 to the experimentally observed phase variations. Second, phase-locking of oscillatory gene expression and developmental events might result from the two processes being truly coupled and/or from one driving the other. In this case, the variations in the two periods would partially explain each other, causing a reduction in the expected phase variation relative to the first scenario (Fig. 3F).
To distinguish between these scenarios, we used error-propagation to calculate the expected error for two independent processes (Methods). Focusing on L2 and L3 stages to exclude any edge effects on period calculation by Hilbert transform, we found that this 5 calculated error consistently exceeded the experimentally observed error (Fig. 3G), for all three reporter genes, for both molt entry and molt exit, and for both larval stages. Thus, our observations agree with the notion that development and oscillatory gene expression are functionally coupled (Fig. 3H), and potentially causally connected. 10 Quantification of amplitude and period behavior over time reveal characteristic systems properties Consistent with the coupling between oscillations and development, both the last larval stage and the period of the last oscillation cycle appeared increased (Fig. 3D), before oscillations ceased.
The characteristics of such a transition from oscillatory to non-oscillatory state, or bifurcation, 15 can be examined in the light of bifurcation theory. Bifurcation, that is, a qualitative change in system behavior, occurs in response to a change in one or more control parameters. Depending on the system's topology, characteristic changes of amplitude and period occur during bifurcation ( Fig. 4A) (Izhikevich, 2000;Salvi et al., 2016;Strogatz, 2015). Thus, transition into a quiescent state through a supercritical Hopf (supH) bifurcation involves a declining, and 20 ultimately undetectable, amplitude and a constant period. By contrast, a Saddle Node on Invariant Cycle (SNIC) bifurcation results in a declining frequency (and thus increasing period) but a stable amplitude.
Hence, to gain a better understanding of the oscillator's bifurcation, we quantified oscillation amplitudes and periods over time. To minimize variations from differences between experiments, we did this for the contiguous 5 -48 hour time course (TC2). This enabled reliable quantification of these features for the last three oscillation cycles, C2 through C4, which begin at 14 h (C2), 20 h (C3) and 27 h (C4), respectively (Fig. 4C). Excluding a small set of 291 genes 5 that exhibited unusual expression trends during the fourth larval stage, i.e., a major change in mean expression levels ( Fig. S4) this analysis revealed a good agreement of amplitudes across stages, and in particular no indication of damping during the last cycle, C4 (Fig. 4B).
We used a Hilbert transform (Pikovsky et al., 2001) to quantify the period over time with high temporal resolution, i.e., at every hour of development. The mean oscillation period thus 10 calculated was approximately seven hours during the first cycles but increased during the fourth cycle (Fig. 4C). This change was also apparent when we reconstructed an oscillation from the mean observed oscillation period and compared it to an oscillation with a constant period of seven hours (Fig. 4D).
In summary, these analyses reveal a sudden loss of oscillation upon transition to 15 adulthood without prior amplitude damping and an oscillator that can maintain a stable amplitude in the presence of period changes. These are features of an oscillator operating near a SNIC rather than a supH bifurcation ( Fig. 4A) (Izhikevich, 2000;Salvi et al., 2016;Strogatz, 2015). 20 Arrest of the oscillator in a specific phase upon transition to adulthood SNIC and supH bifurcations differ not only in amplitude and period behavior, but also in the stable state, or fixed point, that the systems adopts. In a supH bifurcation, the system spirals from a limit cycle onto a fixed point, whereas in a SNIC bifurcation, the fixed point emerges on the limit cycle ( Fig. 5A) (Saggio et al., 2017). In other words, a quiescent oscillator near a SNIC bifurcation adopts a state similar to that of a specific phase of the oscillator; the oscillator has become 'arrested'. By contrast, following a supH bifurcation, the oscillator adopts a stable state that is distinct from any phase of the oscillator. Hence, if the C. elegans oscillator entered an 5 arrested state through a SNIC bifurcation, the overall expression profile of the oscillating genes in the adult stage should resemble that seen at some other time point during larval development.
To test this prediction, we analyzed the correlation of oscillating gene expression for adult time points (TP ≥ 37 h) to all other time points of the fused time course (TC3). (In the following, we will use "TPx" to refer to any time point 'x', in hours, after hatching. Technically, 10 this is defined in our experiment as the time after plating synchronized, first larval stage animals on food.) For this analysis (illustrated in Fig. S5), we used the pairwise correlation matrix resulting from the oscillating gene set without the previously excluded genes that changed in expression in the L4 stage ( Fig 5B). This provided two insights. First, correlation coefficients among adult time points all exceeded 0.8 with little change over time, confirming the high 15 similarity of samples TP37 -48 to one another and thus an absence of detectable oscillations.
Second, in addition to one another, TP37 -48 are particularly highly correlated to a specific time and thus phase -of the oscillatory regime, namely TP13 and the 'repetitive' TP19 and TP26/27 ( Fig. 5C, arrows). In other words, expression levels of oscillating genes in the adult resembled a specific larval oscillator phase, providing further support for a SNIC bifurcation.

Phase-specific arrest of the oscillator after hatching
We noticed that the gene expression states of TP37 -48 also correlated well to each of TP1 -5; i.e., the early L1 larval stage (Fig. 5B). To examine this further, we performed the same correlation analysis as described above, but now for TP1 -5. Mirroring the adult situation, correlation coefficients among all these five time points were high and exhibited little change 5 over time, and TP1 -5 exhibited particularly high levels of correlation to TP13 and the 'repetitive' TP19 and TP27 (Fig. 5D). These are the same larval time points to which the adult time points exhibit maximum similarity. We confirmed these two key observations when fusing the independent time course TC4 to TC2 (generating TC5; Fig. S6).
We conclude that also during the first five hours after plating, oscillating genes adopt a 10 stable expression profile that resembles a specific phase of the oscillator. In other words, both the transition into the oscillatory state during L1 and out of it during L4 occur through a SNIC bifurcation. This finding indeed explains our observation ( Fig. 2) that in L1 stage larvae, oscillations exhibit a structure of phase-locked gene expressing patterns as soon as they become detectable: the oscillator initiates from an arrested phase. 15

Initiation of oscillation soon after gastrulation
We wondered how the oscillator entered the arrested state observed in early larvae, i.e., what dynamics the class of larval oscillating genes exhibited in embryos. Hence, we examined single embryo gene expression data from a published time series (Hashimshony et al., 2015). When 20 plotting the embryonic expression patterns of oscillating genes sorted by their peak phase defined in larvae, we observed a dynamic expression pattern with a striking phase signature (

A shared oscillator phase for experimentally induced and naturally occurring bifurcations
The arrested states of the oscillator in both early L1 stage larvae and in adults are highly similar and resemble the oscillator state at TP19/300°. Therefore, we wondered whether this oscillator phase was particularly conducive to state transitions of the system in response to changes in the 20 developmental trajectory. To test this, we examined animals that exited from dauer arrest, a diapause stage that animals enter during larval development under conditions of environmental stress such as heat, crowding, or food shortage. Using a published time course of animals released from dauer arrest after starvation (Hendriks et al., 2014), we found that their expression patterns of oscillating genes correlated highly with those of animals initiating oscillations (TC3) in the L1 stage (Fig. 7A, B). Additionally, gene expression patterns at 1 hour through 5 hours and at 13 hours post-dauer were highly correlated to those of the repetitive TP13, TP19 and TP26/27 during continuous development. Hence, the system state during a period of quiescence 5 during the first five hours after placing animals on food corresponds to an arrest of the oscillation in a phase seen at TP19/300° of the continuous development time course.
We conclude that the system bifurcates in the same manner during continuous, unperturbed development, after hatching, and in response to a perturbation, namely starvationinduced dauer arrest.

Discussion
In this study, we have characterized biological function and behavior of the C. elegans larval oscillator. Our results from single animal-and population-based analyses reveal a close coupling to development, and specifically molting, and imply that processes essential for molting may not 15 be restricted to lethargus. We have observed that oscillations are highly similar during the four cycles ( Fig. 7C, Fig. S8). Yet, oscillations cease and (re-) initiate several times during physiological development, and similar state transitions, or bifurcations, of the system can be induced through an external perturbation (Fig. 7C). In particular, all non-oscillatory states correspond to an arrest of the oscillator in one specific phase. Hence, the observed bifurcations 20 provide a conceptual model of how a developmental checkpoint can operate to halt larval development at a particular, repetitive point of development. Moreover, they constrain possible system architectures and properties as well as the choice and parametrization of mathematical models that can represent the system.

Oscillatory expression of thousands of genes
Our previous work (Hendriks et al., 2014) identified ~2,700 oscillating (i.e., rhythmically 5 expressing) genes, a number that we now increase to 3,739 genes (24% of total expressed genes).
We attribute this improved identification of oscillating genes to the fact that our present analysis focused on the L1, L2 and L3 stages, where a constant oscillation period of ~7 hours facilitates cosine wave fitting. This contrasts with the situation in the previous experiment, which used data from the L3 and L4 stages and thus, as we reveal here, a time of changing period. 10 Even our current estimate is conservative, i.e., the 'non-oscillating' genes contain genes that exhibit oscillatory expression with low amplitude or, potentially, strongly non-sinusoidal shapes. It is possible that such dynamics may play important roles for specific genes and processes and our data provide a resource to identify these in the future. However, here we focused on genes with robust and extensive oscillations to facilitate functional dissection of the 15 oscillator.

A developmental oscillator with functions in and beyond molting and lethargus
The physiological function of the C. elegans oscillator has remained unclear. Here, we have tested a possible connection with molting. By quantifying the periods of both molting and gene 20 expression oscillations simultaneously, in the same individual animals, we revealed their high degree of similarity and showed that the two processes are more closely phase-locked than expected from mere coincidence, i.e., they are coupled. We propose that a function of the oscillator as a developmental clock provides a parsimonious explanation for the coupling, but other models remain possible, e.g., the oscillator may facilitate an efficient molting process by anticipating the time of peak demand for cuticular building blocks and other factors.
Conventionally, molting is subdivided into three distinct steps, namely apolysis (severing of connections between the cuticle and the underlying epidermis), new cuticle synthesis, and 5 ecdysis (cuticle shedding) (Lažetić and Fay, 2017). The first two occur during, and the latter terminates, lethargus, a period of behavioral quiescence. However, we find that the temporal structure imposed by the C. elegans oscillator extends beyond lethargus. Our data reveal two probable causes: occurrence of processes required for molting before lethargus and a temporal organization that extends to processes unrelated to molting. 10 Specifically, we observed initiation of oscillations in embryos, which execute cuticle synthesis but neither apolysis nor ecdysis, at ~380 min into embryo development and thus long before the first signs of cuticle synthesis at ~600 min (Sulston et al., 1983). Instead, this time coincides with formation of an apical extracellular matrix (ECM). Although termed embryonic sheath, we find genes encoding components of this ECM, namely sym-1; fbn-1; noah-1; noah-2 15 (Vuong-Brender et al., 2017), also detectably expressed in larvae, and all four are required for larval molting or proper cuticle formation (Frand et al., 2005;Niwa et al., 2009). Moreover, their mRNA levels oscillate with high amplitudes, and as predicted from the embryonic data, their expression peaks long before lethargus, and in fact shortly after ecdysis, i.e., after a molt has been completed. Hence, to account for all these facts, we propose that molting involve processes 20 that are executed long before the onset of lethargus and that include ECM remodeling.
However, even a substantially more complex molting process may fail to account for the fact that a majority of oscillating genes is phase-locked with the molt but exhibits peak expression outside the molt and lacks any obvious link to molting. We consider it plausible that robust larval development may benefit from a coordination of molting with other physiological or developmental processes, as previously postulated for skin cell proliferation (Ruaud and Bessereau, 2006). Similarly, as lethargus involves cessation of food-uptake, oscillatory gene expression may serve to anticipate this event, consistent with a large number of intestinally 5 expressed oscillating genes. Nonetheless, even for this class, a broad dispersion of peak expression phases may suggest additional functions, yet to be uncovered. Whatever the benefit, it is evident that the oscillator imposes a temporal structure of gene expression that extends far beyond lethargus. 10

Oscillatory state transitions and developmental checkpoints
We have observed a loss of oscillations under three distinct conditions, in early L1 stage larvae, dauer arrested animals, and adults. The similarity of the oscillator states under all three conditions is striking and involves an arrest in the same specific phase.
Formally, for the L1 arrest, we cannot distinguish between perturbation-induced or 15 naturally occurring arrest, as the sequencing experiments required animal synchronization by hatching animals in the absence of food, causing a transient arrest of development. However, the fact that the L1 stage is extended also in animals hatched into food suggests that they may adopt a similar arrested state even in the presence of food, perhaps because the nutritional resources in the egg (i.e., egg yolk) have become depleted by the time that hatching occurs. In other words, 20 synchronization of L1 animals by hatching them in the absence of food may propagate a preexisting transient developmental and oscillator arrest.
Irrespective of this interpretation, a key feature of the arrests that we observe under different conditions is that they always occur in the same phase. This is a behavior one would predict for a repetitive developmental checkpoint. Such a checkpoint has indeed been found to operate shortly after each larval molt exit, arresting development in response to a lack of food (Schindler et al., 2014). Importantly, developmental arrest does not result from an acute shortage 5 of resources. Rather, it is a genetically encoded, presumably adaptive, response to nutritionally poor conditions, as demonstrated by the fact that mutations in the daf-2/IGFR signaling pathway causes animals to continue development under the same food-depleted conditions (Baugh, 2013;Schindler et al., 2014).
Within the limits of our resolution, the phase of the arrested oscillator corresponds to the 10 phase seen around ecdysis. Hence, oscillations and development are synchronously arrested, and we propose that signals related to food sensing, metabolism, or nutritional state of the animal help to control the state of the oscillatory system and thereby developmental progression. An oscillator operating near a SNIC bifurcation appears ideally suited to processing such information, because it acts as a signal integrator, i.e., it becomes active when a signal threshold 15 is surpassed (Forger, 2017;Izhikevich, 2000). This contrasts with the behavior of oscillators operating near a supH bifurcation, which function as resonators, i.e., they respond most strongly to an incoming signal of a preferred frequency. Hence, both the phase-specific arrest and the integrator function as characteristics of an oscillator operating in the vicinity of a SNIC bifurcation are physiologically relevant features of this C. elegans oscillator. 20

Insights into oscillator architecture and constraints for mathematical modelling
What mechanisms determine the type and behavior of the C. elegans oscillator? Although the nature and wiring of the 'core oscillator', i.e., the machinery that drives the pervasive gene expression oscillation, remains to be established, the behavior of the oscillator that we characterized here provides clear constraints. Thus, a change in period without a noticeable change in amplitude, as seen in L4 stage larvae, is a feature of rigid oscillators (Abraham et al., An exception are so-called amplified negative feedback oscillators, which operate near a SNIC 10 bifurcation and rely on interlinked negative and positive feedback loops. Beyond constraining possible oscillator architectures, our experimental observations will help to guide mathematical modelling of the C. elegans oscillator. Modelling is needed because the nonlinear dynamic behaviors of oscillators are difficult to grasp intuitively. However, it is usually difficult to ensure the relevance of a given model, because both its formulation and its 15 parametrization determine whether oscillations occur and which behaviors the resulting oscillator model displays. Amplified negative feedback oscillators are a case in point as they can also operate near a supH bifurcation; operation near a SNIC bifurcation occurs only in a certain parameter space (Conrad et al., 2008;Guantes and Poyatos, 2006). The experimental characterization of the system's bifurcation that we have achieved here will therefore provide 20 crucial constraints to exclude invalid mathematical models.
We do note that mathematical models of somitogenesis clocks, inspired by mechanistic knowledge about the identity of individual oscillator components and their wiring, tend to represent oscillators operating near a supH bifurcation (Jensen et al., 2010;Webb et al., 2016). (Webb et al., 2016). At the same time it contrasts with changes in both amplitude and period that were shown to occur in zebrafish embryos during somite formation and prior to cessation of oscillation (Shih et al., 2015). Thus, and because an analysis of bifurcation behavior of somitogenesis clocks in vivo is 5 challenging due to a complex space-dependence of oscillation features (Soroldoni et al., 2014), it remains to be answered whether and to what extent the C. elegans oscillator and the somitogenesis clocks share specific properties. However, whatever the answer, a comparison of the similarities and differences in behaviors, architectures and topologies will help to reveal whether and to what extent diverse developmental oscillators follow common design principles.

Conflict of Interests:
The authors declare that they have no conflict of interest. (B) Scatter plot identifying genes with oscillatory expression (henceforth termed oscillating genes, blue) based on amplitude and 99% confidence interval (99%-CI) of a cosine fitting of their expression quantified in TC2 (Methods). A lower CI-boundary ≥ 0, i.e. p-value ≤ 0.01, and a 5 log2(amplitude) ≥ 0.5, which corresponds to a 2-fold change from peak to trough, were used as cut-offs. Genes below either cut-off were included in the 'not oscillating' group (black). Fig. S1C reveals gene distributions in a density scatter plot.
(C) Gene expression heatmaps of oscillating genes as classified in Fig. 1B and S1C. Oscillating  (D) Scatter plot of calculated oscillating gene peak phase (Fig. 1) over the time of occurrence of the first expression peak in L1 larvae, observed in TC4. Peak detection was performed using a 10 spline analysis. As visual inspection did not reveal peaks in the heatmap during the first three hours, and as the first cycle ends at 13 h, we performed this analysis for t = 3 h to t = 13 h to reduce noise.
(E) Scatterplot comparing experimentally identified first peaks of gene expression (as in D) of the two early time course replicates, TC1 and TC4. 15 All analyses for oscillating genes identified in Fig. 1 with detectable expression (n=3,680).  Horizontal gray bars indicate oscillation cycles C1 through C4 as in Fig. 1C.  (C) Schematic depiction of behavior of the C. elegans oscillator from embryo to adult. A phasespecific arrest (red dot) is observed in hatch, early L1, young adults, and dauer-arrested animals.
See Fig. S8 for additional data supporting four similar oscillation cycles during L1 through L4.

C. elegans strains
The Bristol N2 strain was used as wild type. The following transgenic strains were used: inserting them into NheI-digested pYPH0.14 or pMM001. PCR primers and resulting plasmids are listed in the Table S3. Third, we obtained transgenic worms by single-copy integration into EG8079 worms, containing the universal ttTi5605 locus on chromosome II by following the  (Fig. 2), RNA-seq data were mapped to the C. elegans ce10 genome using STAR with default parameters (version 2.7.0f) and reads were counted using htseq-count (version = 0.11.2).
This step was omitted in the early time courses because many genes start robust expressing only As the density of the genes strongly decreased around 0.5 (Fig. S1C) we used amplitude ≥ 0.5 as a first classifier. We propagated the standard error of the coefficients A and B to the amplitude using Taylor expansion in the 'propagate' function (expr=expression(sqrt(((A^2)+(B^2))), ntype = 'stat', do.sim=FALSE, alpha=0.01) from the package 'propagate' (version 1.0-6) (Spiess, 2018) in R. We obtained a 99% confidence interval (99%-CI) for each gene. As 99%-CI that 5 does not include 0 is significant (p-value=0.01), we used the lower boundary (0.5%) of the CI as a second classifier. Thus, we classified genes with an amplitude ≥ 0.5 and lower CI-boundary ≥ 0 as 'oscillating' and genes with an amplitude < 0.5 or a lower CI-boundary < 0 were classified as 'not oscillating' (Fig. 1B, S1C). Every gene thus has an amplitude and a peak phase. A peak phase of 0° is arbitrarily chosen, and thus current peak phases are expected to differ 10 systematically from the previously assigned peak phases (Hendriks et al., 2014). TP20-TP27 and C4: TP27-TP36 developmental time. We applied cosine fitting to C2, C3 and C4 as described above to the expression of oscillating genes in TC2, excluding genes with deviating 15 mean expression in L4 as described above. We excluded C1, because amplitudes were sometimes difficult to call reliably. We used a fixed period of 7 h for C2-C3 and 8.5 h for C4 as determined by quantification of the oscillation period (Fig. 4C). We applied a linear regression using the function 'lm' of the package 'stats' in R to find the relationship between the amplitudes across different stages, i.e. the slope. The correlation coefficient, r, was determined using the 20 'cor' function (method=pearson) of the package 'stats' in R.

Quantification of oscillation period
For a temporally resolved quantification of the oscillation period, we filtered the mean- identified the correlation of particular time points (e.g. TP14-TP19) with their corresponding time points in the next oscillation cycle. Thereby, we were able to identify cycle time points as described in the results section. We defined the first cycle time point, e.g. TP14 of cycle 2, as 0°, and the last unique one, TP19, as 300°. TP14 (0° of cycle 2) is also 360° of cycle 1. Note that a sampling interval of 1 hour means that a TP in one cycle may correlate equally well to two adjacent TPs in another cycle, as seen for instance in the correlation of TP13 to TP26 and TP27. 5 The spline interpolation places the peak of correlation in the middle of these time points at ~TP26.5. The spline analysis thus annotates cycle points correctly even in C4 which has an extended period.
We performed correlation analyses without mean normalization of expression data, hence correlation values cannot be negative but remain between 0 and 1. We made this decision 10 because a correlation analysis using mean-centered data, where correlations can vary between -1 and +1, requires specific assumptions on which time points to include or exclude for mean normalization, and because it is sensitive to gene expression trends. However, we confirmed, as a proof of principle, the expected negative correlation of time points that are in antiphase when using mean-centered data (Fig. S9)  gene's highest and second highest tissue expression and the ratio of is reported. We selected tissue specific genes based on a ratio > 5 and a qvalue < 0.05 (these criteria reduced the number of genes to investigate). Using this list of genes we calculated the percentage of tissues present in both, all genes and oscillating genes using the function "Counter" from "collections" in python (labels, values = zip(*Counter(tissue_info_thr["max.tissue"]).items()) ). In order to obtain the 30 enrichment of tissues, we divided the percentage of tissue X among oscillating genes in the tissue enriched dataset by the percentage of tissue X among all genes in the tissue enriched data set and plotted the resulting values. The list of tissue specific oscillating genes was further used to investigate the peak phases within one tissue by plotting a density plot of the peak phase (from Fig. 1) for every tissue. As we lack data below 0 degree and above 360 degree, density values at these borders are distorted as the density is calculated over a moving window. Since we are 5 confronted with cyclical data and thus 0 degree corresponds to 360 degree, we added and subtracted 360 degree to each phase value, thus creating data that ranged from -360 degree to 720 degree which allowed us to plot the correct density at the borders 0 and 360 degree. We used python (pandas) to plot this data using the following command: data_tissue ["Phase"].plot(kind="kde", linewidth=5, alpha=0.5, bw=0.1) 10 Identification of first gene expression peaks in L1 larvae To identify the first peak of oscillating genes, we used a spline analysis in Python ("from scipy.interpolate import InterpolatedUnivariateSpline") from TP3 -TP13. We chose these time points to remove false positives in the beginning due to slightly higher noise for the first 2 time 15 points as well as not to identify the second peak which occurred at ≥TP14 for some very early

Time-lapse imaging of single animals
Single worm imaging was done by adapting a previous protocol (Turek et al., 2015), and is 20 similar to the method reported in (Gritti et al., 2016). Specifically, we replaced the previous 3.5cm dishes with a "sandwich-like" system: The bottom consisted of a glass cover slip onto which two silicone isolators (GRACE Bio-Labs, SKU: 666103) with a hole in the middle were placed on top of each other and glued onto the glass cover slip. We then placed single eggs inside the single OP50 containing chambers, which were made of 4.5% agarose in S-basal. The chambers 25 including worms were then flipped 180 degree and placed onto the glass cover slip with the silicone isolators, so that worms faced the cover slip. Low melt agarose (3% in S-basal) was used to seal the agarose with the chambers to prevent drying out or drifts of the agarose chambers during imaging. The sandwich-like system was then covered with a glass slide on the top of the silicone isolators to close the system.
We used a 2x sCMOS camera model (T2) CSU_W1 Yokogawa microscope with 20x air objective, NA = 0.8 in combination with a 50µm disk unit to obtain images of single worms. For a high throughput, we motorized the stage positioning and the exchange between confocal and brightfield. We used a red LED light to combine brightfield with fluorescence without closing the shutter. Additionally, we used a motorized z-drive with 2 µm step size and 23 images per z-5 stack. The 488nm laser power for GFP imaging was set to 70% and a binning of 2 was used.
To facilitate detection of transgene expression and oscillation, we generated reporters using the promoters of genes that exhibited high transcript levels and amplitudes, and where GFP was concentrated in the nucleus and destabilized through fusion to PEST::H2B (see strain list above). We placed embryos into chambers containing food (concentrated bacteria HT115 with 10 L4440 vector) and imaged every worm with a z-stack in time intervals of 10 min during larval development in a room kept at ~21°C, using a double camera setting to acquire brightfield images in parallel with the fluorescent images. We exploited the availability of matching fluorescent and brightfield images to identify worms by machine learning. After identification, Time lapse images were analyzed using a customized KNIME workflow (Supplemental File 1). We analyzed every worm over time using the same algorithm. First, we identified the 20 brightest focal planes per time point by calculating the mean intensity from all focal planes per time point and selecting the focal planes that had a higher intensity than the mean. Then we maximum-projected the GFP images over Z per time point and blurred the DIC image and also max projected over Z (blurring the DIC improved the machine learning process later on). All images per worm over time were analyzed by Ilastik machine learning in order to identify the 25 worm in the image. The probability map from Ilastik was used to select a threshold that selected worms of a particular experiment best. (The threshold might change slightly as DIC images can look slightly different due to differences in the sample prep amongst experiments.) Using a customized ImageJ plugin, we straightened the worm. The straightened GFP worm image was then max projected over Y which resulted in a single pixel line representing the GFP intensities 30 in a worm and after stacking up all the single pixel lines in Y direction, we obtained the 47 kymographs. In order to remove noise coming from the head and tail regions of the worm due to inaccuracy of the machine learning, we measured mean GFP intensities per time point ranging from 20% until 80% of the worms anteriorposterior axis. For background subtraction we exploited the fact that only the nuclei were GFP positive and thus subtracted the minimum intensity value between GFP nuclei from their intensity values. 5 After the KNIME workflow, we imported the measured GFP intensities into Python and analyzed the traces using a butterworth filter and Hilbert transform analysis (both from Scipy (Jones et al., 2001)). We used the butterworth bandpass filter using b, a = butter(order =1, [low,high], btype="band") with low=1/14 and high=1/5, corresponding to 14 hour and 5 hour periods respectively. We then filtered using filtfilt(b, a, data, padtype='constant') to linearly filter 10 backwards and forwards.
For individual time points where the worm could not be identified by the Ilastik machine learning algorithm, we linearly interpolated (using interpolation from pandas (McKinney, 2010)) using "pandas.series.interpolate(method = 'linear', axis =0, limit = 60, limit_direction = 'backward'", between the neighboring time points to obtain a continuous time series needed for 15 the Hilbert transform analysis. Using Hilbert transform, we extracted the phase of the oscillating traces for each time point and specifically investigated the phase at molt entry and molt exit for our different reporter strains. In order to determine time points in which worms are in lethargus, we investigated pumping behavior. As the z-stack of an individual time point gives a short representation of a 20 moving worm, it is possible to determine whether animals pump (feeding, corresponds to intermolt) or not (lethargus / molt). Additionally to the pumping behavior, we used two further requirements that needed to be true in order to assign the lethargus time span: First, worms needed to be quiescent (not moving, and straight line) and second, a cuticle needed to be shed at the end of lethargus. Usually worms start pumping one to two time points before they shed the 25 cuticle. This analysis was done manually with the software ImageJ, and results were recorded in an excel file, where for every time point, the worms' behavior was denoted as 1 for pumping and as 0 for non-pumping.
To determine a possible connection between oscillations and development, we applied error propagation, assuming normal distribution of the measured phases and larval stage 30 durations. Thereby, we exploited the inherent variation of the oscillation periods and developmental rates among worms, rather than experimental perturbation, to probe for such a connection. We define the phase at either molt exit or entry as ≡ 2 * ~ ( , 2 ) with ~ (µ , 2 ) being the period of oscillation and ~ (µ , 2 ) the intermolt duration (for phase at molt entry) or larval stage duration (for phase at molt exit), resulting in a phase with mean and a standard deviation . Should the two processes be coupled as in scenario 2, we 5 would expect < .
To calculate the phase at molt entry and molt exit with error propagation we used the "uncertainties" package (Lebigot, n.d.) in python. The larval stage duration as well as intermolt duration and period were treated as ufloat numbers, representing the distributions coming from 10 our measurement (e.g. 7.5 +/-0.2). These distributions were then used to calculate the expected phase at molt entry (using the intermolt duration) and molt exit (using the larval stage duration) using: ℎ ( . ) = 2 * * . This resulted in the phase being represented by an ufloat number and thus a distribution which we used for plotting after normalizing for the mean to compare the variation of the data.