Canalization of C. elegans Vulva Induction against Anatomical Variability.

It is a fundamental open question as to how embryos develop into complex adult organisms with astounding reproducibility, particularly because cells are inherently variable on the molecular level. During C. elegans vulva induction, the anchor cell induces cell fate in the vulva precursor cells in a distance-dependent manner. Surprisingly, we found that initial anchor cell position was highly variable and caused variability in cell fate induction. However, we observed that vulva induction was "canalized," i.e., the variability in anchor cell position and cell fate was progressively reduced, resulting in an invariant spatial pattern of cell fates at the end of induction. To understand the mechanism of canalization, we quantified induction dynamics as a function of anchor cell position during the canalization process. Our experiments, combined with mathematical modeling, showed that canalization required a specific combination of long-range induction, lateral inhibition, and cell migration that is also found in other developmental systems.

Introduction

In a developing embryo, each cell has to assume the correct cell fate in order to give rise to a viable adult organism. Embryonic development is highly reproducible on the organismal level, as is evident in the strikingly similar appearance of identical twins. Yet, it has become clear that many fundamental biological processes, such as gene expression and cell signaling, are stochastic, causing significant variability between otherwise identical cells (Raj and van Oudenaarden, 2008). How highly complex bodies can be built in such a robust manner (i.e., without errors) out of cells that individually show strong intrinsic variability is a fundamental and unsolved question.

An important concept to explain the robustness of development is canalization (Waddington, 1942, Waddington, 1957). Whereas there exist many ways in which biological systems exhibit robustness (Stelling et al., 2004), canalization refers to a specific form of robustness, namely to variability in the initial conditions. In this picture, canalization is an active mechanism that causes a decrease in the variability of developmental processes as development progresses, thereby ensuring an identical outcome despite initial variation between individuals. This decrease in variability over time is thought to be due to the action of the gene regulatory networks that underlie development (Waddington, 1957, Ferrell, 2012, Huang, 2012). However, it is an open question whether variability in initial conditions is indeed a strong source of variability in developmental systems and, hence, to what extent canalization is responsible for the exceeding robustness of development. Moreover, due to the limited number of systems where canalization has been studied on the molecular level (Manu et al., 2009a, Manu et al., 2009b, Gavin-Smyth et al., 2013), it is not known what molecular and cellular mechanisms give rise to it.

A prime example of a robust developmental process is vulva induction in the nematode C. elegans (Sternberg, 2005). The C. elegans vulva forms from a row of vulva precursor cells (VPCs), labeled P3.p–P8.p. During vulva induction, a specific spatial pattern of cell fates, distinguishable by lineage, is induced in a manner that depends only on the relative distance of each VPC to the anchor cell (AC) (Figure 1A): P6.p, the cell adjacent to the AC, assumes 1° fate, the more distant P5.p and P7.p cells assume 2° fate, and the remaining P3.p, P4.p, and P8.p cells assume 3° fate. The resulting cell fate pattern is highly robust. Under standard laboratory conditions, <1% of animals show minor deviations (Braendle and Félix, 2008, Grimbert and Braendle, 2014). Such deviations include centering shifts, in which the correct 2°-1°-2° fate pattern is still induced and results in a functional vulva, but with 1° fate assumed by P5.p or P7.p rather than P6.p. How this observed robustness of the vulva cell fate pattern arises is much-studied but not yet understood (Félix and Barkoulas, 2012).

Figure 1

Correction of AC Misplacement

(A) Overview of vulva induction. LIN-3 from the AC (magenta arrows) induces 1° fate and expression of Notch ligands lag-2 and apx-1 in P6.p. Subsequent Notch signaling (green arrows) inhibits 1° fate and induces 2° fate in P5.p and P7.p.

(B) Definition of the relative AC position R. Cell positions are given by the position of the cell nucleus. Nuclei are visualized by DAPI. Green lines show the approximate outlines of the AC and VPC bodies and magenta circles indicate the VPC nuclei. The other nuclei inside the outline of the VPC belong to adjacent neurons.

(C) Examples of the different observed classes of AC position: correctly placed (red, ), mildly misplaced (green, ) and severely misplaced (blue, ).

(D–F) Relative AC position R as a function of time during vulva induction in (D) wild-type animals, (E) dig-1 mutants, and (F) dig-1;lin-3(e1417) mutants. Each marker corresponds to the relative AC position in an individual animal and the color of the marker corresponds to the degree of AC misplacement as defined in (C). The black line indicates the moving average of relative AC position and the gray area indicates the SD, both with a window size of 2 hr.

(G–I) The fraction of animals with correctly placed AC (red), mildly misplaced AC (green), and severely misplaced AC (blue) as a function of time for (G) wild-type animals. (H) dig-1 animals and (I) dig-1;lin-3(e1417) animals. The difference between dig-1 and dig-1;lin-3 animals is significant (Kolmogorov-Smirnov test, p < 10−5).

The signaling network that controls VPC fate induction is well-characterized, making it uniquely suited to study the molecular mechanisms underlying its robustness (Sternberg, 2005). VPC fate is determined by a LIN-3/epidermal growth factor (EGF) signal produced in the AC, in combination with Notch signaling between neighboring VPCs. It is thought that LIN-3/EGF forms a long-range spatial gradient, which induces Ras signaling in VPCs in a graded manner, with Ras activation strongest in P6.p and weaker in the neighboring P5.p and P7.p cells (Yoo et al., 2004, van Zon et al., 2015). Subsequently, Ras signaling induces expression of Notch ligands (Chen and Greenwald, 2004, Zhang and Greenwald, 2011, van Zon et al., 2015). These activate Notch receptors in neighboring VPCs, causing an inhibition of Ras signaling. Hence, P6.p, the VPC that receives the strongest EGF input, will come to fully inhibit the response to EGF in P5.p and P7.p. As a result, the EGF gradient is amplified into an all-or-none difference in signaling and cell fate between VPCs, with high Ras activity in P6.p (1° fate) and high Notch activity in P5.p and P7.p (2° fate).

In this model, the position of the AC is crucial for the establishment of the VPC fate pattern: it is thought that P6.p assumes 1° fate and induces 2° fate in P5.p and P7.p because it is located directly adjacent to the AC. Surprisingly, we found that at the start of vulva induction the position of the AC relative to the VPCs showed strong animal-to-animal variability, with the AC often positioned in between P5.p and P6.p. In addition, we found that variability in AC position gave rise to significant variability in expression patterns of 1° cell fate markers between animals, with adjacent VPCs often showing similar 1° fate induction levels in animals with a misplaced AC. Yet, we found that vulva induction was canalized: the initial variability in AC position and 1° fate induction decreased over time, ultimately resulting in the same stereotypical cell fate pattern for all animals, with a single 1° fate cell directly adjacent to the AC. By quantifying Notch ligand expression as a measure of 1° fate induction and Notch signaling, and by comparing our results to mathematical models of the induction process, we identified the combined action of (1) graded EGF signaling, (2) lateral Notch inhibition, and (3) EGF-induced VPC movement toward the AC as the key requirements for the observed canalization of AC misplacement and 1° fate induction.

Results

Variability in AC Position and 1° Fate Induction Is Corrected during Vulva Induction

We measured the position of the AC relative to the VPCs in fixed wild-type (N2) animals using DAPI staining to visualize cell nuclei. We defined the relative AC position as R = 2ΔA/ΔL, where ΔA is the distance between the AC and the closest VPC and ΔL the distance between two closest VPCs (Figure 1B). Distances were measured along the body axis and the cell position was defined as the center of its nucleus. A relative AC position of R = 0 corresponded to an AC correctly positioned adjacent to the closest VPC, whereas R = 1 corresponded to a maximally misplaced AC, i.e., equidistant to the two closest VPCs (Figure 1C). We choose a definition of R that does not depend explicitly on the position of P6.p, as the variability in AC position caused the exact identity of the closest and second closest VPC to vary between individuals, particularly in the mutants studied further below. For each fixed animal, we determined the stage of vulva induction by measuring the gonad length, defined as the distance between the two distal tip cells (DTCs). Gonad length increases during vulva induction (Kimble and Hirsh, 1979) and serves as a reliable measure of the time relative to the start of the L2 larval stage (van Zon et al., 2015) (see the STAR Methods). We found that in wild-type animals during the early induction stage (<3 hr into the L2 stage) the AC position was highly variable with a significant fraction of animals showing a severely misplaced AC (; Figures 1D, 1G, and 2B). In those cases, the AC was often positioned between P5.p and P6.p (Figure S1). However, during induction the variability in AC position and the average degree of AC misplacement decreased, with all animals having at the end of induction (Figures 1D and 1G). These observations suggested that the process of vulva induction is able to correct for significant deviations in the AC position. Similar correction of AC misplacement was observed recently in Grimbert et al. (2016), both in C. elegans and other nematode species, suggesting it is crucial for vulva induction.

