Precise temporal control of neuroblast migration through combined regulation and feedback of a Wnt receptor

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

587 views

Version of Record published: June 7, 2023 (This version)
Accepted Manuscript published: May 15, 2023 (Go to version)
Accepted: May 12, 2023
Preprint posted: October 26, 2022 (Go to version)
Received: August 12, 2022

1. Of interest
Beta-cell intrinsic dynamics rather than gap junction structure dictates subpopulations in the islet functional network

Jennifer K Briggs, Anne Gresch ... Richard KP Benninger

Research Article Nov 29, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

Many developmental processes depend on precise temporal control of gene expression. We have previously established a theoretical framework for regulatory strategies that can govern such high temporal precision, but experimental validation of these predictions was still lacking. Here, we use the time-dependent expression of a Wnt receptor that controls neuroblast migration in Caenorhabditis elegans as a tractable system to study a robust, cell-intrinsic timing mechanism in vivo. Single-molecule mRNA quantification showed that the expression of the receptor increases non-linearly, a dynamic that is predicted to enhance timing precision over an unregulated, linear increase in timekeeper abundance. We show that this upregulation depends on transcriptional activation, providing in vivo evidence for a model in which the timing of receptor expression is regulated through an accumulating activator that triggers expression when a specific threshold is reached. This timing mechanism acts across a cell division that occurs in the neuroblast lineage and is influenced by the asymmetry of the division. Finally, we show that positive feedback of receptor expression through the canonical Wnt pathway enhances temporal precision. We conclude that robust cell-intrinsic timing can be achieved by combining regulation and feedback of the timekeeper gene.

Editor's evaluation

This paper deals with an important unsolved problem in developmental biology: how cells execute their dynamics at the right time. The study combines compelling quantitative single-cell and single-transcript experiments with genetic perturbations and computational modelling and provides important insights into how the timing of transcription is regulated.

https://doi.org/10.7554/eLife.82675.sa0

Introduction

Timing plays a central role in development, coordinating processes ranging from cell division (Tyson and Novak, 2008) and differentiation (Ray et al., 2022), to the complex segmentation (Dequéant and Pourquié, 2008) and limb development of vertebrate embryos (Pickering et al., 2018). To measure time, biological clocks generally utilize a component that increases or decreases in activity or abundance (a timer) or cycles through a high and low state (an oscillator) and triggers a response when a specific threshold is crossed (Gliech and Holland, 2020). An example of such a biological clock is the time-dependent expression of Neuropilin-1 in retinal ganglion cells (RGCs), which influences the trajectory of the outgrowing RGC axons (Baudet et al., 2011; Campbell et al., 2001). Here, the clock mechanism consists of a transcriptional repressor (CoREST) that gradually decreases in abundance as a result of miRNA dependent inhibition and releases Neuropilin-1 expression once it reaches a specific threshold. With its dependence on protein production, degradation and molecular interactions between proteins, such a timekeeping mechanism is inherently noisy, raising questions on the regulatory strategies that biological clocks have evolved to increase precision.

In clocks that act across multiple cells, robustness can be increased through averaging, synchronization through intercellular communication (Dequéant and Pourquié, 2008; Oates, 2020) or entrainment using external signals (Aydogan et al., 2020; Bell-Pedersen et al., 2005). The cell intrinsic mechanisms that control the timekeeping mechanism itself are, however, still poorly understood. Mathematical modeling has provided insight into regulatory circuits that can enhance timing precision. One mechanism for temporal regulation is the constant production of a stable timer molecule, which leads to a gradual accumulation and activation of the time-dependent response once a specific abundance threshold is reached. Modeling has shown that such a linear increase is more precise than a strategy where the timer molecule positively regulates its own expression, as such autoregulation also amplifies the noise (Ghusinga et al., 2017). Using a first-passage time approach, we have recently shown that a regulated strategy where the timekeeping molecule increases non-linearly over time can increase precision (Gupta et al., 2018). Such regulation can be mediated through an accumulating transcriptional activator or a gradually decreasing transcriptional repressor that controls the expression of the timekeeper molecule in a threshold-dependent manner, and we showed this to be robust to additional effects such as cell division and bursts in transcription (Gupta et al., 2018). Moreover, in contrast to autoregulation on its own, we found that precision could be further increased when the regulator was combined with autoregulation of the timekeeper (Gupta et al., 2020). It is not known, however, if this predicted regulatory strategy is utilized by biological clocks in vivo.

The migration of the Caenorhabditis elegans QR neuroblast descendants is controlled in a time-dependent fashion (Mentink et al., 2014) and provides a tractable system to study a biological clock at the single-cell level in vivo. QR is born in the mid-body as a sister cell of the hypodermal seam cell V5 and divides at the beginning of the first larval stage into an anterior daughter cell called QR.a and a posterior daughter QR.p (Figure 1A). Both daughter cells migrate a relatively long distance toward the anterior and divide at highly stereotypic positions to generate descendants that differentiate into neurons (Sulston and Horvitz, 1977). While QR.a divides only once to generate QR.ap and an apoptotic sister cell QR.aa, QR.p first divides into QR.pa and an apoptotic sister cell QR.pp. QR.pa continues with a short-range migration toward the anterior but then stops to divide into QR.paa and QR.pap, which migrate dorsally and ventrally to occupy their final positions, where they differentiate into a mechanosensory neuron and an interneuron, respectively. Compared to other neurons that undergo long-range migration in C. elegans, the final position of the QR.p descendants was found to be robust to stochastic variation (Dubois et al., 2021).

Figure 1 with 1 supplement see all

Download asset Open asset

The time-dependent expression of *mig-1* is independent of QR.p division.

(A) Schematic representation of the QR lineage and the migration of the QR descendants. V1–V6 are seam cells, used as landmarks to determine the position of QR and its descendants. For clarity, the QR.a branch of the QR lineage is not represented. (B) Single-molecule FISH (smFISH) quantification of *mig-1* expression in QR and its descendants relative to the position of the seam cells H2 to V5. n=200 from 11 independent experiments, with three replicates each. (C) DNA content in division-blocked QR.p (QR.p*), relative to 2C (Vn.p) seam cells. The DNA content of QR.p* does not significantly deviate from 4C (p=0.31, n=32, from five independent experiments, with three replicates each), indicating the cells do not undergo an additional S-phase after the division is blocked. Whiskers represent 1.5x interquartile range. (D) Representative images of *mig-1* smFISH spots in *cdk-1*(*hu277*[*AID::cdk-1*]). Top: control without auxin and bottom: QR.p division is blocked in the presence of auxin. Scale bar is 10 µm. Green Fluorescent Protein (GFP) is expressed in the seam cells and the Q neuroblast descendants (*heIs63*; Wildwater et al., 2011). (E) smFISH quantification of *mig-1* expression in QR and its descendants relative to the position of the seam cells H2 to V5 in *cdk-1*(*hu277*[*AID::cdk-1*]) in the absence (left) and presence (right) of auxin. Data are from seven independent experiments, with three replicates each.

The anterior migration of QR.p and QR.pa is dependent on Wnt signaling (Mentink et al., 2014; Rella et al., 2021). Wnt proteins are an evolutionarily conserved family of extracellular ligands that can signal through a canonical, β-catenin-dependent pathway to induce the expression of specific target genes, or through β-catenin-independent, non-canonical pathways that directly influence cytoskeletal dynamics (Angers and Moon, 2009; Clevers and Nusse, 2012). We have previously shown that the anterior migration of QR.p and QR.pa is mediated through parallel acting non-canonical Wnt pathways that separately control speed and polarity (Mentink et al., 2014). Once QR.pa reaches its final position, migration is stopped through activation of canonical, BAR-1/β-catenin-dependent Wnt signaling, which induces the expression of target genes that inhibit migration (Rella et al., 2021). This switch in signaling response is mediated through the Wnt receptor MIG-1/Frizzled (Mentink et al., 2014). Quantitative in vivo measurement of mig-1 expression using single-molecule FISH (smFISH) showed that mig-1 expression is rapidly upregulated in QR.pa, and transgenic rescue experiments demonstrated that mig-1 is both necessary and sufficient for migration termination (Mentink et al., 2014). Importantly, we found that this upregulation is not dependent on positional information but that the expression of mig-1 is temporally regulated. Thus, in mutants in which the QR descendants migrate at a slower speed, mig-1 is still expressed at the same time but at a more posterior position, while in animals with faster migrating QR descendants, mig-1 is expressed at a more anterior position. These results demonstrate that mig-1 expression is regulated through a timing mechanism and that the final position of the QR descendants is determined by the speed of migration and the time at which mig-1 expression is induced. Consistent with such a timing mechanism, the final QR descendants localize to a more anterior position in mutants with a shorter body size and at a more posterior position in animals with a longer body size, although it should be noted that some compensation takes place at the level of migration speed (Dubois et al., 2021).

Here, we combine mutant analysis and mathematical modeling to examine the timekeeping mechanism that controls mig-1 expression in the QR lineage. We show that the rapid upregulation of mig-1 expression in QR.pa is dependent on transcriptional activation, supporting a model in which the non-linear increase in mig-1 expression is mediated through threshold-crossing of an accumulating transcriptional activator. We find that this timing mechanism acts across the division of QR.p but is influenced by the asymmetry of the division. Finally, we show that positive feedback of mig-1 expression through the canonical Wnt pathway decreases temporal variability. We conclude that robust timing of gene expression can be achieved by combining an accumulating transcriptional activator with feedback regulation of the timekeeper molecule.

Results

Temporal regulation of mig-1 expression in the QR lineage is independent of cell division and cell cycle progression

During the anterior migration of the QR descendants, mig-1 expression is low in QR.p but is strongly upregulated in its daughter QR.pa (Mentink et al., 2014; Figure 1B), where it triggers the canonical Wnt signaling response that is necessary to stop migration (Rella et al., 2021). The presence of a cell division prior to the upregulation of mig-1 expression – the timing of which is tightly controlled as part of the invariant cell lineage of C. elegans development (Sulston and Horvitz, 1977) – raises the question whether the two are mechanistically linked. To investigate whether the upregulation of mig-1 expression is dependent on the division of QR.p (or subsequent cell cycle reentry in QR.pa), we blocked QR.p mitosis by conditionally depleting the M-phase regulator CDK-1 using the auxin inducible protein degradation system (Zhang et al., 2015). In this system, proteins tagged with an auxin inducible degron (AID) sequence are degraded in the presence of auxin and the F-box protein TIR1 (Zhang et al., 2015). We endogenously tagged cdk-1 with the AID sequence (AID::CDK-1) using CRISPR/Cas9-mediated genome editing and specifically expressed TIR1 in the Q neuroblast lineage using the egl-17 promoter (Branda and Stern, 2000). Since continuous exposure to auxin would block all divisions in the QR lineage, auxin was only applied from the stage at which QR has completed its division (300–345 min after hatching). We found that this efficiently inhibited the subsequent division of QR.p and its sister cell QR.a, as judged by the absence of their descendants (the apoptotic QR.pp and QR.aa cells and the neuronal QR.paa, QR.pap, and QR.ap cells).

Studies in yeast and mammalian cells have shown that loss of CDK1 not only inhibits mitosis but also cell cycle reentry and DNA replication (Coudreuse and Nurse, 2010). To investigate if the cell cycle is similarly blocked in the undivided QR.p cells, we examined DNA content by measuring total 4’,6-diamidino-2-phenylindole (DAPI) fluorescence at a time point at which the daughter of QR.p (QR.pa) would normally have divided into the final descendants QR.paa and QR.pap (8–9 hr after hatching) and normalizing to total DAPI fluorescence of the Vn.p (seam) cell nuclei, which have two unreplicated sets of chromosomes (2C) at this stage in L1 larval development (Hedgecock and White, 1985) and have a similarly sized nucleus. When mitosis and subsequent cell cycle progression is blocked by CDK-1 depletion, the DNA content of the replicated but non-divided QR.p nucleus is expected to be 4C. However, if cycling and DNA replication continues, a DNA content of 8C or higher is expected. As shown in Figure 1C, we observed no significant deviation from a DNA content of 4C. These results show that CDK-1 depletion inhibits mitosis as well as cell cycle reentry in the undivided QR.p cells.