Figure 2

Examples of AC Position and Notch Ligand Expression during Vulva Induction

(A and B) Representative examples of expression patterns of the Notch ligands apx-1 (green) and lag-2 (red) in VPCs for wild-type animals with (A) a correctly placed AC and (B) a severely misplaced AC, at different stages of vulva induction. Single apx-1 and lag-2 mRNA molecules are visualized as diffraction-limited spots using single-molecule fluorescence in situ hybridization (smFISH) and cell nuclei are stained by DAPI (blue). Time corresponds to hours after start of L2. For the late induction stage (>9 hr), we did not observe animals with a severely misplaced AC. Notch ligands are not only visible in VPCs, but also in other cells in the gonad, including the AC.

(C and D) Typical expression patterns in dig-1 mutants with (C) a correctly positioned AC and (D) a severely misplaced AC.

(E) Typical expression patterns in dig-1;lin-3(e1417) mutants with a severely misplaced AC. Due to the lack of vulva induction in this mutant, apx-1 and lag-2 expression in VPCs is strongly reduced and hence not visible.

To examine the mechanism underlying the correction of AC misplacement, we quantified the expression level of the Notch ligands apx-1 and lag-2 in the VPCs as a function of the degree of AC misplacement. Expression of these two Notch ligands is induced by the EGF/Ras signaling pathway and is a measure of 1° fate induction (Chen and Greenwald, 2004, Zhang and Greenwald, 2011, van Zon et al., 2015). We quantified Notch ligand expression using single-molecule fluorescence in situ hybridization (smFISH) to visualize and count individual mRNA molecules in fixed animals (Raj et al., 2008, Ji and van Oudenaarden, 2012). In general, we observed different patterns of Notch ligand expression depending on the degree of AC misplacement. For animals with a correctly placed AC, R ≈ 0, we saw apx-1 expressed predominantly in P6.p and at low levels (Figure 2A, left and middle panel). At this early stage, lag-2 expression was very low. At the late induction stage both apx-1 and lag-2 were expressed in P6.p at much higher levels (Figure 2A, right panel). We observed a strikingly different expression pattern in animals with a severely misplaced AC, R ≈ 1. Here, at the early induction stage, we frequently observed animals with nearly equal levels of apx-1 expression in the two closest VPCs (Figure 2B, left panel), likely reflecting that they received similar levels of EGF input. However, in older animals with R ≈ 1, we observed that apx-1 expression was restricted to only one of the two VPCs (Figure 2B, right panel). Finally, we did not observe older animals with misplaced ACs (; Figures 1D and 1G).

During wild-type development, severe AC misplacement occurred in a significant but still limited fraction of animals. To analyze the response of vulva induction to severe AC misplacement, we therefore sought to increase the variability in AC position. In most dig-1(n1321) animals, the gonad, which contains the AC, is shifted anteriorly, with the AC typically positioned most closely to P5.p (Thomas et al., 1990). In such cases, the vulva is properly induced, but often centered on P5.p instead of P6.p (Figure S1). We found that AC position was highly variable in dig-1 mutants during early induction, with many animals that have (Figures 1E and 1H). In dig-1 animals with P5.p correctly aligned with the AC, apx-1 and lag-2 were typically expressed only in P5.p, with dynamics similar to that observed in P6.p in wild-type animals (Figure S2). We found that in dig-1 mutants, despite the increased frequency of animals with severe AC misplacement at early induction, the degree of AC misplacement decreased as vulva induction progressed, similar to wild-type animals (Figures 1E and 1H). In particular, we found no animals with at the end of induction (>10 hr after start of L2). Additionally, we found for dig-1 animals the same progression of Notch ligand expression patterns as a function of AC misplacement that we observed in wild-type animals (Figures 2C and 2D). Therefore, we analyzed the response to AC misplacement in dig-1 rather than wild-type animals, focusing on apx-1, the earliest expressed Notch ligand.

To systematically examine the time dynamics of apx-1 expression as a function of the relative AC position R, we quantified the relative apx-1 expression level E, defined as E = (M1 − M2)/(M1 + M2). Here, Mi indicates the number of apx-1 mRNAs in the closest (M1) and second closest VPC (M2) with respect to the AC. With this definition, E = 1 corresponds to an animal in which all apx-1 mRNAs are expressed in the closest VPC, E = −1 corresponds to all apx-1 mRNAs expressed in the more distant VPC, and E = 0 corresponds to equal apx-1 mRNA levels. In Figure 3A, we plot for each animal the relative apx-1 expression level as a function of time, with the color of the marker indicating the degree of AC misplacement as shown in Figure 1C. For animals at the start of induction we frequently observed equal levels of apx-1 expression in both closest VPCs (, Figure 3A). This symmetric apx-1 expression pattern correlated with R, occurring more frequently in animals with severely misplaced AC (, Figure 3C). However, the fraction of animals exhibiting apx-1 expression in both VPCs decreased rapidly over the course of induction, with apx-1 expression restricted to the closest VPC in all animals at the end of induction (, Figures 3A and 3B).

Figure 3

Notch Ligand Expression Dynamics as a Function of AC Position

(A) Relative apx-1 expression level E as a function of time in dig-1 animals. Each marker represents the apx-1 expression pattern in an individual animal and is colored according to the degree of AC misplacement as shown in Figure 1C.

(B) Fraction of animals exhibiting the different classes of apx-1 expression patterns indicated schematically in (A): most apx-1 expression in the closest VPC (cyan, ), symmetric apx-1 expression in the two closest VPCs (purple, ) and most apx-1 expression in the more distant VPC (magenta, ), as a function of the time of induction.

(C–E) Number of animals with relative apx-1 expression level E and relative AC position R for the early (C, 0–3 hr), intermediate (D, 3–9 hr), and late (E, 9–12 hr) induction stage.

We compared the relative timing of the decrease in AC misplacement with that of the restriction of apx-1 expression to a single VPC, by plotting the fraction of animals found for each combination of R and E at the early, middle and late stage of vulva induction (Figures 3C–3E and S3). In general, restriction of apx-1 expression appeared to precede the correction of AC misplacement, particularly for animals with misplaced ACs (R ≈ 1, E ≈ 0, Figures 3C and 3D), with full correction of AC misplacement only observed subsequently at the late induction stage (Figure 3E). However, we also observed a significant simultaneous decrease in AC misplacement from the early to the middle induction stage (Figures 1H, 3C, and 3D), suggesting that the correction of AC position itself could also contribute to canalization of 1° fate induction, by bringing one VPC closer to the AC and thereby increasing the amount of LIN-3/EGF signal it receives.

Induced 1° Fate VPCs Move toward the AC

The expression dynamics in Figure 3 suggested that active correction of AC position could be important for the canalization of the observed variability in 1° fate induction. However, animals grow significantly during vulva induction and the observed change in AC position could instead be due to divisions and rearrangements of cells surrounding the AC, independent of the induction process. To test whether the correction of AC position depended on vulva induction, we measured the relative AC position in dig-1;lin-3(e1417) mutant animals. In lin-3(e1417) mutants, lin-3/EGF expression is reduced specifically in the AC, leading to the loss of 1° and 2° fate induction (Hwang and Sternberg, 2004) and to the absence of apx-1 and lag-2 expression in the VPCs (Figure S2). In dig-1;lin-3(e1417) animals, we observed the same distribution of R during early induction as in dig-1 animals, but we did not see a significant reduction in AC misplacement over time (Figures 1F and 1I). In particular, we found many dig-1;lin-3(e1417) animals with , at the end of induction (Figure 2E). Hence, the observed correction of AC misplacement was the result of an active process that depended on the presence of the LIN-3/EGF signal. This link between AC position and LIN-3/EGF signaling was also found independently in Grimbert et al. (2016). In addition, our experiments indicated that this is controlled upstream of the Ras target LIN-1 (Figure S4).