Next, we measured mig-1 expression in the AID::CDK-1 strain using smFISH (Ji et al., 2013). In the absence of auxin, mig-1 expression was not significantly different from animals containing wild type CDK-1 (Figure 1B, D and E). Importantly, however, we found that in the presence of auxin, upregulation of mig-1 expression in undivided QR.p cells occurs at a similar time and range of expression as in the control animals (Figure 1D and E). We conclude that the temporal regulation of mig-1 expression is independent of the cell cycle and QR.p division.

The early and late phases of mig-1 expression are independently regulated

The expression of mig-1 in the QR lineage occurs in two phases (Mentink et al., 2014; Figure 1B): an early phase in which it is expressed in the QR neuroblast, followed by a rapid decrease during QR polarization (Ji et al., 2013) and low expression in QR.p until mig-1 is upregulated again in QR.pa (the late phase). To gain insight into the transcriptional regulation of mig-1, we made the assumption that cis-regulatory elements that control the expression of mig-1 are under stabilizing selection and therefore conserved across species (Gordân et al., 2010; Harbison et al., 2004; Nitta et al., 2015). By comparing the upstream region and first intron of mig-1 in 23 Caenorhabditis species, we found eight conserved 25 bp motifs that are located in four closely linked pairs (Figure 2—figure supplement 1, Figure 3—figure supplement 1).

Two pairs of motifs are located in the first intron (motif pair A and B; Figure 2A, Figure 2—figure supplement 1). To test if these motifs are required for the expression of mig-1, we independently deleted each of the two pairs using CRISPR/Cas9-mediated genome editing and quantified the number of mig-1 transcripts in QR and its descendants using smFISH. While deletion of motif pair B had no effect on mig-1 expression, deletion of motif pair A strongly reduced the early expression of mig-1 in the QR neuroblast (Figure 2B and C). The late expression of mig-1 in QR.pa was, however, not affected by deletion of this motif pair. Moreover, using the final position of the QR.pa daughter cell QR.pap as a measure of total migration distance, we found that the final position of the QR descendants was not significantly different from control animals (Figure 2—figure supplement 2). This is consistent with the observation that only the late phase of mig-1 expression is required for correct QR descendant migration (Mentink et al., 2014). Taken together, these results show that intron motif pair A is necessary for the early but not the late expression of mig-1 and support the notion that the two phases are independently regulated.

Figure 2 with 3 supplements see all

Download asset Open asset

A conserved region in the first intron of *mig-1* controls the early but not the late phase of *mig-1* expression.

(A) Schematic representation of the two conserved intron regions that were deleted from the *mig-1* locus. Green rectangles are exons, the gray dashed line is the upstream region, and black solid lines are introns. The orange rectangles represent the deleted regions that contain intron motif pair A (Δintron A) and intron motif pair B (Δintron B). Data are from three independent experiments, with three replicates each. (B) Single-molecule FISH (smFISH) quantification of *mig-1* expression in QR and its descendants relative to the position of the seam cells H2 to V5 in control and the intron deletion mutants *hu295* (intron deletion A) and *hu299* (intron deletion B). AP, anteroposterior. (C) Quantification of *mig-1* smFISH spots in QR (***p<0.0001, n=30, Welch’s t-test). Whiskers represent 1.5x interquartile range. (D) Mathematical model for the *mig-1* dynamics $m (t)$ regulated by an increasing activator (left) or decreasing repressor (right), fit to the control data (see Materials and methods and Figure 2-figure supplement 3). The sequences of the intron motifs in 23 different *Caenorhabditis* species are in Figure 2—figure supplement 1. Figure 2—figure supplement 2 shows that deletion of the intron motif pairs does not significantly change the final position of QR.pap.

A dynamical model of the temporal regulation of mig-1 expression

Taking into account that the late phase of mig-1 expression is independent of cell division and the early phase of expression, we developed a dynamical model of a cell-intrinsic timer mechanism driving mig-1 transcription. Since the speed of QR.p migration – which crosses most of the distance covered by the QR descendants – is roughly constant (Mentink et al., 2014), we treated distance along the anteroposterior axis as proportional to time. To account for the decrease (QR), followed by the increase (QR.pa), in the mig-1 amount over time, we hypothesized that mig-1 mRNA is continuously degraded and upregulated by a component with its own dynamics: either a transcriptional activator that increases in time or a repressor that decreases in time (Gupta et al., 2018). Thus, the dynamics of mig-1 mRNA number $m$ are described by $d m / d t = - ν m + F (t)$ , where $ν$ is the degradation rate, and $F (t)$ accounts for the regulation. For the activator model, we set $F (t) = {α a}^{H} / (a^{H} + K^{H})$ , where the activator molecule number increases with rate $k$ according to $a (t) = k t$ . For the repressor model, we set $F (t) = {α K}^{H} / (r^{H} + K^{H})$ , where the repressor molecule number decreases with rate μ from initial value $r_{0}$ according to $r (t) = r_{0} e^{- μ t}$ . In both models, $α$ is the maximal production rate of mig-1, $H$ is the Hill coefficient, and $K$ is the half-maximal value of the regulator. Fitting either model to the mig-1 expression data (see Materials and methods) gave good agreement (Figure 2D, Figure 2—figure supplement 3). Consistent with the intron motif pair A deletion data, we found that lowering the initial mig-1 amount ( $m (0) = 0$ ) does not affect the later upregulation dynamics (Figure 2D). Therefore, the model captures the observation that early and late mig-1 expression dynamics are independently regulated.

The late phase of mig-1 expression is dependent on two conserved, upstream activating cis-regulatory elements

The increasing transcriptional activator model and the decreasing transcriptional repressor model make contrasting predictions if the regulator is removed. Removing the activator does not influence mig-1 expression at early timepoints and only prevents upregulation of mig-1 at later timepoints (Figure 3A). Conversely, removing the repressor causes premature upregulation at early timepoints (Figure 3A). To gain insight into the transcriptional regulation of the late phase of mig-1 expression and to test the predictions of the activator and repressor models, we focused on conserved motifs in the mig-1 upstream region.

Figure 3 with 2 supplements see all

Download asset Open asset

Two conserved regions upstream of *mig-1* control the late but not the early phase of *mig-1* expression.

(A) Mathematical model with decreasing activator production rates $k = [6.6,5.6,4.6,0]$ (left) or decreasing initial repressor numbers $r_{0} = [19,13,7, 0]$ (right). (B) Schematic representation of the two conserved regions in the *mig-1* upstream sequence. Green rectangles are exons, the gray dashed line is the upstream region, and black solid lines are introns. The orange rectangles represent the deleted regions that contain upstream motif pair #1 (Δupstream 1) and upstream motif pair #2 (Δupstream 2). (C) Single-molecule FISH (smFISH) quantification of *mig-1* expression in QR and its descendants relative to the position of the seam cells H2 to V5 in control and deletion mutants of the two regions upstream of *mig-1*. None of the deletions affect the early phase of *mig-1* expression in the QR neuroblast (*hu314* Δupstream 1), p=0.59, n=22; *hu315* (Δupstream 2), p=0.07, n=40; *hu335* (Δupstream 1–2), p=0.89, n=22, Welch’s t-test. The upregulation of *mig-1* expression in QR.pa is not significantly changed in Δupstream 1 (p=0.44, n=51) but is significantly reduced in Δupstream 2 (p<0.0001, n=83) and abolished in the double Δupstream 1–2 deletion (p<0.0001, n=29). Data are from eight independent experiments with three replicates for Δupstream 1 and 2 and three independent experiments with three replicates for the double Δupstream 1–2 deletion. The sequences of the upstream motifs in 23 different *Caenorhabditis* species are in Figure 3—figure supplement 1. Figure 3—figure supplement 2 shows the effect of the upstream motif pair deletions on the final position of QR.pap.