The observed alignment between a single 1° fate VPC and AC at the late induction stage raised the question whether this alignment depended specifically on 1° fate induction and if so, whether the 1° fate VPC moved toward the AC or, instead, the AC moved toward the 1° fate VPC. To address this question, we examined the positions of 1° fate VPCs with respect to the AC in a lin-12(0) mutant, where multiple VPCs assume 1° fate (Greenwald et al., 1983). In this mutant, the LIN-12/Notch receptor that mediates Notch signaling between adjacent VPCs is not functional. In addition, as lin-12 is also involved in AC specification, lin-12(0) mutants typically have two ACs (Greenwald et al., 1983), likely leading to increased LIN-3/EGF levels. These combined effects cause not only P6.p but also P5.p and ∼60% of P7.p cells to assume 1° fate (Sternberg and Horvitz, 1989). We measured the position of P(5–7).p along the body axis as a function of time, in both wild-type and in lin-12(0) animals (Figure 4C). We observed that during the late vulva induction stage, ∼90% of P5.p cells and ∼60% of P7.p cells were located significantly closer to the ACs (>9 hr, Figures 4A–4C) and P6.p (Figure S5G) compared with wild-type animals. We examined whether this displacement of P5.p and P7.p depended on 1° fate induction. Indeed, we found that displacement toward the AC correlated strongly with the expression level of the 1° fate marker apx-1 (Figure 4I). This correlation was particularly striking for P7.p: cells with low apx-1 expression were often observed at wild-type distance to the AC (Figures 4A and 4I), while cells with high expression were positioned much closer (Figures 4B and 4I). Moreover, in lin-12(0) animals ectopic expression of apx-1 was observed in P5.p hours before that cell moved closer to P6.p, showing that 1° fate induction preceded VPC migration (Figures S4D and S4E).

Figure 4

Movement of 1° Fate VPCs to the AC

(A) Example of a lin-12(0) mutant animal with two induced VPCs at the late vulva induction stage. Expression of apx-1 (green) and lag-2 (red) is visualized using smFISH. Nuclei (blue) are stained by DAPI. lin-12(0) animals have two ACs. The positions of the ACs and VPC nuclei are indicated.

(B) Example of a lin-12(0) mutant animal with three induced VPCs at the late induction stage.

(C) Distance of the P5.p and P7.p cells to the ACs in wild-type animals (n = 287, gray squares for P5.p and gray diamonds for P7.p) and lin-12(0) mutant animals (n = 112, red circle for P5.p and blue circle for P7.p) as a function of the time of induction. Each marker corresponds to an individual animal. Also shown are the moving averages, with a window size of 1 hr, of the distance of the P5.p and P7.p cell to the P6.p cell for wild-type animals (black lines).

(D) Distribution of distance of P5.p (red) and P7.p (blue) to the AC as a function of time in the gradient-sensing (G-S) model without Notch signaling. The distribution was calculated for n = 103 simulations with LIN-3/EGF level p0 = 1.37 and a random initial relative AC position between R0 = 0 (centered on P6.p) and R0 = 1 (positioned equidistant to P5.p and P6.p). Also shown are trajectories for individual simulations with R0 = 0.2 (solid line) and R0 = 0.8 (dashed line) Distances are given in units of the length L0 of undeformed VPCs.

(E and F) Steady-state model configuration for the trajectories shown in (D) with (E) R0 = 0.8 and (F) R0 = 0.2. Also shown are the AC position (black circle) and induction level (indicated in magenta).

(G and H) Distribution for steady-state distance of P5.p (red) and P7.p (blue) to the AC for (G) the lin-12(0) mutant at the late induction stage (>9 hr) and (H) the G-S model without Notch signaling. Dashed lines indicate (G) mean distance of P5.p (red) and P7.p (blue) in wild-type animals (>9 hr) and (H) distance of P5.p (red) and P7.p (blue) in the G-S model with no induction, p0 = 0.

(I) Displacement of the P5.p (red marker) and P7.p cell (blue marker) toward the ACs, as a function of the apx-1 expression level in each VPC. All animals are at the late induction stage (>9 hr, n = 34). The displacement is defined as the distance between each VPC and the ACs minus the average distance for wild-type animals at the same time of induction, as given by the black line in (C).

The general observation that in animals with multiple 1° fate VPCs these cells are located closer to another (Figures 4C and S5G) can only be explained by 1° fate VPCs moving toward the AC, rather than the reverse. This is consistent with previous observations in which the AC showed no movement with respect to seam cells and body muscles during the course of induction (Thomas et al., 1990, Grimbert et al., 2016), and when the proximal VPCs, P(5–7).p, are ablated, distal VPCs such as P8.p can assume 1° fate and move toward the AC (Sternberg and Horvitz, 1986, Grimbert et al., 2016). Surprisingly, we observed that in lin-12(0) animals in which two VPCs were induced, these were often positioned equidistant to the ACs. This shows that VPC movement by itself does not necessarily result in configurations with a single VPC aligned with the AC.

Distinct Spatial VPC Configurations in a Mathematical Model of LIN-3-Induced VPC Movement

The above results suggested that variability in AC position relative to the VPCs is canalized by the movement of 1° fate VPCs toward the AC. To examine whether VPC movement could also contributes to the canalization of variability in 1° fate induction or whether Notch signaling between VPCs is required for this, we constructed a mathematical model of LIN-3/EGF-induced VPC movement without Notch signaling. To implement movement and deformation of VPCs we used a two-dimensional “vertex model,” where cells are represented by edges that are connected in vertices (Figure 5A; Equations 1 and 2 in the STAR Methods). Such models accurately describe movement and rearrangement of cells in epithelial tissues and allow for cell deformation due to internal and external forces (Fletcher et al., 2014).

Figure 5

Mathematical Model of Canalized Vulva Induction

(A) Overview of the model without Notch signaling (see Equations 1, 2, 3, 4, 5, and 6 in the STAR Methods). VPC geometry is determined entirely by the position of the vertices where cell edges meet (blue dots). The forces produced by the VPCs (green arrows) depend on both their position relative to the AC (black circle) and their level of Ras signaling (indicated in magenta). Also shown is the apical-basal axis (Ap-Ba). L0 corresponds to the length of the undeformed VPC body along the A-P axis.

(B) The LIN-3/EGF gradient p (black line), Ras activation ϕ (magenta line), and force f (green line) as function of the distance to the AC. Black markers indicate values for the configuration in (A).

(C) Examples of different steady-state configurations obtained with the model when varying the total LIN-3/EGF level . Configurations 1 and 2 coincide with those observed in the lin-12(0) mutant while configurations 3 and 4 match those in wild-type and lin-3(e1417) animals, respectively.

(D) Diagram representing the final relative AC position Rfinal as a function of the LIN-3/EGF level p0. Lines correspond to different initial AC positions R0 = 0, 0.2, 0.4, 0.6, 0.8, 1.0 and are colored according to the initial degree of AC misplacement as shown in (E). Points labeled 1–4 correspond to the configurations in (C). The dashed lines indicate the transition regions where the final cell configuration depends strongly on the initial AC position.

(E) Similar to (D) but for the model that includes Notch signaling. The system evolves toward the desired VPC configuration with Rfinal ≈ 0 for all initial AC positions and EGF/LIN-3 levels. The colored cell diagrams correspond to the different degrees of the initial AC misplacement R0, as described in Figure 1C.