Two pairs of conserved, closely spaced motifs were identified in a 3 kb region upstream of the mig-1 locus (Figure 3B, Figure 3—figure supplement 1): a distal pair (#1) at 2781–2571 bp and a proximal pair (#2) at 267–137 bp from the start of the mig-1 coding sequence. Using CRISPR/Cas9-mediated genome editing, we individually deleted each of the motif pairs and also created a double mutant in which both pairs are deleted simultaneously. None of the deletions had an effect on the early phase of mig-1 expression in the QR neuroblast (Figure 3C). However, we found that in the double deletion mutant, the late phase of mig-1 expression in QR.pa was completely lost (Figure 3C). A similar but less penetrant reduction was observed in the single motif pair 2 deletion, while in the motif pair 1 deletion a minor fraction of QR.pa cells failed to upregulate mig-1 expression (Figure 3C), a difference that did not result in a statistically significant difference in the mean expression of mig-1. Consistent with the role of the late phase of mig-1 expression in terminating the anterior migration of QR.pa, the average position of QR.pa (p<0.0001 for the single and double deletions) and the total migration of the QR descendants (as determined from the final position of QR.pap) was shifted anteriorly (Figure 3—figure supplement 2). As expected, this effect on total migration distance was most pronounced in the double deletion and the single motif pair 2 deletion, but there was also significant overmigration in the motif pair 1 deletion mutant, indicating that the mild reduction in mig-1 expression observed in this mutant is still sufficient to affect the final position of the QR descendants. Taken together, these results demonstrate that the two motif pairs represent partially redundant cis-regulatory elements that are required for the late phase of mig-1 expression and provide further support for the conclusion that the early and late phases of the expression are independently regulated.

Importantly, the loss of the late phase of mig-1 expression – instead of premature expression earlier in the lineage – is in agreement with the activator model, not the repressor model, of mig-1 regulation (Figure 3A).

Asymmetry of the QR.p division is required for upregulation of mig-1 expression in QR.pa

The division of QR.p that precedes the upregulation of mig-1 expression in QR.pa is asymmetric in the size and fate of its daughter cells: the large anterior daughter cell QR.pa remains a neuroblast, while the small posterior daughter cell QR.pp undergoes apoptosis and is rapidly engulfed and degraded by neighboring cells (Sulston and Horvitz, 1977). Although we found that the division per se is not necessary for the late phase of mig-1 expression, the size difference of the daughter cells and potentially unequal inheritance of the mig-1 activator could contribute to the temporal regulation of mig-1 expression.

To investigate the role of the asymmetric division of QR.p in mig-1 regulation, we first compared mig-1 expression in the two daughter cells using a null mutation in the essential cell death regulator ced-3 to prevent apoptosis of QR.pp (Ellis and Horvitz, 1986). As expected, loss of ced-3 did not affect the upregulation of mig-1 expression in QR.pa (Figure 4A). In contrast, the number of mig-1 transcripts that could be detected in the surviving QR.pp cells remained low. These results indicate that mig-1 is differently regulated in the two daughter cells. Next, we investigated if disrupting the asymmetry of the QR.p division affects the expression of mig-1 in QR.pa and QR.pp. A key regulator of asymmetric neuroblast division is the serine/threonine kinase PIG-1/MELK (Cordes et al., 2006), which controls the size difference of the QR.p daughters by inducing posterior displacement of the mitotic spindle (Chien et al., 2013) and may also contribute to the asymmetric segregation of cell fate determinants (Chien et al., 2013). As expected from these previous studies, loss of pig-1 resulted in a variable defect in cell size asymmetry, with a significant reduction in the mean cell size difference between QR.pa and QR.pp (Figure 4B). smFISH analysis showed that mig-1 expression in QR.pa was significantly decreased under these conditions (Figure 4A). However, there was no corresponding increase in the number of mig-1 transcripts in QR.pp, and no evidence of active transcription – which is visible as one or two bright smFISH spots in the nucleus (Ji et al., 2013) – was observed. To examine if mig-1 expression in QR.pa correlates with the level of cell size asymmetry, we compared the size ratio of QR.pa/QR.pp pairs with mig-1 expression in QR.pa. As shown in Figure 4C, we found a moderate correlation, with QR.pa cells from less asymmetric pairs showing lower mig-1 expression than QR.pa cells from more asymmetric pairs.

Figure 4 with 1 supplement see all

Download asset Open asset

QR.p division asymmetry influences *mig-1* expression dynamics in QR.pa.

(A) Single-molecule FISH (smFISH) quantification of *mig-1* expression in QR and its descendants relative to the position of the seam cells H2 to V5 in control, *ced-3(n717*), and *pig-1(gm344*) mutants. *mig-1* expression in QR.pa is not significantly different between control and *ced-3* mutants (p=0.36, n=34, Welch’s t-test) but is significantly reduced in *pig-1* mutants (p<0.0001, n=56). Data are from four independent experiments with three replicates each. (B) The size ratio of QR.pa and QR.pp in control (*heIs63*) and *pig-1(gm344*) mutants (***p<0.0001, n=48, from three independent experiments with three replicates each). Whiskers represent 1.5x interquartile range. (C) The QR.pa/QR.pp size ratio and *mig-1* expression show a moderately positive correlation in *pig-1(gm344*) mutants (Pearson correlation, R=0.47, n=48). (D) Nuclear size is not significantly different in QR.p and its daughter QR.pa in wild type control animals (p=0.22, n=12, Welch’s t-test). Whiskers represent 1.5x interquartile range. (E) Mathematical model incorporating cell division, fit to control data (see Materials and methods and Figure 4—figure supplement 1) with wild type QR.pa/QR.pp size ratio $ρ = 3.5$ (dark colors) and for lowered size ratios $ρ = [2.7,1.8,1]$ (light colors). Dashed line: time corresponding to moment in QR lineage migration that seam cell V2 is reached. Gray box: time window corresponding to gray region in F. (F) Plot of *mig-1* molecule number in QR.pa vs. QR.pa/QR.pp size ratio in model at times between $t = 3$ and $t = 3.3$ (gray region) and at the midpoint time $t = 3.15$ (black line).

To understand the effect of division asymmetry on mig-1 expression, we returned to our activator-based mathematical model. First, we incorporated QR.p division into the model at a time inferred from the data (see Materials and methods). Because no active transcription of mig-1 was observed in QR.pp experimentally, we set the mig-1 production rate to zero in QR.pp. Despite changes in the overall cell size, we experimentally observed that the nuclear size of QR.pa is not significantly different from that of its parent QR.p (Figure 4D). Since the activator is likely a transcription factor and therefore acts in the nucleus, the nuclear concentration of the activator in QR.pa should thus depend on its molecule number but not the cell size. Therefore, as before, we assumed that the production of mig-1 in QR.pa is dependent on the number of molecules, and not the concentration, of the activator. Finally, because the nuclear envelope breaks down before division, we assumed that at division, both the activator and preexisting mig-1 transcripts are distributed according to the relative sizes of QR.pa and QR.pp. Thus, at the division time, both the activator and mig-1 transcript numbers in QR.p are multiplied by a factor of $f$ in QR.pa and $(1 - f)$ in QR.pp, where $f = ρ / (1 + ρ)$ for a given QR.pa/QR.pp size ratio $ρ$ . The solution to this model with parameters fit to the control data (see Materials and methods and Figure 4—figure supplement 1) and the wild type mean $ρ = 3.5$ is shown in Figure 4E (dark colors). We observed that decreasing the QR.pa/QR.pp size ratio leads to reduced mig-1 levels in the QR.pa dynamics (Figure 4E, light colors), consistent with the experiments (Figure 4A). Correspondingly, plotting the mig-1 level as a function of the size ratio gave a positive dependency in the model (Figure 4F), which is in agreement with the positive correlation in the experiments (Figure 4C). We conclude that the asymmetry of the QR.p division, with the concomitant difference in cell size and the share of the mig-1 activator that is inherited by the two daughter cells, is required for upregulation of mig-1 expression in QR.pa.

Positive feedback regulation enhances timing precision of mig-1 expression

Mathematical modeling of the regulator-based timekeeping mechanism has shown that temporal precision is enhanced when the regulated timekeeper molecule amplifies its own expression through positive feedback (Gupta et al., 2020), but the in vivo relevance of this predicted strategy is not known.

mig-1 encodes a Frizzled Wnt receptor that triggers a BAR-1/β-catenin-dependent signaling pathway in QR.pa that induces the expression of target genes that inhibit migration (Mentink et al., 2014; Rella et al., 2021). Interestingly, mRNA sequencing showed that mig-1 itself is among the genes that are significantly upregulated when BAR-1 signaling is prematurely activated in QR.p using a constitutively active, N-terminally truncated version of BAR-1 (ΔN-BAR-1; Rella et al., 2021). To examine this further, we quantified mig-1 expression in bar-1(ga80) null mutants (Eisenmann et al., 1998) using smFISH. As shown in Figure 5—figure supplement 1, mig-1 expression in QR.pa was significantly reduced in the absence of BAR-1 signaling. These results support the notion that MIG-1 stimulates its own expression through BAR-1/β-catenin-dependent positive feedback.

Next, we investigated whether positive feedback regulation of mig-1 expression influences timing precision. Because BAR-1 signaling inhibits QR.pa migration (Mentink et al., 2014), the mean position of QR.pa is shifted anteriorly in the bar-1 null mutant and posteriorly in ΔN-BAR-1 expressing animals (Figure 5A). We therefore used the coefficient of variation (CV, standard deviation divided by the mean) in position as a relative measure of timing noise across these different conditions. As shown in Figure 5B, the CV in the position at which QR.pa expresses mig-1 at ≥25 transcripts was increased in both the bar-1 null mutant and animals expressing ΔN-BAR-1.

Figure 5 with 2 supplements see all

Download asset Open asset

BAR-1/β-catenin dependent autoregulation of *mig-1* expression contributes to timing precision.

(A) Single-molecule FISH (smFISH) quantification of *mig-1* expression in QR.p and QR.pa relative to the position of the seam cells H2 to V5 in control, *bar-1(ga80*) null mutants and animals expressing a constitutively active, N-terminally truncated version of BAR-1 (ΔN-BAR-1) in a *mab-5*(*gk670*) mutant background. Data are from four independent experiments with three replicates each. (B) Coefficient of variation (CV) squared of the position of QR.pa cells with ≥25 *mig-1* smFISH spots. Brown–Forsythe test for thresholds of 10–25 transcripts gives p=0.07–0.9 for control vs. *bar-1(ga80*), and 3×10^–6 to 0.008 for control vs. ΔN-BAR-1. Inset: CV² in the model in which positive feedback acts on *mig-1* both through the activator and independently of the activator (see Figure 5—figure supplement 2, purple).

To model the positive feedback, we considered three possibilities: the positive feedback (1) acts on mig-1 independently of the activator, (2) acts on mig-1 through the activator, or (3) both (Figure 5—figure supplement 2A). In all cases, we neglected cell division and mig-1 degradation for simplicity, and we calculated the CV in the time at which $m$ reaches a threshold numerically from the master equation following our previous work (Gupta et al., 2018). In case 1, we found that removing the feedback has little effect on the CV because the removal does not affect the dynamics of the activator (Figure 5—figure supplement 2), and the timing precision of an activated species is only marginally improved by autoregulation (Gupta et al., 2020). In case 2, because removing the feedback affects the dynamics of the activator but not the activation function of mig-1, we found that the removal produces mig-1 dynamics that are dramatically different from those observed in the bar-1 null and ΔN-BAR-1 data (Figure 5—figure supplement 2). In case 3, because removing the feedback affects both the dynamics of the activator and the activation function of mig-1, we found that the removal produces linear mig-1 dynamics whose CV is significantly larger than that with the feedback intact (Figure 5B inset and Figure 5—figure supplement 2). Since the bar-1 null and ΔN-BAR-1 data are more linear (Figure 5A) and have a higher CV (Figure 5B) than the control data, we conclude that case 3 best describes the positive feedback. These results provide experimental support for the model prediction that positive feedback of a regulated timekeeper enhances timing precision.

Discussion

Mathematical modeling of biological timekeeping mechanisms has provided insight into regulatory strategies that can achieve hightemporal precision (Gupta et al., 2020; Gupta et al., 2018), but the relevance of these predicted mechanisms for in vivo timekeeping remains to be determined. Using the time-dependent expression of the Wnt receptor gene mig-1/Frizzled in the C. elegans QR neuroblast lineage as a model system, we provide experimental evidence for a mechanism that combines regulation and positive feedback to mediate precise temporal control of gene expression.

In its most basic form, timing can be achieved through the linear increase of a stable timekeeper molecule that induces a response when a specific abundance threshold is reached (Ghusinga et al., 2017). Using mathematical modeling, we have previously shown that the temporal precision at which this threshold is crossed is enhanced when the concentration of the timekeeper molecule increases in a non-linear fashion due to upstream regulation (Gupta et al., 2018). Moreover, we showed that such a non-linear increase can be mediated through either the gradual accumulation of a transcriptional activator or the gradual decrease of a transcriptional repressor, which induces the expression of the timekeeper gene in a threshold-dependent manner (Gupta et al., 2018). As an experimental model to test these predictions, we used the time-dependent expression of the Wnt receptor mig-1 (Mentink et al., 2014; Rella et al., 2021). Quantification of mig-1 expression during QR lineage progression showed that it is sharply upregulated in QR.pa, which is consistent with a non-linear mode of regulation. However, we found that mig-1 is also expressed in the QR founder cell, raising the question whether these early and late phases of mig-1 expression are linked. Examination of evolutionarily conserved cis-regulatory elements that we deleted from the endogenous mig-1 locus using CRSPR/Cas9-mediated genome-editing showed that the two phases are independently regulated. Thus, we identified a small region in the first intron that is required for the early expression of mig-1 in QR but is dispensable for the late phase of expression in QR.pa. Conversely, two partially redundant regions upstream of the mig-1 coding sequence are required for the late phase of expression but have no detectable role in the early expression. The early phase of mig-1 expression could therefore be disregarded in the analysis of the late phase of expression, which was confirmed using a dynamic model that was fit to the mig-1 expression data.

The dynamic model also made predictions on the regulation of mig-1 expression. In an activator-based model, removal of the activator disrupts the upregulation of mig-1 expression in QR.pa, while in the case of a repressor, removal of the repressor results in premature expression of mig-1 in QR.p. We found that mig-1 expression in QR.pa was lost when the presumptive binding sites of the mig-1 regulator were deleted from the mig-1 upstream region, which supports the conclusion that the time-dependent expression of mig-1 depends on transcriptional activation and not on repression.

Another prediction from our mathematical modeling is that the timing precision of a regulated timekeeper system can be enhanced when the timekeeper molecule also stimulates its own expression through positive feedback (Gupta et al., 2020). mig-1 encodes a Frizzled-type Wnt receptor that triggers a canonical, BAR-1/β-catenin-dependent signaling pathway which induces the expression of target genes controlling QR.pa migration (Mentink et al., 2014; Rella et al., 2021). We found that the expression of mig-1 is enhanced through BAR-1/β-catenin signaling, creating a positive feedback loop. Using the relative position at which QR.pa upregulates mig-1 expression as a measure of timing precision, we discovered that timing was more variable when this feedback loop was disrupted through loss or constitutive activation of the BAR-1 pathway. Our modeling suggested that this effect is consistent with feedback on mig-1 that acts both (1) independently of the activator and (2) through the activator pathway; either type alone was insufficient. We, therefore, speculate that BAR-1/β-catenin-dependent positive feedback on mig-1 acts through multiple pathways. Our model assumed that disrupting the feedback disrupts the activator dynamics completely (the activator level becomes static in time). We also assumed that all of the positional variability in the disrupted cases comes from the expression dynamics and not, for example, from altered variability in the migration speed. Investigating feedback-independent contributions to the activator dynamics or perturbations to migration speed are interesting avenues for future study.

An intriguing aspect of the mig-1 timer is that it acts across the division of QR.p into QR.pa and its apoptotic sister cell QR.pp. Inhibition of QR.p mitosis and cell cycle progression through depletion of the M-phase regulator CDK-1 showed that the temporal regulation of mig-1 expression is independent of the division, ruling out a model in which the expression of mig-1 is linked to the invariant developmental timing of C. elegans cell divisions. However, even though the division itself is not necessary for the upregulation of mig-1 expression, we found that the asymmetry of the QR.p division influences the expression dynamics of mig-1. Our mathematical model allowed us to make predictions regarding the effect of cell division on the timing of mig-1 expression upregulation. First, the model assumed that upregulation is governed by the activator molecule number, not concentration. This assumption is justified by our finding that the activator acts transcriptionally in a cell whose nuclear size does not change upon division. However, for regulators that act post-transcriptionally or for cells whose nuclear volume changes upon division, the model would predict that the upregulation is governed by the cellular or nuclear concentration of the regulator, respectively, and not its molecule number. Second, the model assumes that the activator and mig-1 transcripts are partitioned at division according to cell volume. This assumption is critical, particularly for the activator. If instead the activator were sequestered to QR.pa, for example, the molecule number would not change after division, and the dynamics of mig-1 expression upregulation, therefore, would not depend on the QR.pa/QR.pp size ratio. Consequently, we predict that the activator is not sequestered at division and that unequal partitioning of the mig-1 activator between the differently sized daughter cells is necessary for the upregulation of mig-1 expression in QR.pa.

Examination of conserved cis-regulatory elements upstream of mig-1 showed that two small, partially redundantly acting regions are required for the upregulation of mig-1 expression in QR.pa. Chromatin immunoprecipitation followed by deep sequencing experiments (Araya et al., 2014) show that both regions overlap with binding peaks of multiple transcription factors, including homeobox transcription factors that are known to regulate Q neuroblast descendant migration (Clark et al., 1993), providing an entry point for the identification of the transcriptional regulators that control mig-1 expression. Molecular insight into the expression dynamics and activity of these regulators will provide further insight into how robust timing precision can be mediated through a regulator-based timekeeping mechanism.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Gene (Caenorhabditis elegans)	mig-1	Wormbase	Wormbase ID: WBGene00003238
Gene (C. elegans)	cdk-1	Wormbase	Wormbase ID: WBGene00000405
Gene (C. elegans)	ced-3	Wormbase	Wormbase ID: WBGene00000417
Gene (C. elegans)	pig-1	Wormbase	Wormbase ID: WBGene00021012
Gene (C. elegans)	bar-1	Wormbase	Wormbase ID: WBGene00000238
Strain and strain background (C. elegans)	For strains used in this study, see Supplementary file 1b	This paper		Strains are available upon request.
Recombinant DNA reagent	Pegl-17::tir1::TagBFP::unc-54 3’ UTR in pDESTR4-R3	This paper	pKN618	See Materials and methods. Plasmid is available upon request.
Sequence-based reagent	For sgRNA and ssODN sequences used for CRISPR/Cas9-based gene editing, see Supplementary file 1c	This paper
Chemical compound and drug	Auxin (natural indole-3-acetic acid)	Alfa Aesar	#A10556
Software and algorithm	R Studio	https://www.r-project.org/
Software and algorithm	MATLAB	https://www.mathworks.com/

C. elegans strains and culture

Request a detailed protocol

All C. elegans strains were grown at 20°C using standard culture conditions unless otherwise stated. The Bristol N2 strain was used as a wild type control. The alleles and transgenes generated and used in this work are: LGI: mig-1(hu295)[Δintron A]; mig-1(hu299)[Δintron B]; mig-1(hu314)[Δupstream 1]; mig-1(hu315)[Δupstream 2]; mig-1(hu335)[Δupstream 1–2]; LGIII: huIs210[Pegl-17::tir1::TagBFP, Pmyo-2::TdTomato], cdk-1(hu277[AID::cdk-1]); mab-5(gk670); LGIV: huIs179[Pegl-17::ΔN-bar-1, Pmyo-2::mCherry]; pig-1(gm344); ced-3(n717); LGV: heIs63[Pwrt-2::GFP::PH, Pwrt-2::GFP::H2B, Plin-48::mCherry]; LGX: bar-1(ga80); huIs166[Pwrt-2::mCherry::PH, Pwrt-2::mCherry::H2B, dpy-20(+)]. Details on CRISPR/Cas9-mediated genome edits can be found in Supplementary file 1a and a list of strains used in this study is in Supplementary file 1b.

C. elegans expression constructs and transgenesis

Request a detailed protocol

The Q lineage-specific TIR1 expressing transgene huEx744 was generated by microinjection of expression construct pKN618 (Pegl-17::tir1::TagBFP::unc-54 3’ untranslated region [UTR]) in the pDESTR4-R3 backbone plasmid at 20 ng/µl, using 5 ng/µl Pmyo-2::mCherry as co-injection marker. pBluescriptII was added for a final DNA concentration of 150 ng/µl in the microinjection mix. huIs210 was generated by genomic integration of huEx744 using gamma irradiation.

CRISPR/Cas9-mediated genome editing

Request a detailed protocol

The allele cdk-1(hu277[AID::cdk-1]) was generated using CRISPR/Cas9-based genome editing, dpy-10 co-conversion, and single-stranded DNA repair templates as previously described (Paix et al., 2015). Alleles mig-1(hu295) and mig-1(hu299) were generated using CRISPR/Cas9-based genome editing, pha-1 co-conversion, and single-stranded DNA repair templates (Ward, 2015). Alleles mig-1(hu314), mig-1(hu315), and mig-1(hu335) were generated using CRISPR/Cas9-based genome editing as described (Dokshin et al., 2018). PCR fragments containing the T7 promoter and the single guide RNA (sgRNA) sequence of interest were used to synthesize sgRNA in vitro. Repair templates and sgRNAs were co-injected with recombinant SpCas9 (D’Astolfo et al., 2015). sgRNAs and repair template sequences can be found in Supplementary file 1c.

Analysis of QR descendant migration

Request a detailed protocol

The final position of the QR descendant QR.pap was determined using widefield fluorescence microscopy at the late L1 larval stage (between 12 and 16 hr after hatching). In all experiments, transgenic animals expressing GFP in the seam cells and the Q lineage were used (transgene heIs63; Wildwater et al., 2011). The position of QR.pap was determined by measuring the distance of QR.pap to the seam cells V1.p and V2.p. This distance was then normalized to the distance between V1.p and V2.p to determine the relative position of QR.pap to these cells. The average of two position measurements (one using V1.p and one using V2.p) was used in the analyses. Measurements were performed using the FIJI package (Schindelin et al., 2012).

Single-molecule FISH

Request a detailed protocol

The smFISH protocol was performed as described (Ji et al., 2013; Ji and van Oudenaarden, 2012). In all experiments, transgenic animals expressing GFP or mCherry in the seam cells and the Q neuroblast lineage were used (heIs63 and huIs166, respectively; Wildwater et al., 2011). Synchronized L1 larvae were fixed using 4% paraformaldehyde and kept in 70% ethanol. Hybridization was performed overnight at 37°C, in the dark. Oligonucleotide probes used in hybridization were designed using the Stellaris probe designer and chemically coupled to fluorescent dyes Cy5 (mig-1; Ji and van Oudenaarden, 2012) or ATTO565 (larp-1, a positive control for smFISH staining quality). Animals were incubated in buffer containing DAPI for nuclear counterstaining prior to mounting. Z-stacks of 0.2–0.3 µm thickness were collected on a Perkin Elmer Ultraview VoX mounted on a Leica DMI6000 microscope with a 100×/NA 1.47 oil objective, or a Nikon Ti2 microscope with a 100×/NA 1.45 oil objective. Both systems used a Yokogawa X1 spinning disk and Hamamatsu Orca Flash4.0 V3 camera. Images were acquired at 1024×1024 resolution. mRNA spots were quantified manually from the z-stacks. Only spots detected in two or more consecutive focal planes were counted. Image analysis was performed in FIJI (Schindelin et al., 2012).

AID of CDK-1

Request a detailed protocol

The Arabidopsis thaliana TIR1 protein was expressed as an integrated transgene driven by the egl-17 promoter (huIs210; Branda and Stern, 2000). This ensured CDK-1 was only depleted in the Q neuroblast lineage. By transferring L1 larvae to plates containing 0.5 mM natural auxin indole-3-acetic acid (Alfa Aesar, #A10556) during the time window between QR and QR.p division (300–345 min after hatching), QR.p division was specifically blocked.

Quantification of DNA content

Request a detailed protocol

DNA content was determined using images from smFISH experiments as follows. Per cell, the range of focal planes spanning the nucleus was determined. A z-stack of the focal planes with the nucleus was then made using the ‘sum slices’ option in FIJI. A region of interest (ROI) was subsequently drawn around the nucleus, and total intensity was measured in the ROI. Intensity measurements from multiple seam (Vn.p) cells were averaged as the 2C internal control for each animal. The intensity of DAPI staining in QR.p was then divided by this average internal control to produce the DNA content ratio. A ratio of 1 corresponds to 2C, a ratio of 2–4C, and a ratio of 4–8C.

Identification of evolutionarily conserved regions in the mig-1 upstream region and first intron

Request a detailed protocol

The homologous gene of C. elegans mig-1 was first identified in 22 other Caenorhabditis species using the BLASTP algorithm on the following genomes: Caenorhabditis monodelphis (JU1667_v1), Caenorhabditis plicata (SB355_v1), Caenorhabditis virilis (JU1968_v1), Caenorhabditis quiockensis (C. sp. 38, JU2809_v1), Caenorhabditis castelli (JU1956_v1), Caenorhabditis uteleia (C. sp. 31, JU2585_v1), Caenorhabditis afra (JU1286_v1), Caenorhabditis sulstoni (C. sp.32, JU2788_v1), Caenorhabditis waitukubuli (C. sp. 39, NIC564_v1), Caenorhabditis panamensis (C. sp. 28, QG2080_v1), Caenorhabditis becei (C. sp. 29, QG2083_v1), Caenorhabditis nouraguensis (JU2079_v1), Caenorhabditis kamaaina (QG2077_v1), Caenorhabditis inopinata (C. sp. 34, NK74SC_v1), Caenorhabditis tropicalis (JU1373_301), Caenorhabditis doughertyi (JU1771_v1), Caenorhabditis brenneri (C_brenneri-6.0.1b), Caenorhabditis latens (PX534_v1), Caenorhabditis briggsae (CB4), Caenorhabditis sinica (JU800_v1), Caenorhabditis tribulationis (C. sp. 40, JU2818_v1), and Caenorhabditis zanzibari (C. sp. 26, JU2190_v1). BLASTP was performed on this set of species on http://caenorhabditis.org/. We then looked for motif conservation in mig-1 regulatory regions in this set of species using the MEME Suite tools (Bailey et al., 2009) with the following parameter settings: -mod zoops, -nmotifs 4, -minw 8, and -maxw 18 for the first intron and -mod zoops, -nmotifs 4, -minw 10, and -maxw 25 for the upstream region.

Fitting of mathematical model to data

Request a detailed protocol

We find the parameters of either the increasing activator model or the decreasing repressor model (Figure 2D) using a least-squares fit to the control data (Figure 2B). Specifically, to find $m (0)$ and $ν$ , we perform a linear fit to the log of the mig-1 mRNA number for the QR data only, corresponding to exponential decay; we find $m (0) = 31$ mRNAs and $ν^{- 1} = 0.65$ (in the arbitrary time units of Figure 2D). Then, to find the remaining parameters, we fit the numerical solution of the $d m / d t$ equation to all of the control data, minimizing the sum of the squares of both the vertical and horizontal residuals, scaled by the average value of mRNA number and time across all data, respectively. For the activator model, we find $α = 390$ , $K / k = 3.0$ , and $H = 14$ . For the repressor model, we find $α = 130$ , $μ H = 6.7$ , and $H l n (r_{0} / K) = 18$ . The high value of the Hill coefficients here and below is a reflection of the fact that mig-1 expression remains low in QR.p cells and then increases very sharply in QR.pa cells. Because we model this increase with a single regulating species, whereas mig-1 is likely regulated by more than one species in reality, the Hill coefficient can be viewed as an effective parameter rather than a biochemical constant. Furthermore, we see that some parameters need only be determined in particular combinations. The reason is the following. For the activator model, plugging $a (t) = k t$ into $F (t) = α a^{H} / (a^{H} + K^{H})$ obtains $F (t) = α / [1 + {(K / k t)}^{H}]$ , which depends only on the ratio $K / k$ . For the repressor model, plugging $r (t) = r_{0} e^{- μ t}$ into $F (t) = α K^{H} / (r^{H} + K^{H})$ obtains $F (t) = α / [1 + {(r_{0} e^{- μ t} / K)}^{H}] = α / {1 + e x p [H l n (r_{0} / K) - μ H t]}$ , which depends only on the combinations $H l n (r_{0} / K)$ and $μ H$ . Although these parameter combinations are sufficient to specify the mig-1 dynamics, for illustrative purposes, we plot example activator and repressor dynamics consistent with these values in Figure 2D by choosing $K = 20$ for the activator model (giving $k = 6.6$ ), and $K = 5$ and $H = 14$ for the repressor model (giving $r_{0} = 19$ and $μ = 0.48$ ). For the zero initial mig-1 condition (Figure 2D), we set $m (0) = 0$ .

To find the division time from the control data (Figure 4A), we count the cumulative number of QR.p data points as a function of position (moving posteriorly), count the cumulative number of QR.pa data points as a function of position (moving anteriorly), and find the position at which these numbers are equal. This position is $T = 2.32$ in the arbitrary time units of Figure 4E. The QR.pa dynamics from the model with division, using the wild type mean size ratio $ρ = 3.5$ , are fit to the control data using the same procedure as above, yielding the parameters $α = 240$ , $K / k = 2.6$ , and $H = 17$ . This fit is shown for QR.p, QR.pa, and QR.pp in Figure 4E (dark colors). To illustrate the effect of decreasing size ratio $ρ$ , we also plot these dynamics with $ρ = [2.7,1.8,1]$ in Figure 4E (light colors). To illustrate the dependence of the QR.pa mig-1 molecule number $m (t)$ on the size ratio $ρ$ , we plot $m (t)$ vs. $ρ$ at times between $t = 3$ and $t = 3.3$ (gray region) and at the midpoint time $t = 3.15$ (black line) in Figure 4F.

Modeling the effect of feedback on timing precision

Request a detailed protocol

To model the effect of BAR-1/ $β$ -catenin dependent feedback on the noise in mig-1 expression, we consider three possibilities: (1) the feedback acts on mig-1 independently of the activator (Figure 5—figure supplement 2A, left), (2) the feedback acts on mig-1 through the activator (Figure 5—figure supplement 2A, middle), or (3) both (Figure 5—figure supplement 2A, right). In all cases, we neglect cell division and mig-1 degradation for simplicity, and we set $t = 0$ using the rightmost position data point in Figure 5A (ΔN-BAR-1). We calculate the mean $\bar{t}$ and variance $σ_{t}^{2}$ in the time at which $m$ reaches its threshold $m_{*}$ numerically from the master equation following our previous work (Gupta et al., 2018).

Possibility 1 (Figure 5—figure supplement 2A, left) corresponds to the dynamics $a (t) = k t$ and $d m / d t = F (a, m)$ . Both the bar-1(ga80) null mutant and the constitutively active ΔN-BAR-1 break the feedback, corresponding to $F (a, m) = F (a)$ . In these cases, we take $F (a) = α / [1 + {(K / a)}^{H}]$ as before. For the control case, previous work Gupta et al., 2020 found that this topology reduces noise when the feedback is negative for low $a$ and positive for high $a$ . The simplest modification that achieves this property is $F (a, m) = α / [1 + {(K / a)}^{H_{0} m}]$ , corresponding to a regulation function that sharpens with $m$ . Because breaking the feedback can also affect the activation parameters for this topology, we fit the bar-1(ga80) and ΔN-BAR-1 cases separately from the control case. Fitting as above results in $α = 260$ , $K / k = 3.4$ , and $H_{0} = 6.7$ for the control (Figure 5—figure supplement 2B, left, red); $α = 490$ , $K / k = 4.2$ , and $H = 13$ for bar-1(ga80) (Figure 5—figure supplement 2B, middle, red); and $α = 490$ , $K / k = 4.2$ , and $H = 6$ for ΔN-BAR-1 (Figure 5—figure supplement 2B, right, red). Calculation of $\bar{t}$ and $σ_{t}^{2}$ requires specifying $k$ , and we find that for no $k$ value is the $C V = σ_{t} / \bar{t}$ significantly increased upon breaking the feedback (Figure 5—figure supplement 2C, middle); Figure 5—figure supplement 2 shows the results for $k = 30$ and $m_{*} = 25$ .

Possibilities 2 (Figure 5—figure supplement 2A, middle) and 3 (Figure 5—figure supplement 2A, right) both correspond to the dynamics $d a / d t = G (m)$ , but possibility 2 has $d m / d t = F (a)$ , whereas possibility 3 has $d m / d t = F (a, m)$ . Because we found above that breaking autoregulatory feedback has little effect on the $C V$ , we simplify the dynamics in possibility 3 to also read $d m / d t = F (a)$ . Nevertheless, we retain the key topological difference that breaking the feedback cannot change the activation parameters in possibility 2 but can in possibility 3. Thus, in possibility 2, the bar-1(ga80) and ΔN-BAR-1 cases inherit the activation parameters from the control case, whereas in possibility 3, the bar-1(ga80) and ΔN-BAR-1 cases are fit separately from the control case.

For the control case in possibilities 2 and 3, we take $F (a) = α_{0} + α / [1 + {(K / a)}^{H}]$ and analogously $G (m) = β_{0} + β / [1 + {(J / m)}^{H}]$ , where the Hill coefficients are equal for simplicity, and at least one of $α_{0}$ and $β_{0}$ must be nonzero to initiate production from zero molecules. Previous work Gupta et al., 2020 found that this topology reduces noise when the times at which $a$ and $m$ cross their thresholds $K$ and $J$ , respectively, are well separated. We focus on the regime where $a$ crosses first ( $K ≪ J$ ), although we will see below that noise is reduced even for $K > J$ . To reduce the number of fitting parameters, we recognize that the limit $K ≪ J$ corresponds to the simplification $G (m) = β_{0}$ . Therefore we fit in this limit, with $α_{0} = 0$ , and we later check the robustness of our results as we reintroduce $J$ and $β$ . Fitting results in $α = 2800$ , $K / β_{0} = 4.5$ , and $H = 12$ (Figure 5—figure supplement 2B, left, cyan and purple).

For the bar-1(ga80) and ΔN-BAR-1 cases in possibilities 2 and 3, we recognize that, unlike in possibility 1, breaking the feedback disrupts the production of the activator. In these cases, we take $a$ to be a constant $a_{0}$ , and therefore, we have $d m / d t = F (a_{0})$ . As discussed above, in possibility 2, the parameters of $F$ are inherited from the control case. Consequently, the bar-1(ga80) case, which corresponds to a low value of $a_{0}$ , produces a shallow linear increase in mig-1 (Figure 5—figure supplement 2B, middle, cyan); while the ΔN-BAR-1 case, which corresponds to a high value of $a_{0}$ , produces a steep linear increase in mig-1 (Figure 5—figure supplement 2B, right, cyan). Because neither is a reasonable description of the data, we conclude that possibility 2 is not a good model of the feedback, and we do not consider the $C V$ for possibility 2.

In possibility 3, as discussed above, the bar-1(ga80) and ΔN-BAR-1 cases are fit separately from the control case. Fitting results in $F (a_{0}) = 5.1$ for bar-1(ga80) (Figure 5—figure supplement 2B, middle, purple) and $F (a_{0}) = 7.7$ for ΔN-BAR-1 (Figure 5—figure supplement 2B, right, purple). The $C V$ in these cases can be calculated analytically (Gupta et al., 2018) as $C V = m_{*}^{- 1 / 2}$ . The control case requires specifying $β_{0}$ , and we find that for $β_{0} ≳ 10$ , the $C V$ is increased upon breaking the feedback (Figure 5—figure supplement 2C, right) in agreement with the experimental data (Figure 5B and Figure 5—figure supplement 2C, left). Indeed, we find that this property remains true even for $J$ as small as 1 and for $β ≫ β_{0}$ . Figure 5B (inset) and Figure 5—figure supplement 2 show the results for $β_{0} = 25$ , $J = 10$ , $β = 75$ , and $m_{} = 25$ .

Statistical ANOVA data

Request a detailed protocol

To compare the CVs of mig-1 expression in bar-1 mutants (Figure 5B), we performed a Brown–Forsythe ANOVA test. Specifically, we made two comparisons (control vs. bar-1 null and control vs. ΔN-BAR-1) and, therefore, multiplied p-values by two to account for the Bonferroni correction. Position data were included if the mig-1 mRNA molecule number was greater than or equal to a threshold and rescaled by their mean to obtain the CV. We performed the test for thresholds in the range 10–25. For the control vs. bar-1 null comparison, the Bonferroni-corrected p-values ranged from 0.07 to 0.9, with a mean value of 0.4. For the control vs. ΔN-BAR-1comparison, the Bonferroni-corrected p-values ranged from 3×10^–6 to 0.008, with a mean value of 0.001.

Statistical testing and modeling software

Request a detailed protocol

Statistical testing and plot generation was performed using R 4.1.0 and RStudio 1.4.1717, using various packages from the tidyverse. Mathematical modeling was performed in Matlab.

Data availability statement

Request a detailed protocol

Source datafiles containing the numerical data used to generate the figures can be accessed at https://github.com/erikschild/mig1_timer_code, (copy archived at Schild, 2023a) and the code used for modeling can be found at https://github.com/amugler/mig1, (copy archived at Schild, 2023b).

Data availability

All data generated or analysed during this study are included in the manuscript and figures. Source datafiles containing the numerical data used to generate the figures can be accessed at https://github.com/erikschild/mig1_timer_code, (copy archived at Schild, 2023a).

References

1. Angers S
2. Moon RT
(2009) Proximal events in WNT signal transduction
Nature Reviews. Molecular Cell Biology 10:468–477.
https://doi.org/10.1038/nrm2717
- PubMed
- Google Scholar
1. Araya CL
2. Kawli T
3. Kundaje A
4. Jiang L
5. Wu B
6. Vafeados D
7. Terrell R
8. Weissdepp P
9. Gevirtzman L
10. Mace D
11. Niu W
12. Boyle AP
13. Xie D
14. Ma L
15. Murray JI
16. Reinke V
17. Waterston RH
18. Snyder M
(2014) Regulatory analysis of the C. elegans genome with Spatiotemporal resolution
Nature 512:400–405.
https://doi.org/10.1038/nature13497
- PubMed
- Google Scholar
1. Aydogan MG
2. Steinacker TL
3. Mofatteh M
4. Wilmott ZM
5. Zhou FY
6. Gartenmann L
7. Wainman A
8. Saurya S
9. Novak ZA
10. Wong S-S
11. Goriely A
12. Boemo MA
13. Raff JW
(2020) An autonomous oscillation times and executes Centriole Biogenesis
Cell 181:1566–1581.
https://doi.org/10.1016/j.cell.2020.05.018
- PubMed
- Google Scholar
1. Bailey TL
2. Boden M
3. Buske FA
4. Frith M
5. Grant CE
6. Clementi L
7. Ren J
8. Li WW
9. Noble WS
(2009) MEME SUITE: tools for motif discovery and searching
Nucleic Acids Research 37:W202–W208.
https://doi.org/10.1093/nar/gkp335
- Google Scholar
1. Baudet ML
2. Zivraj KH
3. Abreu-Goodger C
4. Muldal A
5. Armisen J
6. Blenkiron C
7. Goldstein LD
8. Miska EA
9. Holt CE
(2011) miR-124 acts through Corest to control onset of Sema3A sensitivity in navigating retinal growth cones
Nature Neuroscience 15:29–38.
https://doi.org/10.1038/nn.2979
- PubMed
- Google Scholar
(2005) Circadian rhythms from multiple Oscillators: lessons from diverse organisms
Nature Reviews. Genetics 6:544–556.
https://doi.org/10.1038/nrg1633
- PubMed
- Google Scholar
1. Branda CS
2. Stern MJ
(2000) Mechanisms controlling sex myoblast migration in Caenorhabditis elegans hermaphrodites
Developmental Biology 226:137–151.
https://doi.org/10.1006/dbio.2000.9853
- PubMed
- Google Scholar
1. Campbell DS
2. Regan AG
3. Lopez JS
4. Tannahill D
5. Harris WA
6. Holt CE
(2001) Semaphorin 3A elicits stage-dependent collapse, turning, and branching in Xenopus retinal growth cones
The Journal of Neuroscience 21:8538–8547.
https://doi.org/10.1523/JNEUROSCI.21-21-08538.2001
- PubMed
- Google Scholar
(2013) Caenorhabditis elegans PIG-1/MELK acts in a conserved PAR-4/Lkb1 polarity pathway to promote asymmetric Neuroblast divisions
Genetics 193:897–909.
https://doi.org/10.1534/genetics.112.148106
- PubMed
- Google Scholar
(1993) Control of cell fates in the central body region of C. elegans by the Homeobox Gene Lin-39
Cell 74:43–55.
https://doi.org/10.1016/0092-8674(93)90293-y
- PubMed
- Google Scholar
1. Clevers H
2. Nusse R
(2012) Wnt/Β-Catenin signaling and disease
Cell 149:1192–1205.
https://doi.org/10.1016/j.cell.2012.05.012
- PubMed
- Google Scholar
(2006) The C. elegans melk ortholog pig-1 regulates cell size asymmetry and daughter cell fate in asymmetric neuroblast divisions
Development 133:2747–2756.
https://doi.org/10.1242/dev.02447
- PubMed
- Google Scholar
1. Coudreuse D
2. Nurse P
(2010) Driving the cell cycle with a minimal CDK control network
Nature 468:1074–1079.
https://doi.org/10.1038/nature09543
- PubMed
- Google Scholar
1. D’Astolfo DS
2. Pagliero RJ
3. Pras A
4. Karthaus WR
5. Clevers H
6. Prasad V
7. Lebbink RJ
8. Rehmann H
9. Geijsen N
(2015) Efficient intracellular delivery of native proteins
Cell 161:674–690.
https://doi.org/10.1016/j.cell.2015.03.028
- PubMed
- Google Scholar
1. Dequéant M-L
2. Pourquié O
(2008) Segmental Patterning of the vertebrate embryonic axis
Nature Reviews. Genetics 9:370–382.
https://doi.org/10.1038/nrg2320
- PubMed
- Google Scholar
(2018) Robust genome editing with short single-stranded and long, partially single-stranded DNA donors in Caenorhabditis elegans
Genetics 210:781–787.
https://doi.org/10.1534/genetics.118.301532
- PubMed
- Google Scholar
1. Dubois C
2. Gupta S
3. Mugler A
4. Félix M-A
(2021) Temporally regulated cell migration is sensitive to variation in body size
Development 148:dev196949.
https://doi.org/10.1242/dev.196949
- PubMed
- Google Scholar
1. Eisenmann DM
2. Maloof JN
3. Simske JS
4. Kenyon C
5. Kim SK
(1998) The Β-Catenin Homolog BAR-1 and LET-60 Ras coordinately regulate the hox gene Lin-39 during Caenorhabditis elegans vulval development
Development 125:3667–3680.
https://doi.org/10.1242/dev.125.18.3667
- PubMed
- Google Scholar
1. Ellis HM
2. Horvitz HR
(1986) Genetic control of programmed cell death in the Nematode C. elegans
Cell 44:817–829.
https://doi.org/10.1016/0092-8674(86)90004-8
- PubMed
- Google Scholar
(2017) First-passage time approach to controlling noise in the timing of intracellular events
PNAS 114:693–698.
https://doi.org/10.1073/pnas.1609012114
- PubMed
- Google Scholar
1. Gliech CR
2. Holland AJ
(2020) Keeping track of time: the fundamentals of cellular clocks
The Journal of Cell Biology 219:e202005136.
https://doi.org/10.1083/jcb.202005136
- PubMed
- Google Scholar
(2010) Finding regulatory DNA motifs using alignment-free evolutionary conservation information
Nucleic Acids Research 38:e90.
https://doi.org/10.1093/nar/gkp1166
- PubMed
- Google Scholar
(2018) Temporal precision of regulated gene expression
PLOS Computational Biology 14:e1006201.
https://doi.org/10.1371/journal.pcbi.1006201
- PubMed
- Google Scholar
(2020) Temporal precision of molecular events with regulation and feedback
Physical Review. E 101:062420.
https://doi.org/10.1103/PhysRevE.101.062420
- PubMed
- Google Scholar
1. Harbison CT
2. Gordon DB
3. Lee TI
4. Rinaldi NJ
5. Macisaac KD
6. Danford TW
7. Hannett NM
8. Tagne J-B
9. Reynolds DB
10. Yoo J
11. Jennings EG
12. Zeitlinger J
13. Pokholok DK
14. Kellis M
15. Rolfe PA
16. Takusagawa KT
17. Lander ES
18. Gifford DK
19. Fraenkel E
20. Young RA
(2004) Transcriptional regulatory code of a Eukaryotic genome
Nature 431:99–104.
https://doi.org/10.1038/nature02800
- PubMed
- Google Scholar
1. Hedgecock EM
2. White JG
(1985) Polyploid tissues in the nematode Caenorhabditis elegans
Developmental Biology 107:128–133.
https://doi.org/10.1016/0012-1606(85)90381-1
- PubMed
- Google Scholar
1. Ji N
2. van Oudenaarden A
(2012) Single molecule fluorescent in situ hybridization (smFISH) of C. elegans worms and embryos
WormBook 1:1–16.
https://doi.org/10.1895/wormbook.1.153.1
- PubMed
- Google Scholar
(2013) Feedback control of gene expression variability in the Caenorhabditis elegans WNT pathway
Cell 155:869–880.
https://doi.org/10.1016/j.cell.2013.09.060
- PubMed
- Google Scholar
1. Mentink RA
2. Middelkoop TC
3. Rella L
4. Ji N
5. Tang CY
6. Betist MC
7. van Oudenaarden A
8. Korswagen HC
(2014) Cell intrinsic modulation of WNT signaling controls Neuroblast migration in C. elegans
Developmental Cell 31:188–201.
https://doi.org/10.1016/j.devcel.2014.08.008
- PubMed
- Google Scholar
1. Nitta KR
2. Jolma A
3. Yin Y
4. Morgunova E
5. Kivioja T
6. Akhtar J
7. Hens K
8. Toivonen J
9. Deplancke B
10. Furlong EEM
11. Taipale J
(2015) Conservation of transcription factor binding Specificities across 600 million years of Bilateria evolution
eLife 4:e04837.
https://doi.org/10.7554/eLife.04837
- PubMed
- Google Scholar
1. Oates AC
(2020) Waiting on the fringe: cell autonomy and signaling delays in Segmentation clocks
Current Opinion in Genetics & Development 63:61–70.
https://doi.org/10.1016/j.gde.2020.04.008
- PubMed
- Google Scholar
(2015) High efficiency, Homology-directed genome editing in Caenorhabditis elegans using CRISPR-Cas9 Ribonucleoprotein complexes
Genetics 201:47–54.
https://doi.org/10.1534/genetics.115.179382
- PubMed
- Google Scholar
1. Pickering J
2. Rich CA
3. Stainton H
4. Aceituno C
5. Chinnaiya K
6. Saiz-Lopez P
7. Ros MA
8. Towers M
(2018) An intrinsic cell cycle Timer TERMINATES limb Bud outgrowth
eLife 7:e37429.
https://doi.org/10.7554/eLife.37429
- PubMed
- Google Scholar
1. Ray A
2. Zhu H
3. Ding A
4. Li X
(2022) Transcriptional and epigenetic regulation of temporal Patterning in neural Progenitors
Developmental Biology 481:116–128.
https://doi.org/10.1016/j.ydbio.2021.10.006
- PubMed
- Google Scholar
(2021) A switch from Noncanonical to Canonical WNT signaling stops Neuroblast migration through a SLT-Robo and RGA-9B/ARHGAP-dependent mechanism
PNAS 118:e2013239118.
https://doi.org/10.1073/pnas.2013239118
- PubMed
- Google Scholar
Software
1. Schild ES
(2023a) Mig1_Timer_Code, version swh:1:rev:fa01c90606c4175f06a3b7f61d4966a8a7800301
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:48ba37ce42a15d029f6d37eb3a21271a4f7d1b53;origin=https://github.com/erikschild/mig1_timer_code;visit=swh:1:snp:6b3c4e6c3be4eab9e34905c92784bbb76a764c73;anchor=swh:1:rev:fa01c90606c4175f06a3b7f61d4966a8a7800301
Software
1. Schild ES
(2023b) Mig1, version swh:1:rev:cfae1d435109365b2bd7bb8866ac4a9de1c9f0ef
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:54587a2f1782b09c6d3e24a558ac0d7198c64a58;origin=https://github.com/amugler/mig1;visit=swh:1:snp:8183239be2fe60d9ebb52f1faafcb939b1977758;anchor=swh:1:rev:cfae1d435109365b2bd7bb8866ac4a9de1c9f0ef
1. Schindelin J
2. Arganda-Carreras I
3. Frise E
4. Kaynig V
5. Longair M
6. Pietzsch T
7. Preibisch S
8. Rueden C
9. Saalfeld S
10. Schmid B
11. Tinevez J-Y
12. White DJ
13. Hartenstein V
14. Eliceiri K
15. Tomancak P
16. Cardona A
(2012) Fiji: an open-source platform for biological-image analysis
Nature Methods 9:676–682.
https://doi.org/10.1038/nmeth.2019
- PubMed
- Google Scholar
1. Sulston JE
2. Horvitz HR
(1977) Post-embryonic cell lineages of the Nematode, Caenorhabditis elegans
Developmental Biology 56:110–156.
https://doi.org/10.1016/0012-1606(77)90158-0
- PubMed
- Google Scholar
1. Tyson JJ
2. Novak B
(2008) Temporal organization of the cell cycle
Current Biology 18:R759–R768.
https://doi.org/10.1016/j.cub.2008.07.001
- PubMed
- Google Scholar
1. Ward JD
(2015) Rapid and precise engineering of the Caenorhabditis elegans genome with lethal Mutation Co-conversion and inactivation of NHEJ repair
Genetics 199:363–377.
https://doi.org/10.1534/genetics.114.172361
- PubMed
- Google Scholar
(2011) Cell shape and WNT signaling Redundantly control the division axis of C. elegans epithelial stem cells
Development 138:4375–4385.
https://doi.org/10.1242/dev.066431
- PubMed
- Google Scholar
1. Zhang L
2. Ward JD
3. Cheng Z
4. Dernburg AF
(2015) The Auxin-inducible degradation (AID) system enables Versatile conditional protein depletion in C. elegans
Development 142:4374–4384.
https://doi.org/10.1242/dev.129635
- PubMed
- Google Scholar

Decision letter

Marianne E Bronner

Senior and Reviewing Editor; California Institute of Technology, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Precise temporal control of neuroblast migration through combined regulation and feedback of a Wnt receptor" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Marianne Bronner as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. Both reviewers raise the lack of integration between the model and data.

2. Secondly, the evidence for the positive feedback loop is weak/lacking. Additional experiments would be optimal. Alternatively, modelling could be used to interpret or strengthen these results, but this would require explicitly including the positive feedback loop and stochastic variability. At least it would be important to substantially rephrase the conclusions, to recognize that the evidence is relatively weak. I refer you to the comments of the reviewers for details.

Reviewer #1 (Recommendations for the authors):

- l. 151. "2C". Given the article's mixed audience of biologists and modelers, it would be good to explain the 2C, 4C, and 8C notation. Currently, it is not easy to deduce the meaning of 2C for the article text, if you do not already know it.

- l. 154. It is not explained why it is expected that CDK1-inhibited neuroblasts are 4C, not 2C. This is presumably because inhibiting CDK1 does not inhibit DNA duplication but does inhibit the subsequent division. It would be good to mention this explicitly.

- l. 189-208. In general, the fitting procedures outlined in this section and the Materials and methods seem fine, but it would be helpful if they would be presented in a way that allows visual comparison of the quality of the fit. Currently, smFISH data is always shown separately from fits, and plots with fitting data are mostly used to highlight predictions of the model. It would be good if somewhere, either in the main text or in the SI, the fitting data is plotted on top of the data, just to see how good the fit actually is. I realize that data and fit have different time axes, but this should not be a fundamental problem in plotting data and fit together if I understand the fitting procedure correctly.

- l. 189-208, l. 508-533. The fitted parameter values are only mentioned in the Methods. The values for the Hill coefficients for both the repressor and activator model are very high, H=14 or H=17, but this is not really commented on. I understand from the rest of the paper that mig-1 likely has some form of positive feedback on its own expression, which could in principle explain high Hill coefficients. But are such high Hill coefficients consistent with such a positive feedback mechanism? I think the authors should discuss the high Hill coefficient and whether such high values are plausible based on their proposed mechanism.

- l. 228, Figure 3A, C. The model predicts that decreasing the activator (hereby mutating activator binding sites) leads to a slower increase in mig-1 level. In the experiments, this is not visible for Δ upstream 1, but here a substantial number of animals appear to fail in inducing mig-1. For Δ upstream 2, there are even more animals like this. This fact is mentioned in the main text, but not interpreted in any way. How do the authors interpret this in light of the activator model?

- l. 245. 'required for the rapid upregulation'. To me, this formulation seems too strong. Looking at Figure 4A, I still see rapid upregulation in pig-1 mutants, as in happening quickly after the division of QR.p, but with lower mig-1 levels reached compared to the control.

- l. 245-298, Figure 4. Related to the above point, it is not clear to me how the authors interpret the pig-1 experiments in relation to the fitting results in Figure 4E. The model predicts that decreasing the size ratio between QR.pa and QR.pp leads to a slower accumulation of mig-1, but with the same mig-1 level ultimately reached, albeit at a later time. The pig-1 data in Figure 4A could be consistent with slower accumulation, but, in contrast to Figure 4E, mig-1 expression never reaches the peak control level. How do the authors interpret this contrasting result? An alternative model to explain the data could be that in pig-1 mutants mig-1 expression is induced at the same time as in control animals, but with a lower 'amplitude' of mig-1 expression. Can the authors rule out such an alternative model? Finally, in the author's model QR.pa migration stops once a specific mig-1 level is reached. In that case, wouldn't one expect a positional defect for pig-1 mutants? I didn't see this mentioned anywhere in the article.

- l. 281-285. I found the explanation here confusing. I think I understand the main point: the cells may have different sizes, but have the same nuclear volume. As the activator is likely in the nucleus, one can therefore just use the same molecular level as before, without taking into account the difference in volume. However, the way it is written here suggested initially to me that a modification was made to the model by focusing on a number of molecules rather than concentration, even though, if I understand correctly, the model was formulated in terms of a number of molecules from the beginning. It would be good to clarify this.

- l. 296-298. I have the same issue here as above: it is not clear to me that division asymmetry is required for rapid upregulation. The proposed role of asymmetry can also be formulated more sharply. If I understand it correctly, according to the model, the only impact of size asymmetry is in the asymmetric inheritance of the activator. In particular, the subsequent rate of production of the activator in Q.pa does not depend on its size, correct? Or am I missing something?

- l. 300-329. I find this section the weakest of the manuscript, both in terms of experiments and modelling. The main aim of this section is to test the model prediction from an earlier paper (Gupta 2020) that positive feedback of mig-1 on its own expression reduces variability in timing. The main experimental weakness is that the authors do not specifically perturb this positive feedback loop, but instead use two Wnt mutants (the loss-of-function bar-1(0) and constitutively active DeltaN-BAR-1) that likely impact many genes, including those responsible for migration (l. 95,96), as also evidenced by the resulting change in average Q.pa position. While the authors find an increase in variability of Q.pa position, this could also be explained by changes in the expression of any of the other BAR-1 targets. The proper experiment, I think, would be to remove the POP-1/TCF sites in the mig-1 promoter, although there might be experimental reasons why this is not feasible. The last sentence of the section (l. 328-329) is a bit careful in drawing strong conclusions from these result, but this caution is absent from the abstract and introduction. I also think that the limits of the BAR-1 experiments should be explicitly discussed.

In terms of modelling, I found it surprising that, in contrast to the rest of the paper, here the experiments are not compared to the model (apart from a brief reference to the main conclusion of an earlier publication). It seems to me that the experimental results raise questions that are not addressed by these earlier publications. For instance, DeltaN-BAR-1 animals show no increase in mig-1 levels, indicating that Wnt signaling is completely saturated. Would the model predict an increase in variability in timing (and position) when Wnt signaling is saturated? That seems to me, at first sight, a different kind of perturbation than removing positive feedback, e.g. by removing mig-1's POP-1/TCF sites. I would think that for saturated Wnt signaling, mig-1 expression becomes less dependent on the activator level, and thus less susceptible to variability in the activator level.

- l. 405-412. Unless there are very many transcription factors that bind these sequences, it seems a bit odd to not mention the identity of these transcription factors.

- l. 516-517. Apparently, some parameters, like H, r0, and K, are only constrained by the experimental data in specific ratios. It would be good if the authors provide some form of explanation for this.

-l. 526. "T=2.32". What does this time correspond to in terms of cell position?

- Figure 1D. What is the GFP reporter?

- Figure 2-supplement 2. What statistical test has been used? I am also asking because from looking at the data it is not so clear that control and Δ intron A are really not significantly different. What is the P value?

- Figure 4E. The figure shows that decreasing the size ratio leads to slower mig-1 induction, but it is not clear how well it fits the data. If you split the data in Figure 4C into brackets of different ratios/rho, would it then fit the model data for the rho of that bracket? Based on Figure 4E, F, it remains open to me how well the model actually fits the data.

- l.823-824. What is this a test of? I am confused because the text mentions transcripts, but the panel is about position.

Reviewer #2 (Recommendations for the authors):

The strength of the manuscript is the overall question, the experimental data, and the simple model. But given the quality of the data, the data analysis, model fitting, and predictions could be improved or at least better presented.

Specifically, I would like to see how much variation exists between different replica experiments. Right now, all the RNA-FISH data is pooled and not separated by replica experiments. One simple way to address this would be to plot the data points separated by different symbols specific to each biological replica experiment. This comment is relevant to all figures with data.

How did the authors determine a threshold in these data sets? Applying a piecewise linear fit to the experimental data without the model enables very precise quantification of the threshold. This should be done for each replica experiment to determine the mean and standard deviation in the threshold. This should be also done for all the WT and mutant data sets to show how the expression levels and threshold are effects by different mutants.

With regard to the model, the authors need to show how the model fits to the data, by overlaying the model and the data in the same plot. The authors state in the manuscript that the spatial position of the cells is proportional to time, so this should be possible. This should be done for each replica to determine the mean and standard deviation in the model parameter. The model parameters and their errors should be added to the manuscript.

By comparing the model fits and predictions using the estimated parameter errors, authors can better determine if the differences seen in the WT and mutant experiment are significantly large to explain this difference by changing specific rates in the model.

Although the authors quantify the position equals temporal variance in the mig1 expression, they do not use their model to make this prediction or even analyze the data with their model. Why is this? Can the model not compute temporal variances? Please clarify.

https://doi.org/10.7554/eLife.82675.sa1

Author response

Essential revisions:

1. Both reviewers raise the lack of integration between the model and data.

We now provide additional supplemental figures in which the quantitative gene expression data is superimposed on the different model fits.

2. Secondly, the evidence for the positive feedback loop is weak/lacking. Additional experiments would be optimal. Alternatively, modelling could be used to interpret or strengthen these results, but this would require explicitly including the positive feedback loop and stochastic variability. At least it would be important to substantially rephrase the conclusions, to recognize that the evidence is relatively weak. I refer you to the comments of the reviewers for details.

We now explicitly model the positive feedback and compare our findings to the experimental data. We find that when mig-1 is subjected to feedback, both through the activator and independently from the activator, then disrupting the feedback either through loss or constitutive activation leads to increased variability in the timing of mig-1 expression, which is consistent with the experiments.

Reviewer #1 (Recommendations for the authors):

- l. 151. "2C". Given the article's mixed audience of biologists and modelers, it would be good to explain the 2C, 4C, and 8C notation. Currently, it is not easy to deduce the meaning of 2C for the article text, if you do not already know it.

We have rewritten this part to explain in more detail the terminology and implications of the differences in DNA content.

- l. 154. It is not explained why it is expected that CDK1-inhibited neuroblasts are 4C, not 2C. This is presumably because inhibiting CDK1 does not inhibit DNA duplication but does inhibit the subsequent division. It would be good to mention this explicitly.

Yes, that is correct. We now explain this point more explicitly.

- l. 189-208. In general, the fitting procedures outlined in this section and the Materials and methods seem fine, but it would be helpful if they would be presented in a way that allows visual comparison of the quality of the fit. Currently, smFISH data is always shown separately from fits, and plots with fitting data are mostly used to highlight predictions of the model. It would be good if somewhere, either in the main text or in the SI, the fitting data is plotted on top of the data, just to see how good the fit actually is. I realize that data and fit have different time axes, but this should not be a fundamental problem in plotting data and fit together if I understand the fitting procedure correctly.

We have added supplemental figures showing combinations of the data and corresponding fitting (see Figure 2, figure supplement 3 and Figure 4, figure supplement 1).

- l. 189-208, l. 508-533. The fitted parameter values are only mentioned in the Methods. The values for the Hill coefficients for both the repressor and activator model are very high, H=14 or H=17, but this is not really commented on. I understand from the rest of the paper that mig-1 likely has some form of positive feedback on its own expression, which could in principle explain high Hill coefficients. But are such high Hill coefficients consistent with such a positive feedback mechanism? I think the authors should discuss the high Hill coefficient and whether such high values are plausible based on their proposed mechanism.

The high value of the Hill coefficient is a reflection of the fact that mig-1 expression remains low in QR.p cells, then increases very sharply in QR.pa cells. Because we model this increase with a single activating species, whereas mig-1 is likely regulated by more than one species in reality, the Hill coefficient can be viewed as an effective parameter rather than a biochemical constant. The reviewer is correct that when we include feedback of mig-1 on its own expression, as discussed below, we find that the Hill coefficient is lowered (from 14 or 17, to 12).

- l. 228, Figure 3A, C. The model predicts that decreasing the activator (hereby mutating activator binding sites) leads to a slower increase in mig-1 level. In the experiments, this is not visible for Δ upstream 1, but here a substantial number of animals appear to fail in inducing mig-1. For Δ upstream 2, there are even more animals like this. This fact is mentioned in the main text, but not interpreted in any way. How do the authors interpret this in light of the activator model?

The expression of mig-1 is dependent on two, partially redundant cis-regulatory elements in the upstream region of the gene. These two elements may act cooperatively in mig-1 activation, and therefore deletion of the individual elements may increase the threshold at which mig-1 expression is induced by the activator. The small population of QR.pa cells that fail to upregulate mig-1 expression in upstream deletion #1 and the larger population in upstream deletion #2 may represent cells in which this threshold is not reached within the time-window in which the cells were examined.

- l. 245. 'required for the rapid upregulation'. To me, this formulation seems too strong. Looking at Figure 4A, I still see rapid upregulation in pig-1 mutants, as in happening quickly after the division of QR.p, but with lower mig-1 levels reached compared to the control.

We agree and have adapted the text.

- l. 245-298, Figure 4. Related to the above point, it is not clear to me how the authors interpret the pig-1 experiments in relation to the fitting results in Figure 4E. The model predicts that decreasing the size ratio between QR.pa and QR.pp leads to a slower accumulation of mig-1, but with the same mig-1 level ultimately reached, albeit at a later time. The pig-1 data in Figure 4A could be consistent with slower accumulation, but, in contrast to Figure 4E, mig-1 expression never reaches the peak control level. How do the authors interpret this contrasting result? An alternative model to explain the data could be that in pig-1 mutants mig-1 expression is induced at the same time as in control animals, but with a lower 'amplitude' of mig-1 expression. Can the authors rule out such an alternative model? Finally, in the author's model QR.pa migration stops once a specific mig-1 level is reached. In that case, wouldn't one expect a positional defect for pig-1 mutants? I didn't see this mentioned anywhere in the article.

Our modeling suggests that a decrease in cell size ratio will lower the amount of activator that is segregated into QR.pa, and thereby delay the activation of mig-1 expression. Since we measure in the same time-window as in the control animals, we see this effect as a reduction in expression. We would need measurements at later time-points to confirm that mig-1 expression will further increase as predicted by the model, but this is not possible because QR.pa divides shortly after our analysis time-window.

The alternative model raised by the reviewer is that pig-1 may directly influence the expression of mig-1. Although we cannot formally exclude this possibility, our observation that mig-1 expression correlates with the degree of cell size asymmetry in pig-1 mutants strongly supports our model that segregation of the activator, and not pig-1 itself, determines the observed effects on mig-1 expression.

As to the question on QR.pa position in pig-1 mutants, we did not examine this because the reduction in QR.pa cell size may affect its ability to migrate, making these results difficult to interpret.

- l. 281-285. I found the explanation here confusing. I think I understand the main point: the cells may have different sizes, but have the same nuclear volume. As the activator is likely in the nucleus, one can therefore just use the same molecular level as before, without taking into account the difference in volume. However, the way it is written here suggested initially to me that a modification was made to the model by focusing on a number of molecules rather than concentration, even though, if I understand correctly, the model was formulated in terms of a number of molecules from the beginning. It would be good to clarify this.

Yes, all the modeling uses molecule number, not concentration. This point has been added to the main text.

- l. 296-298. I have the same issue here as above: it is not clear to me that division asymmetry is required for rapid upregulation. The proposed role of asymmetry can also be formulated more sharply. If I understand it correctly, according to the model, the only impact of size asymmetry is in the asymmetric inheritance of the activator. In particular, the subsequent rate of production of the activator in Q.pa does not depend on its size, correct? Or am I missing something?

We agree and have adapted the text.

- l. 300-329. I find this section the weakest of the manuscript, both in terms of experiments and modelling. The main aim of this section is to test the model prediction from an earlier paper (Gupta 2020) that positive feedback of mig-1 on its own expression reduces variability in timing. The main experimental weakness is that the authors do not specifically perturb this positive feedback loop, but instead use two Wnt mutants (the loss-of-function bar-1(0) and constitutively active DeltaN-BAR-1) that likely impact many genes, including those responsible for migration (l. 95,96), as also evidenced by the resulting change in average Q.pa position. While the authors find an increase in variability of Q.pa position, this could also be explained by changes in the expression of any of the other BAR-1 targets. The proper experiment, I think, would be to remove the POP-1/TCF sites in the mig-1 promoter, although there might be experimental reasons why this is not feasible. The last sentence of the section (l. 328-329) is a bit careful in drawing strong conclusions from these result, but this caution is absent from the abstract and introduction. I also think that the limits of the BAR-1 experiments should be explicitly discussed.

We agree that mutating POP-1/TCF binding sites in the mig-1 promoter would more specifically address this question. Unfortunately, there are 12 predicted POP-1/TCF binding sites within the first 1 kb of upstream sequence, making it challenging to mutate all of these sites in the endogenous locus.

In terms of modelling, I found it surprising that, in contrast to the rest of the paper, here the experiments are not compared to the model (apart from a brief reference to the main conclusion of an earlier publication). It seems to me that the experimental results raise questions that are not addressed by these earlier publications. For instance, DeltaN-BAR-1 animals show no increase in mig-1 levels, indicating that Wnt signaling is completely saturated. Would the model predict an increase in variability in timing (and position) when Wnt signaling is saturated? That seems to me, at first sight, a different kind of perturbation than removing positive feedback, e.g. by removing mig-1's POP-1/TCF sites. I would think that for saturated Wnt signaling, mig-1 expression becomes less dependent on the activator level, and thus less susceptible to variability in the activator level.

We now explicitly model the positive feedback and compare our findings to the experiments. Specifically, we find that when mig-1 is subjected to feedback, both through the activator and independently from the activator, then disrupting the feedback either through loss or constitutive activation leads to increased variability in the timing of mig-1 expression, consistent with the experiments (see the new Figure 5B, Figure 5, figure supplement 2, and Materials and methods section “Modeling the effect of feedback on timing precision”).

These new results account for saturation of Wnt signaling in ΔN-BAR-1 animals by making the activator saturated at a constant, time-independent level in this case. The mig-1 expression level then indeed becomes effectively independent of the activator, and thus not subject to variability in the activator level. However, it also loses the benefit of the activator dynamics to create the sharp increase in mig-1 expression, which acts to reduce variability. The net result, the modeling shows, is that the mig-1 dynamics become more linear and more variable, consistent with the experiments.

- l. 405-412. Unless there are very many transcription factors that bind these sequences, it seems a bit odd to not mention the identity of these transcription factors.

Transcription factors that bind to these regions include homeobox transcription factors such as LIN-39 and MAB-5 that regulate Q neuroblast descendant migration. This has been added to the main text.

- l. 516-517. Apparently, some parameters, like H, r0, and K, are only constrained by the experimental data in specific ratios. It would be good if the authors provide some form of explanation for this.

For the activator model, plugging $a (t) = k t$ into $F (t) = α a^{H} / (a^{H} + K^{H})$ obtains $F (t) = α / [1 + {(K / k t)}^{H}]$ , which depends only on the ratio $K / k$ . For the repressor model, plugging $r (t) = r_{0} e^{- μ t}$ into $F (t) = α K^{H} / (r^{H} + K^{H})$ obtains $F (t) = α / [1 + {(r_{0} e^{- μ t} / K)}^{H}] = α / {1 + e x p [H l n (r_{0} / K) - μ H t]}$ , which depends only on the combinations $H l n (r_{0} / K)$ and $μH$ .

-l. 526. "T=2.32". What does this time correspond to in terms of cell position?

It is the position at which the number of QR.pa data points to the right equals the number of QR.p data points to the left. Quantitatively, it is 36% of the way from the V2 seam cell to the V3 seam cell (because $P = P_{max} - T = 4.68 - 2.32 = 2.36$ in units of the seam cells, where $P_{max} = 4.68$ is the position of the rightmost QR data point).

- Figure 1D. What is the GFP reporter?

The GFP report is expressed in the seam cells and the QR lineage descendants (Pwrt-2::gfp, transgene heIs63). This has been added to the figure legend.

- Figure 2-supplement 2. What statistical test has been used? I am also asking because from looking at the data it is not so clear that control and Δ intron A are really not significantly different. What is the P value?

Welch’s t-test, as in Figure 3, supplemental figure 2. The p-value is p<0.0001. This has been added to the figure legend.

- Figure 4E. The figure shows that decreasing the size ratio leads to slower mig-1 induction, but it is not clear how well it fits the data. If you split the data in Figure 4C into brackets of different ratios/rho, would it then fit the model data for the rho of that bracket? Based on Figure 4E, F, it remains open to me how well the model actually fits the data.

See new Figure 4, supplemental figure 1 for combination of experimental data and model fit.

- l.823-824. What is this a test of? I am confused because the text mentions transcripts, but the panel is about position.

We performed a Brown-Forsythe analysis on the coefficients of variation (CV) at different thresholds of mig-1 expression. See Materials and methods section “Statistical analysis of variance data” for more details.

Reviewer #2 (Recommendations for the authors):

The strength of the manuscript is the overall question, the experimental data, and the simple model. But given the quality of the data, the data analysis, model fitting, and predictions could be improved or at least better presented.

Specifically, I would like to see how much variation exists between different replica experiments. Right now, all the RNA-FISH data is pooled and not separated by replica experiments. One simple way to address this would be to plot the data points separated by different symbols specific to each biological replica experiment. This comment is relevant to all figures with data.

The smFISH measurements were performed in synchronized populations enriched for QR, QR.p or QR.pa, respectively. Combined, these different experiments show the full range of mig-1 expression. The data of the individual replicate experiments of Figure 1A is shown in Figure 1, figure supplement 1.

How did the authors determine a threshold in these data sets? Applying a piecewise linear fit to the experimental data without the model enables very precise quantification of the threshold. This should be done for each replica experiment to determine the mean and standard deviation in the threshold. This should be also done for all the WT and mutant data sets to show how the expression levels and threshold are effects by different mutants.

By threshold we take the reviewer to mean the activator level $K$ above which the mig-1 production rate increases. A piecewise linear fit to the mig-1 data would indeed quantify $K$ precisely. However, a piecewise linear fit would be equivalent to taking the Hill coefficient to $H \to \infty$ , which we do not want to assume. Therefore, we fit $K$ and $H$ simultaneously. Importantly, the value of $K$ fit with $H \to \infty$ would be different than the value of $K$ fit together with a finite $H$ .

With regard to the model, the authors need to show how the model fits to the data, by overlaying the model and the data in the same plot. The authors state in the manuscript that the spatial position of the cells is proportional to time, so this should be possible. This should be done for each replica to determine the mean and standard deviation in the model parameter. The model parameters and their errors should be added to the manuscript.

We have added figures that overlay the experimental data and model fits (Figure 2, figure supplement 3 and Figure 4, figure supplement 1).

By comparing the model fits and predictions using the estimated parameter errors, authors can better determine if the differences seen in the WT and mutant experiment are significantly large to explain this difference by changing specific rates in the model.

In the majority of cases, the mutants produce changes in the data that go beyond the spreads in the data points themselves. Thus, it is unlikely that these changes could be encompassed within parameter uncertainties in the model, even if these uncertainties could be estimated from replicates. This observation justifies treating the model as a phenomenological description whose predictions manifest as one investigates changes in its parameters.

Although the authors quantify the position equals temporal variance in the mig1 expression, they do not use their model to make this prediction or even analyze the data with their model. Why is this? Can the model not compute temporal variances? Please clarify.

https://doi.org/10.7554/eLife.82675.sa2

Article and author information

Author details

Erik S Schild

Hubrecht Institute, Royal Netherlands Academy of Arts and Sciences and University Medical Center Utrecht, Utrecht, Netherlands

Contribution
Conceptualization, Data curation, Software, Formal analysis, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared
Shivam Gupta

Department of Physics and Astronomy, Purdue University, West Lafayette, United States

Contribution
Conceptualization, Investigation

Competing interests
No competing interests declared
Clément Dubois

Institut de Biologie de l'Ecole Normale Supérieure, CNRS, INSERM, Paris, France

Contribution
Investigation, Methodology

Competing interests
No competing interests declared
Euclides E Fernandes Póvoa

Hubrecht Institute, Royal Netherlands Academy of Arts and Sciences and University Medical Center Utrecht, Utrecht, Netherlands

Contribution
Investigation, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-0920-7783
Marie-Anne Félix

Institut de Biologie de l'Ecole Normale Supérieure, CNRS, INSERM, Paris, France

Contribution
Conceptualization, Supervision, Funding acquisition, Project administration, Writing – review and editing

Competing interests
No competing interests declared
Andrew Mugler
1. Department of Physics and Astronomy, Purdue University, West Lafayette, United States
2. Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, United States
Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Project administration, Writing – review and editing

For correspondence
andrew.mugler@pitt.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9367-7026
Hendrik C Korswagen
1. Hubrecht Institute, Royal Netherlands Academy of Arts and Sciences and University Medical Center Utrecht, Utrecht, Netherlands
2. Institute of Biodynamics and Biocomplexity, Department of Biology, Utrecht University, Utrecht, Netherlands
Contribution
Conceptualization, Data curation, Supervision, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing

For correspondence
r.korswagen@hubrecht.eu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7931-4472

Funding

Human Frontier Science Program (RGP0030/2016)

Marie-Anne Félix
Andrew Mugler
Hendrik C Korswagen

National Science Foundation (PHY-1945018)

Shivam Gupta
Andrew Mugler

Simons Foundation (376198)

Shivam Gupta
Andrew Mugler

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

Acknowledgements

We thank the Hubrecht Imaging Center for technical support. This work was funded by the Human Frontier Science Program Grant RGP0030/2016 to AM, MAF, and HCK and grants from the Simons Foundation (376198) and the National Science Foundation (PHY-1945018) to SG and AM. Some strains were provided by the CGC, which is funded by the NIH Office of Research Infrastructure Programs (P40 OD010440). We are grateful to Damien Coudreuse and members of the Korswagen and Galli groups for helpful discussions and critical reading of the manuscript.

Senior and Reviewing Editor

Marianne E Bronner, California Institute of Technology, United States

Version history

Received: August 12, 2022
Preprint posted: October 26, 2022 (view preprint)
Accepted: May 12, 2023
Accepted Manuscript published: May 15, 2023 (version 1)
Version of Record published: June 7, 2023 (version 2)

Copyright

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