Little is known about the molecular mechanisms that control VPC movement. However, two qualitatively different general mechanisms have been proposed to explain cell migration in a chemoattractant gradient: “attractant-maximization” (A-M, Equation 3 in the STAR Methods), where cells change their position and shape to maximize the amount of attractant integrated over their surface (Savill and Hogeweg, 1997) and “gradient-sensing” (G-S, Equation 4 in the STAR Methods), where cells measure the external chemoattractant gradient and move in the direction of the largest increase in attractant concentration (Szabó et al., 2010). To implement A-M and G-S models of VPC movement, we assumed that an exponential LIN-3/EGF gradient p(x), centered at the AC, acts as chemoattractant to provide the directional cue (Figure 5B, black line; Equation 5 in the STAR Methods). In addition, we assumed LIN-3/EGF plays a second role in controlling VPC movement. Following the experimentally observed correlation between VPC movement and 1° fate induction (Figure 4I), we assumed that the propensity to migrate increases with the Ras activation level ϕ (Figure 5B, magenta line; Equation 6 in the STAR Methods), which itself depends on the absolute level of external LIN-3/EGF.

When varying the total LIN-3/EGF level , both models generate steady states that cluster around a small number of distinct VPC configurations, with sudden transitions between configurations as p0 increases (Figures 5C, 5D, and S5). These configurations resemble those observed experimentally in mutants that differ in LIN-3/EGF level: for low p0, there is neither induction nor correction of AC misplacement (Figures 5C, line 4, and 5D), similar to the lin-3(0) mutant. For intermediate p0, a single VPC is induced and aligns with the AC (Figures 5C, line 3, and 5D) as in wild-type animals. For high p0, we observe the two configurations seen in the lin-12(0) mutant with two LIN-3-expressing ACs (Figure 5C, lines 1 and 2), with either P5.p and P6.p induced and equidistant to the AC, or P(5–7).p induced and P6.p aligned with the AC. The steady-state configurations depend on the initial degree of AC misplacement (Figure 5D), reflecting that stronger VPC deformations and, hence, larger forces are required to correct for more severe AC misplacement.

To gain insight into the origin of the different configurations, we calculated the total migration force produced by each VPC as function of distance to the AC. Even though the A-M and G-S models differ strongly in cell shape dynamics (Figure S5), the expressions for their force-distance curves are identical (Equation 5 in the STAR Methods) and have a characteristic shape (Figure 5B): at sufficiently large distance from the AC, the magnitude of the force increases with the Ras activation level ϕ. However, the force peaks once either the anterior or posterior edge of the VPC aligns with the AC. When the VPC is positioned closer to the AC, the two sides of the VPC body are exposed to opposing gradients, causing the force to decrease and ultimately vanish when the VPC is centered with the AC. For intermediate p0, when only P6.p is induced, this decrease in force ensures that the cell comes to rest when aligned with the AC (Movie S1). However, this property of the force-distance curve also leads to stability of the misaligned configuration with two induced cells (Figure 5C, line 2) for larger p0. In this case, even when the AC is positioned much closer to P6.p than P5.p, P6.p will produce a lower migratory force than P5.p due to its proximity to the AC. As a result, P5.p will push P6.p away from the AC until they are equidistant and produce equal but opposite migration forces (Movie S2). For even larger p0, this configuration becomes unstable: as P6.p contracts toward the AC, P7.p is pulled closer and is also induced. P6.p and P7.p together produce sufficient force to push P5.p away from the AC, leading to a configuration with P6.p aligned with the AC and compressed by the opposing forces from P5.p and P7.p (Movie S3). Additional modeling showed that the shape of the force-distance curve in Figure 5B, with the magnitude of the force peaking at a distance to the AC comparable to the size of a VPC, is essential for reproducing the experimentally observed configuration with two VPCs equidistant to the AC (Figure S6).

For LIN-3/EGF levels at the transition between different stable cell configurations, the steady-state configuration depends strongly on the initial AC position (Figure 5D, dashed lines). Interestingly, the VPC induction and migration dynamics in these regimes can provide an appealing explanation for the broad range of VPC configurations we observed in lin-12(0) mutants (Figures 4A–4C). For instance, for p0 = 1.4, the model evolves toward the configuration with two induced VPCs equidistant to the AC for a sufficiently misplaced AC, but generates the configuration with P(5–7).p induced and P6.p aligned with the AC for more correctly positioned ACs (Figures 4D–4F; Movie S4). In lin-12(0) mutants, we also observed that induced P5.p and P7.p cells were more closely positioned to the ACs (Figures 4C and 4I) and P6.p (Figure S5G) than when not induced. Based on this observation, we can rule out the A-M model in favor of the G-S model. In the A-M model, where cells elongate along the A-P axis to maximize overlap with the LIN-3/EGF gradient, induced P5.p and P7.p cells were positioned further apart (Figure S5I). In contrast, in the G-S model P5.p and P7.p moved closer to the AC and P6.p upon induction (Figures 4H and S5H). Moreover, when we assumed that the initial AC position was distributed uniformly between P5.p and P6.p, the G-S model reproduced both the average P5.p and P7.p position and the broad distribution of P7.p positions relative to the AC that we observed experimentally (Figures 4G and 4H), with the main difference that the experimentally observed VPC positions were more widely distributed. Together, these results show that the experimentally observed variability in initial AC position alone could be sufficient to explain the intrinsic variability in VPC configuration observed in the lin-12(0) mutant. In addition, the overall agreement between the gradient-sensing model and the experiments suggests that it captures the essential features of LIN-3-induced VPC migration.

Notch Signaling Is Essential for Correction of Errors in AC Position and 1° Fate Induction

The results of the model showed that in principle LIN-3/EGF-induced VPC movement alone could correct for AC misplacement (Figures 5C, line 3, and 5D). However, this is only achieved for a narrow range of LIN-3/EGF level p0. Notch signaling, by restricting 1° fate to a single VPC, could drastically expand the range of LIN-3/EGF levels over which canalization would occur. Indeed, when we added Notch signaling to our mathematical model (Equations 7 and 8 in the STAR Methods), we found that now the correct expression pattern and position of the AC was realized for the full range of p0 (Figure 5E). To examine the contribution of Notch signaling to canalization, we sought to inhibit the Notch pathway without affecting AC specification, as the additional ACs in the lin-12(0) mutant likely increase the amount of secreted LIN-3/EGF. We shifted a dig-1 mutant carrying a temperature-sensitive lin-12 mutation, dig-1;lin-12(n676n930ts), to 25°C at the start of vulva induction. Indeed, these animals possessed a single AC (Figure 6B), even though lin-12(ts) animals often have multiple ACs when grown at 25°C during the L1-L2 larval stage (Sundaram and Greenwald, 1993). The model made a strong prediction for the effect of loss of Notch signaling: if LIN-3-induced VPC movement alone were sufficient to correct for AC misplacement and variability in 1° fate induction, corresponding to the model at low LIN-3 level, then dig-1;lin-12(ts) animals with inhibited Notch signaling would correct this as efficiently as dig-1 animals. However, for higher LIN-3 levels, the model predicted that inhibition of Notch signaling would result in multiple 1° fate cells with aberrant spatial configurations, similar to the lin-12(0) mutant (Figure 4).

Figure 6

Vulva Induction in Absence of Notch Signaling

(A) Relative AC position R as a function of time during vulva induction in dig-1;lin-12(ts) animals shifted to 25°C after AC specification. The difference with dig-1 animals in Figure 1E is significant (Kolmogorov-Smirnov test, p = 10−5).

(B) dig-1;lin-12(ts) animal with misplaced AC and Notch ligand expression in two VPCs. Expression of apx-1 (green) and lag-2 (red) is visualized using smFISH. Nuclei (blue) are stained with DAPI.

(C) Number of animals with relative apx-1 expression level E and relative AC position R for the late (9–12 hr) induction stage. The difference with dig-1 animals in Figure 3E is significant (Kolmogorov-Smirnov test, p < 10−5).

(D) Steady-state relation between relative AC position R and relative induction level E for the model with Notch signaling (green line) and without Notch signaling (purple). The relative induction level is given by E = (ϕ1 − ϕ2)/( ϕ1 + ϕ2), where ϕ1 and ϕ2 are the induction level in the closest and second-closest VPC. Both models were evaluated at LIN-3/EGF level p0 = 1.1 and with initial relative AC position R varying from equidistant to P5.p and P6.p (R = 1, indicated by the circle) to centered on P6.p (R = 0, square).

(E) Distance between the closest and second-closest VPC as function of the relative apx-1 expression level, both for dig-1 (green) and dig-1;lin-12(ts) animals (purple) at the late induction stage (9–12 hr). The dashed line corresponds to the average distance in the dig-1;lin-3(e1417) mutant without induction. Only animals in which the AC is positioned between P4.p and P5.p are shown.

(F) Distance between the closest and second-closest VPC as function of the relative induction level, both for the model with Notch signaling (green) and without (purple), for the same range of initial AC positions as in (D). The dashed line corresponds to the distance between VPCs in the absence of induction, p0 = 0.

We found that during early induction (0–3 hr, Figure 6A) the relative AC position in dig-1;lin-12(ts) at 25°C showed a wide distribution, similar to dig-1 animals (Figure 1E). Whereas dig-1 animals efficiently canalized the variability in AC position and 1° fate induction, we found, instead, that dig-1;lin-12(ts) animals at 25°C showed a wide range of relative AC positions and apx-1 expression patterns (Figures 6A–6C). Moreover, animals with a misplaced AC also typically showed equal expression of Notch ligands in the closest two VPCs (E ≈ 0, Figures 6B, 6C, and S3). At the same time, dig-1;lin-12(ts) animals at 15°C corrected AC misplacement and restricted Notch ligand expression almost as well as dig-1 animals (Figure S3). These observations differ from those in Grimbert et al. (2016) where correction of AC misplacement was not impacted by RNAi knockdown of lin-12. This could be due both to a weaker effect of RNAi compared to the lin-12(ts) mutant and to that it was performed in wild-type animals, where severe AC misplacement is infrequent compared to dig-1 mutants: we also observed a high fraction of animals with a single aligned 1° fate VPC at the end of induction in the lin-12(ts) single mutant at 25°C (Figure S7).

We found that the wide distribution of configurations with either one or two induced VPCs observed in dig-1;lin-12(ts) animals resembled the steady states generated by the model without Notch signaling at the transition regime at p0 ≈ 1.1 (Figure 5D). Specifically, the model reproduced the experimental observation that in configurations with two induced VPCs, the cells are equidistant to the AC (E ≈ 0, R ≈ 1, Figures 6C and 6D), whereas in configurations with a single induced VPC, the cell is aligned with the AC (E ≈ 0, R ≈ 0). The model also predicted a correlation between induction pattern and VPC deformation: as two induced VPCs will migrate toward another, the distance between them should be smaller than that of a single induced VPC to its neighbors (Figure 6F). Indeed, we confirmed experimentally that VPCs were positioned closer together if they exhibited more similar apx-1 expression levels (E ≈ 0) when compared to VPC positions both in dig-1 animals that have a single 1° fate VPC and lin-3(e1417) animals with no induced VPCs (Figure 6E). Taken together, these results show that for wild-type animals VPC movement alone is insufficient to produce the correct 1° fate pattern and that Notch signaling between VPCs is essential to restrict 1° fate induction and VPC migration to a single VPC.

VPC Movement Helps Notch Signaling Correct Errors in 1° Fate Induction

However, given that VPC movement alone can generate the correct 1° fate pattern for the full range of initial AC positions but a narrow range of low LIN-3/EGF levels, 0.7 < p0 < 1 (Figures 5C, line 3, and 5D), we examined whether VPC movement could support Notch signaling in restricting 1° fate to a single VPC for p0 = 1.1, the higher level consistent with the dig-1;lin-12(ts) mutant. Indeed, we found that adding a small amount of VPC movement to the model with Notch signaling dramatically increased the difference in steady-state induction level between the closest and second closest VPC (Figure 7A). Increasing the amount of VPC movement further did not impact the steady-state induction level, but did significantly increase the speed of 1° fate restriction (Figure 7B). While the model dynamics for weak Notch signaling and strong VPC movement was inconsistent with the experiments (Figures 3C–3E and 7D), for strong Notch signaling the model reproduced the experimentally observed dynamics also in the limit where VPC movement strongly contributes to the speed of 1° fate restriction (Figure 7C). These results show that, during 1° fate induction, VPC movement and Notch signaling might be highly intertwined: in this picture, Notch signaling is required for correction of AC misplacement, by restricting VPC movement to a single cell, but VPC movement in turn helps Notch signaling by bringing one cell closer to the AC and thereby increasing its level of Ras activation and Notch ligand expression. As a result, even though in the model VPC movement (Figure 5D) and Notch signaling (Figure 7A) alone can produce the correct 1° fate induction pattern, albeit for a limited range of parameters, the combination of both mechanisms canalizes variability in AC position and 1° fate induction much more efficiently and over a much wider range of biochemical parameters such as the LIN-3 level.

Figure 7

Cooperation between Notch Signaling and Cell Migration in Canalization of 1° Fate Induction

(A) Restriction of 1° fate induction to a single VPC for the model with Notch signaling (Equations 1, 2, 3, 4, 5, 6, and 7 in the STAR Methods), as function of the Notch inhibition strength KS and cell migration strength f0. Simulations were started with two VPCs approximately equidistant to the AC, R0 = 0.9. Color indicates the quality of restriction at steady state, given by ϕ2/ϕ1. Here, ϕ1 and ϕ2 are the induction level in the closest and second closest VPC and complete restriction of 1° fate to the closest VPC corresponds to ϕ2/ϕ1 = 0. The dashed line corresponds to a steady-state relative induction level E = 0.95, where E = (ϕ1 − ϕ2)/( ϕ1 + ϕ2). Red markers correspond to the simulations highlighted in (C) and (D).

(B) Time required to restrict 1° fate induction to a single VPC. Simulations were performed as in (A) and the time was measured until the relative induction level reached a value E = 0.95. The color indicates the time normalized to the longest time encountered for range of KS and f0 examined here. For the gray region below the dashed black line, E < 0.95 in steady state and hence the time was not recorded. Red markers correspond to the simulations highlighted in (C) and (D).

(C and D) Time dynamics of relative AC position R and relative induction level E for (C) strong and (D) weak Notch signaling and different levels of cell migration. Each trajectory corresponds to a different initial AC position ranging from R0 = 0.1 and R0 = 0.9. Note that in this deterministic model the relative induction level is always positive. Numbers correspond to the markers in (A) and (B).

Discussion

How development always results in the same adult structures despite strong genetic, environmental, and intrinsic variability is a fundamental unsolved question. In C. elegans vulva induction, the AC induces vulva cell fate in the VPCs in a distance-dependent manner, resulting in an invariant cell fate pattern with 1° fate only induced in the closest VPC, P6.p. Surprisingly, we found that the initial position of the AC with respect to the VPCs showed strong variability, with the AC located equidistant to two VPCs, typically P5.p and P6.p, in ∼20% of wild-type animals. A strong initial variability, correlated with AC position, was also observed in 1° fate induction, as measured by expression of the Notch ligand apx-1. Specifically, if the AC was equidistant to two VPCs, 1° fate was often initially induced at equal levels in both VPCs (Figure 3). However, we observed that vulva induction was canalized, i.e., both the variability in AC position and 1° fate induction decreased in time, resulting in the same configuration, with a single 1° fate cell aligned with the AC, in all animals (Figure 3).

We found that AC misplacement was corrected by movement of the 1° fate VPC that depended on LIN-3/EGF from the AC (Figures 1 and 4). Canalization of 1° fate induction was achieved by the progressive restriction of 1° fate to a single VPC that occurred simultaneously with correction of AC misplacement (Figure 3). A mathematical model of LIN-3-induced movement of VPCs toward the AC showed that VPC movement alone, by bringing one VPC close to the AC, could be sufficient to restrict 1° fate to a single VPC (Figure 5). However, we found that Notch signaling between VPCs was crucial for 1° fate restriction: upon inhibiting Notch signaling, most animals failed to canalize variability in AC position and 1° fate induction and instead exhibited a range of configurations with one or two 1° fate VPCs at varying positions to the AC (Figure 6). Yet, modeling showed that VPC movement can significantly improve 1° fate restriction by Notch signaling, leading to stronger and faster amplification of differences in 1° fate induction between adjacent VPCs (Figure 7). This shows that even though both VPC movement (Figure 5D) and Notch signaling (Figure 7) alone can generate the correct 1° fate induction pattern, both mechanism combined do so in a highly improved manner and for a much wider range of parameters.

The above mechanism acts by restricting 1° fate to the closest VPC rather than to P6.p specifically, raising the question how 1° fate is invariably induced in P6.p. We found that in wild-type animals, severely misplaced ACs can be shifted toward either P5.p or P7.p (Figure S2). Yet, wild-type animals under normal growth conditions only rarely exhibit a VPC fate pattern that is shifted so that P5.p or P7.p assume 1° fate instead of P6.p (Braendle and Félix, 2008, Grimbert and Braendle, 2014). As a possible solution to this paradox, however, we observed no wild-type animals where the AC was positioned closer to P5.p or P7.p than to P6.p (Figure S2). This suggests the AC might be positioned just accurately enough that 1° fate is restricted to P6.p in almost all animals, even when using a mechanism that purely selects the closest VPC.

The observed alignment of the AC and the 1° fate VPC is likely also important for many aspects of vulva cell fate patterning and morphogenesis that occur after 1° fate induction. First, the AC is responsible for induction of VulE and VulF fate in P6.p descendants by LIN-3/EGF signaling (Wang and Sternberg, 2000). Second, a Wnt signal originating from the AC controls the proximal-distal orientation of the asymmetric divisions of the 1° and 2° fate lineages, with the daughter cell closest to the AC assuming a different fate than the more distant daughter (Green et al., 2008). In both these cases, a misplaced AC could result in incorrect cell fate assignment. Finally, the AC is also instrumental in patterning the ventral uterus and connecting the uterus to the vulva (Newman et al., 1996, Sherwood and Sternberg, 2003), the success of which likely requires the precise relative alignment of the uterine cells, the AC and the cells of the P6.p lineage.

We have constructed the first model of vulva induction that takes into account VPC movement and deformation. The model not only generated the qualitatively different induction patterns observed in wild-type animals and mutants with different numbers of 1° fate VPCs (Figure 5) but also correctly reproduced the quantitative changes in VPC configurations observed in these mutants (Figures 4 and 6). In general, the model shows that LIN-3-induced movement is a potent force for establishing robust patterns, even in absence of Notch signaling. As such, it might have implications for vulva induction in other nematode species. For instance, in the Panagrolaimidae both P6.p and P7.p assume 1° fate in an AC-dependent manner and are positioned equidistant to the AC (Félix and Sternberg, 1997, Félix et al., 2000), a stable configuration naturally generated by our model (Figure 5C, line 2). Moreover, it shows that nematode species with the same 1° fate pattern as C. elegans might vary substantially in the relative importance of VPC movement and Notch signaling in restricting 1° fate to a single VPC. In general, the combination of long-range induction followed by Notch inhibition between and migration of the induced cells occurs more widely, for instance during development of the Drosophila tracheal system (Sutherland et al., 1996, Ikeya and Hayashi, 1999, Ghabrial and Krasnow, 2006) and blood vessel formation in vertebrates (Ochoa-Espinosa and Affolter, 2012). Our results suggest that these common mechanisms might provide robustness, in particular to intrinsic variability in the position of the source of the inductive signal relative to the induced cells.

Finally, our observation that Notch signaling is required to canalize variability in 1° fate induction, caused in turn by variation in AC position, represents a significantly different role for Notch signaling than has been assumed so far. Specifically, this result could clarify a long-standing debate on the role of Notch signaling in vulva induction, namely whether 1° fate is induced in multiple VPCs and Notch signaling is required to restrict 1° fate to a single VPC (the “graded” model) (Yoo et al., 2004), or whether 1° fate is induced exclusively in P6.p and Notch signaling is only required subsequently to induce 2° fate in P5.p and P7.p (the “sequential” model) (Simske and Kim, 1995). When the AC is initially correctly placed relative to P6.p, the LIN-3/EGF gradient might be narrow enough to induce 1° fate only in P6.p, even without Notch signaling. Indeed, we often observe 1° fate induction in a single VPC in the absence of Notch signaling in lin-12(ts) mutants (Figure S7), where the initial variability in AC position is reduced compared to dig-1;lin-12(ts) mutants. However, in animals where the AC is severely misplaced, Notch signaling is indeed essential to restrict 1° fate to a single VPC (Figure 6). Hence, this role of Notch signaling might exist as an insurance policy against relatively infrequent cases of AC misplacement. Thus, it suggests that some aspects of signaling networks make sense as adaptations only when viewed in light of the variability encountered during development.

STAR★Methods

Key Resources Table

Contact for Reagent and Resource Sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Jeroen van Zon (j.v.zon@amolf.nl).

Experimental Model and Subject Details

Strains

Wild-type nematodes were strain N2. The following mutants were used in this study: LGIII: lin-12(n941) (Greenwald et al., 1983), dig-1(n1321) (Thomas et al., 1990), unc-32(e189) lin-12(n676n930) (Sundaram and Greenwald, 1993), LGIV: lin-3(e1417) (Hwang and Sternberg, 2004), lin-1(n1790) (Jacobs et al., 1998). All strains were handled according to the standard protocol (Brenner, 1974). Briefly, animals were grown on agar plates containing Nematode Growth Medium (NGM) and E. coli strain OP50 as a food source. Unless indicated otherwise, all strains were grown at 20°C. To inhibit Notch signaling without producing additional ACs, dig-1(n1321);unc-32(e189) lin-12(n676n930) eggs were placed on NGM plates at the permissive temperature (15°C) and were allowed to hatch. After 24 to 34 hr of plating the eggs, larvae were shifted to the restrictive temperature (25°C). Animals with two ACs were excluded. We also excluded animals carrying the dig-1(n1321) allele that display a dorsal gonad. The unc-32(e189) allele, which on its own has no vulva development defects, is used to follow the lin-12(n676n930) mutation to which it is closely linked.

Method Details

Single Molecule Fluorescence In Situ Hybridization

To visualize mRNA transcripts, smFISH hybridization was performed as previously described (Raj et al., 2008, Ji and van Oudenaarden, 2012). Probes for smFISH were designed for optimal GC content using a web-based program (http://singlemoleculefish.com) and were coupled to Cy5 (GE Amersham) or Alexa594 (Invitrogen). The sequences of the oligonucleotide probes used in this study are listed in Table S1. Animals were collected by washing plates with M9 and were fixed in 4% formaldehyde in 1XPBS for 45 min at room temperature. Fixed animals were permeabilized in 70% ethanol overnight at 4°C. Subsequently, animals were incubated with the smFISH probes overnight at 30°C in hybridization solution containing 10% formamide. The next day, animals were washed twice with 10% formamide and 2X SSC, each time followed by an incubation for 30 min at 30°C. To visualize cell nuclei, DAPI was added at 5 μg mL-1 at the last wash step. Microscopy images were acquired with a Nikon Ti-E inverted fluorescence microscope, equipped with a 100X plan-apochromat oil-immersion objective and an Andor Ikon-M CCD camera controlled by μManager software (Edelstein et al., 2014). Exact three-dimensional positions of smFISH fluorescent spots in each animal were detected using a custom MATLAB (The Mathworks) script, based on a previously published algorithm (Raj et al., 2008). Briefly, we first convolved smFISH microscopy images with a Gaussian filter to increase the brightness of spots of the correct size and suppress the background signal. Next, we select candidate spots by thresholding, using a manually determined threshold. We further refined the candidate spots by finding regional intensity maxima within each candidate spot, to separate smFISH spots whose fluorescence signals are partially overlapping. Finally, the resulting smFISH spots were manually assigned to individual VPCs. We converted gonad length , measured as the distance between the two DTCs, to time using the expression for and for , where is in units of , in units of hours and corresponds to the start of the L2 larval stage (van Zon et al., 2015). Distances between DTCs as well as between the AC and VPCs were all measured along the body axis of the animal.

Experimental Design

In general, we did not perform a systematic replication of our experimental results. However, two key datasets, on Notch ligand expression in wild-type animals and dig-1 mutants, combine data gathered by different researchers on different microscopy set-ups which showed excellent agreement. During data analysis, no randomization and blinding strategies were followed. Exclusion of animals only occurred of mutant animals with gross morphological deviations that rendered them irrelevant to our studies, such as dig-1 mutants with a dorsally located gonad.

Mathematical Model of Vulva Induction

Overview

In the model, the position and shape of VPCs is determined by the position of basal vertices and apical vertices (Figure S5A). In absence of migration, this is implemented by assuming that movement of vertices minimizes an energy function of the form (Fletcher et al., 2014):where is the surface area of cell , is the cell elasticity, is the length of edge and its line tension. For simplicity, we ignore cell growth and hence each cell has a constant preferred cell area . This leads to the following equation of motion for the vertex model that includes migration:where we assumed the movement of vertices is overdamped with the mobility, is the elastic cell area force, is the line tension force and is an additional force that is not derived from the energy function in Equation 1 and describes the cell migration in response to the LIN-3/EGF gradient. In the attractant-maximization (A-M) model, cells seek to maximize the amount of attractant integrated over the cell surface (Savill and Hogeweg, 1997) leading to:

In the gradient-sensing (G-S) model, cells generate internal polarity based on the difference in chemoattractant concentration measured over the cell surface (Szabó et al., 2010). Here, for simplicity, we assume a linear relationship between polarity and chemoattractant gradient, leading to:

For both the A-M and G-S model, is a coupling constant that describes how strongly cell induces cell migration in response to the LIN-3/EGF signal and is the the local LIN-3/EGF concentration given by:where is the total amount of LIN-3/EGF excreted by the AC, is the decay length of the LIN-3/EGF gradient and is the AC position. Furthermore, for both models we assume that the LIN-3/EGF gradient is sensed and migration forces are produced at the basal VPC surface. We can calculate , the migration force produced by cell , as . For both the M-A and G-S model, this yields the same expression:where and are the anterior and posterior basal vertex of cell . Based on our experimental results in Figures 1 and 4, we assume that depends on the Ras activation level as , where is the VPC migration strength. Because it is not known how LIN-3/EGF or the Ras pathway controls cell migration, we assume that Ras activation depends cooperatively on the external LIN-3/EGF signal through a standard Hill curve:where is the dissociation constant, the Hill coefficient and is the center of the basal surface of the cell. Notch signaling from neighboring VPCs inhibits LIN-3/EGF induction (Yoo et al., 2004, Berset et al., 2001) and this is implemented in Equation 6 by the term , with indicating the amount of Notch signal received by cell and corresponding to a model without Notch signaling. While expression of Notch ligands is induced by LIN-3/EGF signaling, all VPCs express Notch receptors, which can be activated by Notch ligands expressed in neighboring VPCs (Sternberg, 2005). Therefore, we assume that the activation of Notch signaling in cell depends on the level of LIN-3/EGF induction in the neighboring VPCs and :where is the rate of spontaneous Notch signal deactivation. In this picture, corresponds to the fraction of activated Notch receptors and to the Notch signaling strength, i.e., the amount of inhibition of Ras activation per activated Notch receptor. The full equations of motion of the A-M and G-S models as well as values of all parameters are discussed next.

Vertex Model of Vulva Precursor Cell Shape without Migration

In so-called “vertex models,” cells are represented by polygons with edges corresponding to the cell membrane (Fletcher et al., 2014). In such models, cell deformation and cell movement are described entirely by movement of vertices. Following the standard approach to modeling epithelial sheets, we reduce the complexity of the model by treating the cells as two-dimensional entities, thereby considering only dynamics along the antero-posterior (A-P) and apical-basal (Ap-Ba) axis (Figure S5A). In this picture, the shape and position of vulva precursor cells are determined by vertices , where and correspond to the basal and apical vertices, respectively. In this description, the center of mass of cell is given by:

To derive the elastic cell area force from the energy function in Equation 1, we calculated the cell area from the positions of the vertices that constitute the cell, as follows:where . Using this, we find:where . To derive the line tension force , we used the following expression for the total line tension calculated:where and are the line tension in the Ap-Ba and A-P direction, respectively. Using this, we arrive at:

To implement constrains in deformation due to the presence of the gonad on the basal side of vulva precursor cells, we assume that basal vertices, , are constrained to move along the A-P axis, corresponding to x axis in our model:where is the unit vector in the x axis. Finally, we assume that the vulva precursor cells are attached to fixed outer boundaries, reflecting the constraints to movement imposed by the rest of the animal’s body. This is implemented by imposing that the vertices at the outer A-P boundaries do not move:

Vertex Model of Vulva Precursor Cell Migration by Attractant-Maximization

In tissue modeling approaches that rely on energy functions, such as vertex model and cellular Potts models, cell migration in a gradient of chemoattractant is often implemented by adding a term to the energy function in Equation 1 that decreases the energy with increasing concentration of the chemoattractant:where the energy is proportional to the chemoattractant concentration integrated over the basal surface of the vulva precursor cell and is a coupling constant that indicates how strongly each individual cell responds to the chemoattractant signal. In such a description, cells change their shape and position to maximize the chemoattractant concentration integrated over their cell surface or volume, causing cells to move in the direction of increasing chemoattractant concentration (Savill and Hogeweg, 1997). The migration force on vertex is then given by . For simplicity, we assume that the chemoattractant concentration only depends on the distance along the A-P axis, corresponding here to the x axis. In addition, we assume that migration forces are only produced on the basal surface, corresponding to the vertices . This results in the following expression for the migration force :which reduces to the expression for the migration force in Equation 3. The total migration force produced by a cell is defined as , giving rise to Equation 6. Equation 3 shows that for the attractant-maximization mechanism, movement of vertices at the interface between cells and requires that . In practice, this means that in the presence of fixed outer boundaries (Equation 17), movement and deformation of cells only occurs if varies between cells. In addition, Equation 3 can be interpreted as the sum of two forces acting on vertex , one produced by cell and one by cell . Comparing the sign of the forces shows that they always point outward from the cell. As a consequence, if cell responds more strongly to the chemoattractant than its neighbors, , it will by definition expand along the direction of the gradient while contracting in the orthogonal direction in order to maximize the overlap between the cell surface and the attractant.

Vertex Model of Vulva Precursor Migration by Gradient-Sensing

An alternative approach to modeling migration in a chemoattractant gradient is gradient-sensing (Szabó et al., 2010). In this picture, cells measure the external attractant gradient over the surface of the cell, e.g., by comparing the attractant concentration at the front and the back of the cell. Such measurements then lead to an internal polarization of signaling proteins, ultimately giving rise to migration toward high chemoattractant concentration. Here, we do not explicitly take into account the internal polarization machinery, but make the simplest assumption that the migratory force is proportional to the change in concentration, integrated over the basal surface of the cell, and pointing in the direction of steepest increase in attractant (Szabó et al., 2010), resulting in Equation 4. The expression for the total migration force produced by a cell , defined as is identical to that of the attractant-maximization model (Equation 6). As in the attractant-maximization model, Equation 4 shows that the force on vertex can again be interpreted as the sum of two migration forces, produced by cells and . In contrast to the attractant-maximization model, Equation 4 shows that the migration force produced by each cell has the same direction and magnitude for each vertex in that cell. As a consequence, cells do not expand along the direction of the gradient in response to the chemoattractant.

Expressions for Chemoattractant Gradient and Coupling Constant

For both the attractant-maximization and gradient-sensing models, we assume that the chemoattractant forms an exponential concentration gradient that depends only on the distance to the anchor cell along the A-P axis (Equation 5). For the attractant-maximization model, we already showed that needs to vary between cells to result in VPC movement, even for a non-uniform attractant profile . In the gradient-sensing model, movement also arises for constant when is non-uniform. However, to reproduce our experimental data we found that for both models we have to make dependent on the distance to the AC (Figure S6). Our experiments show that the propensity to migrate depends on the presence of LIN-3/EGF (Figures 1F and 1I). Moreover, we find a correlation between VPC displacement and Notch ligand expression levels (Figures 4A–4C), with Notch ligand expression itself being a measure of the level of LIN-3/EGF induction (van Zon et al., 2015). The simplest picture consistent with these observations is one in which corresponds to the LIN-3/EGF gradient and is proportional to amount of Ras signaling in cell at position . Specifically, we assume that , where is a parameter that describes how much force is produced for a given level of Ras activation . For simplicity, we assume that depends on the external LIN-3/EGF concentration though a Hill curve (Equation 7). We find similar dynamics when is proportional to integral of over the basal surface of the cell (data not shown). We incorporated Notch signaling by adding an inhibitory term to Equation 7 that corresponds to the strength of the Notch signal received by cell .

Parameters and Simulation

We simulate N = 9 cells, which would in our picture correspond P(2-10).p, with only minor deformations seen in the most distal VPCs for parameters used here. To accurately capture the shape of VPCs, we estimated the aspect ratio of the cell body in the A-P and Ap-Ba direction from the limited experimental data available (Whitfield et al., 1999, Skorobogata et al., 2014, Haag et al., 2014) and found a value of , meaning the cell dimension in the A-P direction was times larger than the dimension in the Ap-Ba direction. We fixed the preferred VPC area to . Then, , the undeformed VPC length in the A-P direction, is given by . Furthermore, for the line tensions and of edges in the Ap-Ba and A-P direction respectively, we assume that . This assures that the VPCs have the correct aspect ratio without being stretched by their connection to the fixed outer boundary: if this connection would be severed, i.e., Equation 17 would not hold, VPCs would maintain the correct aspect ratio in absence of a chemoattractant cue. The remaining mechanical parameters, and , were tuned to reproduce the qualitative VPC configurations in Figure 4C and the shape of the diagram in Figure 4D. In particular, we choose so that the deviation in cell area from was less than 3% through the simulation. Otherwise, varying these parameters most strongly impacted the magnitude of the deformations upon induction.

The LIN-3/EGF level and the decay length play analogous roles, i.e., increasing either one of them increases the range of induction. We choose here to vary since the main difference between wild-type animals and the lin-3 and lin-12 mutants studied is likely the LIN-3/EGF level rather than the decay length. We choose a decay length comparable to the length of the VPC body along the A-P axis, which is for our choice of parameters. The parameters , and were chosen to produce a force-distance curve capable of generating the configuration of two induced VPCs equidistant to the AC (Figures 5 and S6). For that, we have to assume a high degree of cooperativity, . However, the resulting induction profile is consistent with that measured experimentally using Notch ligand expression as a marker of induction by LIN-3/EGF (van Zon et al., 2015). Here, for simplicity we use a Hill curve to generate this induction profile. However, qualitatively similar induction curve can be generated with lower degrees of cooperativity, e.g., by taking into account saturation of EGF receptors (van Zon et al., 2015) or different shapes of the LIN-3/EGF gradient. The parameter determines how much force is produced for a given amount of EGF/Ras activation and magnitude of LIN-3/EGF difference measured over the VPC surface and, together with the line tensions and , the main parameters controlling the amount of VPC deformation. The Notch signaling parameter and were chosen so that the model dynamics reproduced that observed experimentally (Figures 3C–3F).

At the start of the simulation, the undeformed VPCs are defined by the vertices . We then calculate the AC position so that , where is the relative AC position and the AC is positioned between VPCs and . We solve the differential equations using the scipy.integrate.odeint function in Python. For the parameters used here, the dynamics of deformation is sufficiently slow and the deformations themselves sufficiently small that no topological transitions occur during the simulations.

We generally use the same parameters for the A-M and G-S models. However, due to the difference in cell shape dynamics between the two models, some parameters had to be chosen differently to reproduce the experimental data. Below, we give an overview of all parameter values used, with those parameters in the A-M model that differ from the G-S model indicated by an asterisk.

Mechanical parameters

EGF/Ras induction parameters

Notch signaling parameters

Quantification and Statistical Analysis

TwoDimensional KolmogorovSmirnov Test by Bootstrapping

To compare the two-dimensional distributions in Figures 1D–1F, 2A, 2C–2E, 6A, and 6C between different strains, we calculated P values as follows: given two samples and , with respective sizes and , and whose elements consist of points defined in a two-dimensional space, e.g., relative AC position versus relative expression level, we want to test the null hypothesis that both samples were drawn from the same unknown distribution. To do this, we used the two-dimensional Kolmogorov-Smirnov (K-S) test combined with bootstrapping strategy, following Press et al. (1996). We will briefly outline our approach below.

The K-S test relies on a statistic that serves as a measure of the difference between two measured distributions and . is generally defined as the maximum absolute difference between their two cumulative distributions. In one dimension, the cumulative distribution of a sample is defined as , which represents the probability that a variable drawn from the sample will have a value that is less or equal than . In this case, . In two dimensions, there are four ways of defining the cumulative distribution of a sample, namely , , , and , where is a random variable drawn from the sample. In this case, is calculated as:

Practically, we compute for the two samples and as follows. For each data point in , we calculate . We then calculate the maximum difference between the two cumulative distributions over all data points, . Similarly, we calculate for each data point and subsequently the maximum value over all data point , . The two-dimensional K-S statistic is now given by .

To obtain the significance level or p value we use the following bootstrapping strategy. We concatenated both samples to obtain a pooled dataset of length , i.e., We then resample with replacement times, obtaining a bootstrapped sample at iteration . Next, is divided into two samples and , and we compute the respective K-S statistic , as outlined above. Then, the significance level is given by the fraction of bootstrapping iterations in which the resulting statistic exceeds the statistic from the original samples. That is,where is the Heaviside step function. For all reported P values we used bootstrap samples. The p values reported for the comparison between Figures 3C–3E and 6C in the main text are calculated for the underlying distributions shown in Figure S3.

Author Contributions

G.H.-P. and J.S.v.Z. designed the experiments. G.H.-P. conducted the experiments. G.H.-P. and J.S.v.Z. constructed the mathematical model and performed simulations. G.H.-P. and J.S.v.Z. wrote the manuscript.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, Peptides, and Recombinant Proteins
lag-2 smFISH probes	Biosearch technologies (van Zon et al., 2015)	N/A
apx-1 smFISH probes	Biosearch technologies (van Zon et al., 2015)	N/A
Dextran sulfate	Sigma-Aldrich	D8906-5G
E. coli tRNA	Sigma-Aldrich	R1753-500UN
Vanadyl ribonucleoside complex	New England Biolabs	S1402S
Formamide	ThermoFisher Scientific	AM9342
Nuclease free water	ThermoFisher Scientific	AM9932
Nuclease free 20x SSC	ThermoFisher Scientific	AM9770
Nuclease free 10x PBS	ThermoFisher Scientific	AM9624
RNase free BSA	ThermoFisher Scientific	AM2616
Tris-HCl, pH 8.0	Sigma-Aldrich	T2694-100ML
Catalase	Sigma-Aldrich	C3515-25MG
Glucose Oxidase from Aspergillus niger	Sigma-Aldrich	G2133-10KU
TE buffer	Sigma-Aldrich	93283-100ML
Experimental Models: Organisms/Strains
lin-12(n941)	Prof. R. Horvitz (Greenwald et al., 1983)	MT16472
dig-1(n1321)	Caenorhabditis Genetics Center (Thomas et al., 1990)	MT2840
unc-32(e189) lin-12(n676n930)	Caenorhabditis Genetics Center (Sundaram and Greenwald, 1993)	GS60
lin-3(e1417)	Caenorhabditis Genetics Center (Hwang and Sternberg, 2004)	CB1417
lin-1(n1790)	Caenorhabditis Genetics Center (Jacobs et al., 1998)	WU125
Software and Algorithms
smFISH Analysis Software	(van Zon et al., 2015)	N/A

Parameter	Definition	Values
N	Number of cells	9
ρ	VPC aspect ratio	10
K	Area elastic modulus	1
A^0	Preferred VPC area	1
γ_∥	Line tension for A-P edges	0.01
μ	Mobility coefficient	1

Parameter	Definition	Values
p_0	Total excreted LIN-3 level	0−1.6, 0−10∗
λ	LIN-3 gradient decay length	2.5, 1.9∗
K_p	Activation coefficient	0.1
n	Hill coefficient	5
f_0	Migration strength	0.15, 0.05∗

Parameter	Definition	Values
K_S	Notch signaling strength	10^−2
b_S	Notch signal deactivation rate	10^−3