# Evolution of developmental bias explains divergent patterns of phenotypic evolution in two nematode clades

Joao Picao-Osorio, Charlotte Bouleau, Pablo M Gonzalez de la Rosa, Lewis Stevens, Nina Fekonja, Mark Blaxter, Christian Braendle, Marie-Anne Félix

► To cite this version:

Joao Picao-Osorio, Charlotte Bouleau, Pablo M Gonzalez de la Rosa, Lewis Stevens, Nina Fekonja, et al.. Evolution of developmental bias explains divergent patterns of phenotypic evolution in two nematode clades. Proceedings of the National Academy of Sciences of the United States of America, 2025, 122 (34), 10.1101/2025.02.04.636455v1. hal-05336432

HAL Id: hal-05336432

https://hal.science/hal-05336432v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Evolution of developmental bias explains divergent patterns of phenotypic evolution
in two nematode clades

Joao Picao-Osorio¹†*, Charlotte Bouleau², Pablo M. Gonzalez de la Rosa³, Lewis Stevens³, Nina Fekonja¹, Mark Blaxter³, Christian Braendle²* and Marie-Anne Félix¹*

¹ Institut de Biologie de l'École Normale Supérieure, CNRS, Inserm, ENS, PSL; 46 rue d'Ulm, 75005 Paris, France.

² Institut de Biologie Valrose, Université Côte d'Azur, CNRS, Inserm; 06000 Nice, France.

³ Tree of Life, Wellcome Sanger Institute; Cambridge, UK.

*Corresponding authors. jgosorio@ciencias.ulisboa.pt; christian.braendle@univ-cotedazur.fr; felix@bio.ens.psl.eu

† Present address: Centre for Ecology, Evolution and Environmental Changes, CHANGE, Department of Animal Biology, Faculdade de Ciências, Universidade de Lisboa; 1749-016 Lisboa, Portugal.

Preprint: https://www.biorxiv.org/content/10.1101/2025.02.04.636455v1

Classification: Biological Sciences - Evolution

Keywords: evolution, developmental bias, mutational variance, *C. elegans*

Short title: Evolution of developmental bias
34

**Abstract**

35

36 Rates of phenotypic evolution vary across traits, and these evolutionary patterns themselves 37 evolve. Understanding how development contributes to such patterns remains a challenge because 38 it requires large-scale measurement of phenotypic variation resulting from random mutations 39 across multiple species. Using the experimentally tractable system of nematode vulval 40 development, we quantified the evolution and mutational sensitivity of six cell fates across two 41 clades. First, we find that the cell whose fate evolves fastest is distinct for each clade: P3.p in 42 *Caenorhabditis*, P4.p in *Oscheius*. Second, we show that this evolutionary pattern matches the 43 pattern of variation within species of each clade. Third, we find that, within each clade, the 44 differential effect on each cell of random mutation, quantified as their mutational phenotypic 45 variance, is sufficient to explain the observed phenotypic evolution. Finally, we propose that the 46 difference between clades is explained by a simple spatial shift across cells of the sensitive region 47 in a Wnt dose-response curve. These findings underscore the importance of integrating the 48 evolving structure of the genotype-phenotype map and measures of developmental sensitivity into 49 the understanding of phenotypic change at both micro- and macroevolutionary scales. 50

51

52

53

54

**Significance statement**

55

56 Phenotypic traits, such as vertebrate teeth size, evolve at different rates, and these rates may differ 57 among groups of organisms. For example, molar size may change in one vertebrate clade while 58 canine size changes in another. These evolutionary patterns are shaped by selection, but also by 59 how readily the underlying developmental system can be modified. Such developmental bias can 60 be probed by assaying the phenotypic effect of random mutations. Here, we use nematode vulval 61 cell fates to measure developmental bias in two clades in which a different cell undergoes fast 62 evolution in fate. We find that in each clade, this cell fate is also most sensitive to mutation. 63 Evolution of developmental bias can thus explain divergent patterns of phenotypic evolution.

63
64 Main Text

65
66 Phenotypic traits evolve at different rates, but it remains under active debate how different
67
68 mechanisms—such as natural selection, genetic drift, mutation and development—explain these
69 rates (1, 2). As noted by the palaeontologist G.G. Simpson, a prominent figure of the mid-20th
70
different clades of organisms may display distinct directions of evolutionary change (3). Natural
71
selection is a central force in the evolution of phenotypes, yet it only operates on available
72
phenotypic variation. Simpson remarked, that in the space of phenotypes, "It is not only
73
improbable but also inconceivable that mutations in every imaginable direction occur with equal
74
frequency," and called "to determine whether or not mutations (in the broadest sense) are
75
directional and parallel to the actual direction of evolution" (3). Simpson's question translates to
76
quantifying the spectrum of phenotypes produced by random mutation of the genotype: in the
77
framework of quantitative genetics, the mutational variance (V<sub>M</sub>) (4). Given that development
78
mediates the relationship between genotype and phenotype (5-7), the anisotropic pattern of
79
phenotypic variation generated by random mutation is called developmental bias (8-10).
80
Determining the importance of developmental bias in phenotypic evolution has long been
81
recognized as the central focus of evolutionary developmental biology (5-7, 10-15). In line with
82
the view that developmental bias plays a central role in shaping phenotypic variation and evolution,
83
David Houle's seminal study (16) demonstrated a strong positive relationship between mutational
84
variance and within-species heritable phenotypic variation (standing genetic variance, V<sub>G</sub>) in
85
*Drosophila melanogaster*'s wing dimensions, bristle number and life history traits. Furthermore,
86
subsequent studies on dipteran wing shape (17), nematode mitotic spindle traits (18) and vulval
87
precursor cells (19) showed that mutational variance is a sufficient predictor of both genetic
88
variance within a species and divergence among species.

89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
Developmental bias itself may evolve. Specifically, as Simpson suggested, developmental
bias may shift over macroevolutionary timescales. Here, we use the homologous cellular
framework of nematode vulval development to test whether evolutionary change in developmental
biases, as measured by the effect of random mutation, can explain the differences in phenotypic
evolution between two clades. In a wide range of nematode species, 12 ventral epidermal precursor
cells (called Pn.p cells, n=1 to 12) are positioned along the anterior-posterior axis of early larvae
(Fig. 1A). Of those, the six cells numbered P3.p to P8.p form the vulva competence group in
*Caenorhabditis elegans* (20). These six cells undergo a series of fate specifications and divisions
under the influence of several signalling pathways, including Wnt, EGF/Ras, and Notch. The
central cells P5.p, P6.p, P7.p normally form the vulva proper, required for oviposition and
copulation with males (Fig. 1B, Fig. S1A), but can be replaced by P3.p, P4.p and P8.p. These latter
three cells usually divide once at the third larval stage and fuse to the epidermal syncytium.
Remarkably, in *C. elegans* cell fates of P4.p to P8.p are robust to environmental and genetic
perturbations (19, 21, 22), and remain largely unchanged across the *Caenorhabditis* genus (23-
25). In contrast, P3.p has a stochastically variable cell fate in *C. elegans*: P3.p can either divide
once like P4.p or fuse directly with the syncytium during the L2 larval stage without division (21,
26). The frequency of this binary cell fate decision (L2 stage fusion versus division) evolves rapidly within and among *Caenorhabditis* species, and is highly sensitive to environmental and genetic perturbations (19, 22, 23, 25).

In other genera of the family Rhabditidae, such as *Oscheius* (24), the same six cells remain unfused during early larval development (27) and the vulva is also formed from P5.p, P6.p and P7.p progeny. However, in *Oscheius*, P3.p almost never divides and loses competence (28). In contrast, P4.p and P8.p in most *Oscheius* species can adopt one of three distinct cell division patterns (no, one, or two cell divisions), with the frequency of these patterns varying significantly both within and among *Oscheius* species and in response to environmental variation (23).

These observations support a compelling and testable hypothesis: divergent patterns of phenotypic evolution arise as a consequence of the evolution of developmental biases (Fig. S1B). Here we test this hypothesis using the nematode vulva precursor cell fates. The overall experimental design is shown in Fig. 1D, assaying variation in phenotypes at three levels, among species, within species and in response to random mutation. Altogether, we report here on the phenotype of over 80,000 individual animals. We use a quantitative genetic framework to analyze the data. We represent results on the scale of phenotypes in Fig. 2 and of estimated variance in Fig. 3. For readers unfamiliar with the quantitative genetic framework, consider these variances as measures of phenotypic variation for each trait (the fate of each cell) due to genetic variation at three evolutionary scales: in response to random mutation (M), and in intraspecific (G) and interspecific (R) variation (Fig. 1D).

## Divergent rates of cell fate evolution across clades

We started by determining the fates of the six vulval precursor cells using Nomarski optics in 53 *Caenorhabditis* and 24 *Oscheius* species (Fig. 1). To infer the phylogenetic relationships of these species (Figs S3, S4), we used 560 orthologues obtained from assembled genomes and transcriptomes. We calibrated the phylogenetic tree to estimate the divergence time between the *Caenorhabditis* and *Oscheius* genera at approximately 160 million years ago (MYA; Fig. 1C, Fig. S5). We then mapped the vulval cell fate changes onto this phylogeny (Fig. S6). The division frequency of P3.p showed the greatest variation in *Caenorhabditis*, whereas in *Oscheius*, P4.p exhibited the most variation, producing between one and four cells that all fuse to the epidermal syncytium ('S' fate, Figs. 1C, 2A, Figs S5, S6, S7). For simplicity, we consider for further data analysis only the most common pattern in *Oscheius*, where two divisions produce four cells, versus all alternative patterns. We used phylogenetic mixed models (PMM) (29) to quantify rates of phenotypic evolution by estimating the 'rate among species genetic variance per million year' ($V_R$). We thus find that, in *Caenorhabditis*, the fate of P3.p evolved approximately a thousand times faster than the fate of P4.p (Fig. 3A, 4A right panel). Conversely, in *Oscheius*, the fate of P4.p evolved about a hundred times faster than that of P3.p (Fig. 3A, 4A right panel). The phylogenetic relationships within each clade explained a large proportion of the phenotypic variance of the fastest-evolving cells (Fig. S8), as measured by phylogenetic heritability (30). This ensures that the source of the corresponding variation is not environmental. Thus, we conclude that the fates of
P3.p and P4.p cells display highly contrasting patterns of phenotypic evolution between *Caenorhabditis* and *Oscheius* species.

**Alignment of intraspecific variation with evolution in each genus**

Next, we asked whether we could detect corresponding patterns within species. For this, we scored in parallel wild isolates of *C. elegans*, *C. briggsae*, *O. tipulae* and *O. onirici* to estimate the genetic (*V*<sub>G</sub>) variances of the six traits. These four species were chosen because they reproduce primarily by selfing and, therefore, each strain is composed of isogenic individuals. Phenotyping followed a block design and blinded approach with four observers (see Methods).

Similar to results for interspecific evolution, P3.p within each *Caenorhabditis* species exhibited much greater phenotypic variation in division frequency compared to any other cell (Fig. 2B; Fig. S9). Specifically, *C. briggsae* and *C. elegans* displayed approximately ten and one hundred times higher *V*<sub>G</sub> for P3.p than for P4.p, respectively (Figs. 3B, 4A middle panel, Figs S10A, S11A). In contrast, we found that within the two *Oscheius* species, P4.p showed the greatest phenotypic variation (Fig. 2B, Fig. S9). In *O. tipulae* and *O. onirici*, *V*<sub>G</sub> of P4.p was a hundred times greater than that of P3.p (Figs. 3B, 4A middle panel, Figs S10B, S11B). Although three of the four species hardly show any variation for one of the phenotypic axes, *C. briggsae* exhibits a more complex pattern, displaying variation along both axes, yet with a stronger effect for P3.p fate. In conclusion, the rate of cell fate evolution differs significantly between the two genera and the axis of variation within species aligns with that observed among species of the genus.

**Evolutionary change in developmental bias**

We next asked whether the contrasting phenotypic evolutionary rates of the six precursor cells could be explained by the evolution of developmental bias between *Caenorhabditis* and *Oscheius*. To this end, we created eight panels of lines with *de novo* mutations induced by chemical mutagenesis in wild isolates of both genera, allowing us to quantify the mutational variance (*V*<sub>M</sub>) of cell fates (Fig. S2A). To account for genotype-specific effects on mutational variances, we performed the chemical mutagenesis on two different wild isolates for each species and two species per clade (Fig. 1D bottom panel). To account for environmental and other non-genetic effects, we phenotyped the ancestral isogenic wild isolates multiple times in parallel with the mutant lines (Datasets S9, S10). The phenotyping was also performed in parallel with the wild isolates scored to estimate the standing genetic variances (*V*<sub>G</sub>) mentioned above.

We found that, in each genus, vulval phenotypic variation generated by *de novo* mutation is anisotropic in phenotype space (Fig. 2C, Fig. S12), thus revealing a developmental bias. Most importantly, this developmental bias strongly differs between the two genera. In both *Caenorhabditis* species, the principal axis of phenotypic variation resulting from random mutation concerns P3.p fate: *C. briggsae* and *C. elegans* showed over ten times greater variation upon random mutation (*V*<sub>M</sub>) for P3.p compared to P4.p (Figs. 3C, 4A left panel), in line with our previous study using spontaneous mutation accumulation lines (19). The difference between the two cells is more nuanced in *C. briggsae* than in *C. elegans*, indicating evolution within the genus. In contrast, the two *Oscheius* species showed significantly greater phenotypic variation in P4.p fate
in response to new mutations (Fig. 2C, Fig. S12): P4.p presented a $V_M$ that was a thousand times greater than that of P3.p in both *O. tipulae* and *O. onirici* (Fig. 3C, 4A left panel).

The variances in fates of the two cells at the three levels (mutation, polymorphism, divergence) can be expressed as a ratio, which is robust to differences in mutation efficacy, species genetic diversity and genus age. These results demonstrate evolutionary change in developmental bias between the *Caenorhabditis* and *Oscheius* clades (Fig. 4A). Altogether, our results reveal a strong alignment between the axes of heritable phenotypic variation resulting from *de novo* mutations and the axes of natural variation within each species and genus, though they contrast between genera. Thus, the change in developmental bias matches the axis of phenotypic evolution within each genus (Fig. 4B).

## Discussion and Conclusions

The mechanisms underlying the tempo and mode of evolution have been a subject of active debate for a long time, aiming to bridge micro- and macro-evolution (3, 31-34). Our results imply that the evolutionary change in mutational variances between two nematode clades can explain patterns of phenotypic evolution of the vulva precursor cell fates. This has important implications for the role of development in shaping the tempo and mode of phenotypic evolution. Over short timescales within a lineage, species with a similar developmental bias are likely to exhibit shared patterns of phenotypic evolution (18, 35-37). Over longer timescales, lineages with distinct developmental biases may exhibit divergent patterns of phenotypic evolution across different traits. These observations are consistent with a mutational constraint scenario (38), in which developmental bias within a given clade shapes evolutionary divergence along the 'genetic lines of least resistance' (12, 13). Alternatively, theoretical studies indicate that natural selection can shape the structure of the mutational (co)variances (39, 40), thereby affecting the stability and evolution of standing genetic (co)variances (41, 42). This process could result in an alignment of mutational, genetic and divergence variances that reflects the evolutionary history of the population, species, or lineage. Thus, our results, together with those of others (Houle et al. 2017), on the triple alignment of variances within a clade may reflect mutational constraints influenced by the selection surface.

Our study was not aimed to address evolutionary mechanisms (selection or drift) underlying the macro-evolutionary change in bias. Regarding selection on these cell fates, the vulva is required for production of a large number of progeny, while variation in P3.p and P4.p fates may have little impact on organismal biology because these cells (or their progeny) will fuse to the epidermal syncytium regardless of whether they divide. We consider this scenario of apparently low selection on peripheral cells of the competence group to increase the potential for developmental bias to manifest.

The division of P3.p, P4.p and P8.p may be selectively maintained if they are used as a backup system in environments that increase deviant fates for the central vulval cells (22). Indeed, for cells to be competent, they must remain unfused at the L2 larval stage. Selection pressure for the maintenance of a large competence group appears weak, since some nematodes restrict competence to only those cells that actually form the vulva (24, 43, 44). Hence, the observed
227 evolution of mutational and standing genetic variances in these cell fate traits likely results from
228 selection on other traits, particularly those involving the same signalling pathways, such as the
229 highly pleiotropic Wnt pathway (45).

230 Mechanistically, we propose that the evolution of developmental bias in nematode vulva
231 cell fate patterning is due to distinct developmental sensitivities of the cells, based on their anterior-
232 posterior positioning in a Wnt gradient (Fig. 4C). Past research shows that in *C. elegans*, altering
233 the doses of either of two genes encoding Wnt ligands, *egl-20* and *cwn-1*, significantly affects the
234 P3.p division frequency more than that of P4.p (46, 47) (Fig. 4C, D). Specifically, reducing the
235 Wnt dose leads to a lower P3.p division frequency. The two Wnt genes, especially *egl-20*, are
236 transcribed in a posterior body region, forming a protein gradient that extends toward the anterior
237 (48, 49), placing P3.p at the fading end of this gradient. Indeed, interindividual variability in
238 downstream Wnt signalling through β-catenin/BAR-1 can be detected in P3.p (but not in P4.p),
239 and is linked to its stochastic cell fate outcome (50). These observations strongly suggest that the
240 anterior position of P3.p contributes to its heightened developmental sensitivity to stochastic,
241 environmental, and genetic variation (22). In line with this, we previously showed that the high
242 mutational variance of P3.p in *Caenorhabditis* arises from a broad mutational target, where
243 alterations in various biochemical processes affect its cell fate (51). We previously showed that
244 *Oscheius tipulae* displays an overall conserved pattern of Wnt expression (52). Modulation of *egl-
245 20/Wnt* dose influences the division frequency of P4.p but not that of P3.p (Fig. 4D). These results
246 suggest that the evolution of developmental bias corresponds to the evolution of developmental
247 sensitivity of the cells, which itself derives from quantitative variation in the Wnt gradient or
248 response.

249 This mode of phenotypic evolution, mediated by the evolution of nonlinear dose-response
250 curves, resembles the theoretical framework of Wright (53) and Rendel (54), describing how gene
251 activity links to phenotypes through, in this case, a single threshold. Fates of cells close to the
252 threshold exhibit high sensitivity to mutational inputs and evolve fast. Conversely, cells at the low
253 and high end of the curve exhibit low variance in response to a given perturbation, that is,
254 canalisation (55) or robustness (56). Accordingly, we propose a mechanistic model of phenotypic
255 evolution in which the developmental sensitivity of precursor cells channels phenotypic variation
256 during the evolutionary process. We expect that the distinct positions of the cells in the Wnt dose-
257 responses curve of *Caenorhabditis* and *Oscheius* are due to quantitative changes in the Wnt
258 signalling pathway or in factors acting in parallel to it. In this scenario, reduced Wnt signalling in
259 *Oscheius* vulva precursor cells would position P4.p at the fading end of the activity gradient,
260 making it the most developmentally sensitive cell, while P3.p would remain outside the sensitive
261 zone, rendering it insensitive to incoming variation. Importantly, in this model, there is no
262 discontinuity between the situations in *Caenorhabditis* and *Oscheius*, and *C. briggsae* is located
263 at a crucial transition point where P4.p starts to be affected (Fig. 3). The strength of this model is
264 that a single dose-response curve can account for evolution across different scales, within each
265 clade as well as between the two clades.

266 Altogether, our work shows how the evolution of the variational properties of cell fate
267 decisions explains the evolution of evolutionary rates for different traits and can be explained by
268 a simple shift of cells in a dose-response curve of a signalling pathway. Developmental bias

269 therefore emerges as a central mechanism connecting micro- and macro-evolutionary processes.

270 These results emphasise the importance of revealing the evolving structure of the genotype-

271 phenotype map and of probing dose-response curves in the underlying developmental system to

272 understand phenotypic variability and its evolution across evolutionary time scales.
273 **Materials and Methods**

274 **Overall design of the study**

275 We phenotyped vulval precursor cell (VPC) fates (Fig. S1) at three evolutionary scales in the same 276 environment: among species (V<sub>R</sub>, rate among species genetic variance), within species (V<sub>G</sub>, 277 standing genetic variance), and in *de novo* Random Mutant Lines (RML) together with their 278 respective ancestral control lines (V<sub>M</sub>, mutational variance) (Fig. S2). The V<sub>G</sub> and V<sub>M</sub> datasets were 279 acquired in parallel in a block design.

280 To determine the rate of phenotypic evolution in two clades, we first determined vulval 281 fates in available species of the *Caenorhabditis* and *Oscheius* genera (53 *Caenorhabditis* species 282 and all available 24 *Oscheius* species), for which we built a phylogeny from genomic and 283 transcriptomic data. For many *Caenorhabditis* and *Oscheius* species, only a single wild isolate is 284 available and we used the reference strain of each species.

285 To estimate the amount of heritable phenotypic variation within species (standing genetic 286 variance, V<sub>G</sub>, expressed in square units of phenotype), we assayed wild isolates of four species. 287 These four species were chosen as representative of the two genera because of their androdioecious 288 reproductive system (hermaphrodites and males) and large available collections of wild isolates. 289 We chose 50 isogenic wild isolates of *C. elegans*, *C. briggsae* and *O. tipulae* based on their broad 290 geographic distributions and genetic distances (57-59); and all 16 available *O. onirici* wild isolates. 291 The isogenic wild isolates were previously established by self-fertilisation of single individuals 292 isolated from the wild.

293 To estimate the amount of heritable phenotypic variation resulting from random *de novo* 294 mutation (mutational variance, V<sub>M</sub>), we first generated random mutation lines (see below and Fig. 295 S2A) in two wild isolate background of each of the four species. We assayed them in parallel with 296 the wild isolates to assess the mutationally accessible phenotypic spectrum for the six cells P3.p 297 to P8.p, in these two clades (4, 5).

298 The data were analyzed using a quantitative genetic framework accounting for non-genetic 299 factors (see below). We then used a similar quantitative genetic framework to estimate the 300 evolutionary rate among species within each genus. In this way, we account for environmental 301 variance to estimate a "rate among species genetic variance" V<sub>R</sub> (29, 35).

302 For readers unfamiliar with the quantitative genetic framework, consider these variances 303 as measures of phenotypic variation for each trait (the fate of each cell) at three evolutionary 304 timescales: in response to random mutation (M), and in intraspecific (G) and interspecific (R) 305 variation. The variance in fates of the six cells are directly comparable on the same scale.

306 Taken together, we present here findings based on the phenotypic analysis of over 80,000 307 individual animals, grown under identical laboratory conditions and examined using Nomarski 308 microscopy. The raw data for each individual is shown in Datasets S3 and S10. We represent 309 results on the scale of phenotypes in Fig. 2 and of estimated variances in Fig. 3. Supplementary 310 figures display various analyses of the data.

311 **Nematode strains and culture**

312 All nematode species and strains were maintained in standard laboratory conditions at 20°C on 313 NGM (nematode growth medium) plates with *Escherichia coli* OP50 as a food source (60). A 314 complete list of strains used in this study is provided in Dataset S1.

315 **Phylogeny**
To generate a phylogeny for Rhabditina (clade V) (Fig. S3), we used BUSCO 5.2.2(61) to analyze the genomes of 23 nematode clade V species and a clade IV outgroup (first sheet of Dataset S2). We used busco2 fasta (available at https://github.com/lstevens17/busco2 fasta) to identify and extract the protein sequences of 464 BUSCO genes (337,234 sites) that were single-copy in at least 20 (85%) of the 23 genomes. We aligned the protein sequences using MAFFT v7.49 (62) and used IQ-TREE 2.2.0.3 (63) to infer a gene tree for each BUSCO gene, allowing the best-fitting substitution model to be automatically selected (64). We provided the resulting gene trees to ASTRAL-III 5.7.8 (65) to infer the species tree. We concatenated the alignments of each protein sequence into a supermatrix using cat fasta2 phyml v1.1.0 (available at https://github.com/nylander/cat fasta2 phyml) and estimated branch lengths in substitutions per site for the ASTRAL-III species tree using IQ-TREE under the LG substitution model (66) with gamma-distributed rate variation among sites (Dataset S4).

To infer a phylogeny for *Caenorhabditis* and *Oscheius* (Fig. 1A, Fig. S4), we ran BUSCO 5.2.2 (61) on the genomes and transcriptomes of 56 *Caenorhabditis* species and 25 *Oscheius* species (second sheet in Dataset S2). We used busco2 fasta to identify and extract the protein sequences of 560 BUSCO genes (477,102 sites) that were single copy in at least 74 (91%) of the 81 genomes. We inferred a species tree as previously described: we aligned the protein sequences using MAFFT, inferred gene trees using IQ-TREE, species tree using ASTRAL-III, and branch lengths in substitutions per site using catfast2 phyml and IQ-TREE, as previously described (Dataset S5).

This considerably extends the *Caenorhabditis* phylogeny (67, 68) and of the *Oscheius* genus (24).

The set of *Oscheius* species appears somewhat less diverse at this genomic level than the set of *Caenorhabditis* species, albeit not greatly (Fig. S4).

## Dating and chronogram

We time calibrated the phylogenies in order to estimate phenotypic evolutionary rate ($V_R$, see below). We calibrated the chronogram of *Caenorhabditis* and *Oscheius* species (Fig. 1C, Figs S5, S6) by constraining the tree with three *Caenorhabditis* divergence times previously estimated (69) with the following corrections: a 14-day generation time and a corrected mutation rate ($\mu$) of 2.19x10⁻⁹ per site, per generation. This mutation rate was obtained by the weighted average of four recent studies that estimated *C. elegans* (N2) mutation rate using spontaneous mutation accumulation lines (MAL): 1.84 x10⁻⁹ with 10 MALs (70), 2.33x10⁻⁹ with 32 MALs (71), 2.7x10⁻⁹ with 10 MALs (72), and 1.33x10⁻⁹ with 7 MALs (73). These corrections resulted in the following estimated times of divergence: *C. briggsae* and *C. sinica* between 16.97 to 19.15 million years ago (MYA); *C. tropicalis* and *C. brenneri* between 20.22 to 22.36 MYA; and *C. elegans* and *C. inopinata* between 20.14 and 22.48 MYA. We estimated an ultrametric chronogram using the semiparametric method penalized likelihood (PL) with the Truncated-Newton (TN) algorithm for optimising the likelihood of branch lengths and an additive penalty implemented in the program r8s v.1.7(74) (method=pl; algorithm=tn; penalty=add). The PL method models different substitution rates among branches with a penalty for rates changing abruptly from nearby branches that is determined by a 'smoothing parameter'. We estimated the best smoothing parameter as 1 using cross-validation of a wide range of smooth parameters from 0.01 to 10,000 (crossv = yes; cvstart = -2; cvinc = 0.5; cvnum = 15). This low smoothing value indicates a strong deviation from a strict molecular clock. Moreover, we used the gradient function to check for multiple solutions using different initial starting points (checkGradient=yes, num_time_guessed=3). We obtained an estimated divergence time of 160.53 MYA between the *Caenorhabditis* and *Oscheius*
genera, 81.58 MYA for the *Oscheius* genus (*O. tipulae* and *O. myriophilus*), and 106.30 MYA for the *Caenorhabditis* genus (*C*. sp. 41 and *C. krikudae*) (see Dataset S6, Fig. 1 and Fig. S5).

We estimated the Rhabditina chronogram (Dataset S7; see Fig. S3B) in a similar manner with a best smoothing parameter of 0.1, using the following fixed divergences times from the above analysis:160.53 MYA between the *Caenorhabditis* and *Oscheius* genera; 81.58 MYA for the *O. tipulae* and *O. myriophilus*; 56.47 MYA for the *O. tipulae* and *O. dolichura*; 38.67 MYA for *C. elegans* and *C. japonica*; and 69.48 MYA for *C. japonica* and *C. angaria*.

The set of *Oscheius* species roots somewhat later than the set of *Caenorhabditis* species, between the branching dates of *C. parvicauda* and *C. plicata* (Fig. S5).

## Generation of mutation lines

We generated eight panels of random mutation lines (RML) from hermaphroditic isogenic wild isolates of two *Caenorhabditis* and two *Oscheius* species: JU1200 and PB306 for *C. elegans*, AF16 and PB800 for *C. briggsae*, CEW1 and JU178 for *O. tipulae*, and PS2068 and JU77 for *O. onirici*. The wild isolates AF16, CEW1 and PS2068 were chosen for being reference strains of the respective species (75-77). JU1200 was chosen for being the wild isolate closest genetically to the *C. elegans* reference strain N2 (78). We did not use the N2 strain because it harbours many known laboratory-adapted polymorphisms (79) that could have an impact on our study. The wild isolates PB800, JU178 and JU77 were selected for their genetic and/or phenotypic distance to the respective reference strains (23, 25, 58, 59); and PB306 because it was used in a previous study on vulval mutational variance with spontaneous mutation accumulation (MA) lines. The RMLs were generated by chemical mutagenesis (ENU) following the approach of (80) (see Fig. S2) to ensure genome-wide mutagenised lines. We used N-ethyl-N-nitrosourea (ENU) as mutagen to obtain a broader distribution of molecular lesions compared to ethylmethanesulfonate (EMS) (80), and avoid the laborious construction of spontaneous MA lines. Importantly, the molecular mutational distribution is similar between ENU and MA lines (80) and our results compare well with those in (19) for the effect of spontaneous mutation on VPC fates in *C. elegans* and *C. briggsae* mutation accumulation lines.

Briefly, P0 larvae undergoing germline proliferation and differentiation at the 3rd/4th instar were incubated with 5 mM ENU for four hours at room temperature in a rotary shaker. For the *Oscheius* strains before adding the mutagen, nematodes were incubated in 0.25% SDS for 2 minutes to increase the tissue penetrance(81) to the mutagen. Several mutagenised P0s were placed on each plate and F1 animals were screened and selected for the dominant twitching phenotype (*unc-22[-/+]) in 1% nicotine (82). This F1 selection was conducted to ensure that gametes had been effectively mutagenised since *unc-22* is a very long protein-coding gene in *C. elegans* as well as the other species and its mutant phenotype is highly specific. Next, F2 individuals were counter-selected for the twitching phenotype to generate *unc-22[+],* isogenised by 10 generations of selfing hermaphrodites and cryopreserved using standard methods (60). At each generation, a randomly early larva was picked to establish the next generation, thus avoiding selection biases. Importantly, for a given ancestral genotype, only RMLs derived from different P0 plates were phenotyped to guarantee their genetic lineage independence.

Before the mutagenesis, we generated two Ancestral Pseudo-Lines (APL) for each of the eight ancestral wild isolates (e.g. JU1200-A and JU1200-B). This was accomplished by establishing derived lines of the wild isolates from a single individual, allowed the population to expand and cryopreserved several copies of an APL as before (see below). All RMLs and APLs were generated in the laboratory in Paris.
412 Block structure and vulval phenotyping

413 The nematode phenotyping was conducted by four experimental observers: J.P.-O., N.F. and M.-
414 A.F. in Paris, and C. Bouleau in Nice.

415 The phenotyping assays to estimate the within species standing genetic variance (V<sub>G</sub>) and
416 mutational variance (V<sub>M</sub>) were performed together for V<sub>G</sub> and V<sub>M</sub> in a blind replicated block design
417 by the four experimental observers. The block design accounts for environmental effects within-
418 and among-lines and across experimental phenotyping blocks, and observer bias. All the blocks
419 had control ancestor APLs to account for environmental and non-mutational effects (4, 83) (p.
420 332) (18). In the beginning of each block, strains were given a random number for the observers
421 to phenotype strains blindly to the genotype, i.e. wild isolates, RMLs and control APLs. In total,
422 we performed one common pilot block in the two laboratories (P_PILOT and N_Pilot blocks), 10
423 blocks in Paris (P1 to P10) and 12 blocks in Nice (N1 to N12). Prior to starting the experimental
424 blocks, an intense training session on vulval phenotyping was conducted in Paris to establish
425 uniformity in ascribing vulval cell fates. Next, the pilot block was conducted to ensure phenotyping
426 consistency among observers, where the same RMLs and APLs from JU1200 (C. elegans) and
427 CEW1 (O. tipulae) were phenotyped by the four observers. In the next blocks, RMLs and wild
428 isolates were randomly allocated to blocks in the two laboratories. In the P blocks, the three
429 observers randomly phenotyped strains/replicates; while in the N blocks, all strains and replicates
430 were phenotyped by C. Bouleau. A block consisted of strains that were thawed together from APLs
431 and RMLs of a given ancestor/species, and wild isolates of the same species (from the 3rd block
432 onwards), cultured in the same batches of plates and phenotyped in sequential days. The mode
433 value of APLs, RMLs and Wild Isolates per block for a given species was two, four and three,
434 respectively. Most blocks (P1-P8, N1-6/11) had the four species; and only two blocks had species
435 from one genus. Detailed information on the block structure is in Dataset S9.
436

436 Nematode culturing and phenotyping of each block was conducted as in (19) but with two
437 replicates (see Fig. S2B). Briefly, a single individual from each thawed line (G0) was arbitrarily
438 selected to initiate the experimental population in order to eliminate potential genetic variation.
439 After two generations of population growth (G1-2), 30-50 gravid adult hermaphrodites per line
440 were cleaned of microbial contaminants by hypochlorite treatment (60). Following two discrete
441 generations of culture expansion, 15-30 adult hermaphrodites (G4) were transferred to a new plate
442 to lay the offspring to be phenotyped for replicate A. When the majority of their progeny had
443 reached the L4 stage (G5; after 2-6 days depending on the line, species or genus), vulval cell fates
444 were scored on early to mid L4 larval stages (see below). In parallel, nematode cultures were
445 maintained in standard growing conditions for 3 discrete generations (G4-6') by agar chunking,
446 then 15-30 adult hermaphrodites at G7' were transferred to a new plate to lay the offspring to be
447 phenotyped for replicate B (G8) as before. In sum, the control APLs, RMLs and the wild isolates
448 were treated in the same way and blindly, underwent approximately 5 and 8 generations for the
449 two replicates, respectively, between thawing and phenotyping on NGM plates at 20°C and low
450 density. To avoid morphological selection bias for the phenotyped specimens, nematodes were
451 washed from the plate with M9 buffer, spun-down, mounted onto an agar pad, and the first 50
452 early to mid L4 larvae were phenotyped with Nomarski microscopy with a 100x objective. The
453 pads were made of 6% noble agar and 3 mM of sodium azide (anaesthetic) (84) with a Vinyl
454 Record (85) to keep the nematodes straight to facilitate the low-throughput phenotyping. Standard
455 criteria were used to infer vulval cell fates based on the arrangement, topology, size and number
456 of cells (86-88). P(5-7).p adopt vulval fates and divide to form the adult vulva. The three other
457 cells are competent to form the vulva but normally fuse with the epidermal syncytium hyp7 either
458 before division (L2 stage) or after division (L3 stage). In *Oscheius tipulae*, the most common cell
lineage is two rounds of divisions (four cells, "SSSS") but variants with 3, 2 or 1 cell are also seen.
Here we simplify the dataset for *Oscheius* and code the phenotype as binary: two divisions versus
other fates. For further information on the rhabditid family and outgroup, see (24, 43).

The culturing and phenotyping of *Caenorhabditis* and *Oscheius* species were performed independently by two observers. In general, each species was phenotyped in two consecutive replicates of 50 specimens each. Detailed information on the number of strains and specimens phenotyped are in Tables S1, S2 and S3 below. Detailed phenotyping of all individuals is available in Datasets S3 and S10.

Quantitative genetic analysis

The analysis of the phenotyping data was performed using R (v. 4.2.3) (89) in RStudio (v. 1.4.1106) (90). To estimate quantitative genetic (QG) parameters of the Pn.p vulva precursor cells (VPC) at the three evolutionary scales ($V_R, V_G$ and $V_M$), we used the phenotypic data among lines to run Bayesian generalised linear mixed models, using the R package MCMCglmm (v. 2.36) (91).
We focus here on univariate models for each cell. The univariate models were constructed following the course notes of MCMCglmm and the tutorials of the R package QGglmm (92). We considered the six vulval traits (P3.p to P8.p) to follow a binomial distribution, where the most common fate for a given VPC was assigned a value of 1 and all others a value of 0. That is, in *Caenorhabditis* species the most common division fate for P3.p, P4.p and P8.p is one division (fate 'SS'), while in *Oscheius* species it is two divisions (fate 'SSSS'). For the central VPCs P5.p to P7.p, we considered their canonical 2º, 1º, 2º fates, respectively, as 'wt' with a value of 1. Given the binomial distribution of all traits, we used the threshold model (53) with a probit link-function that has a better mixing of the Markov Chain Monte Carlo (MCMC) chains (93). To ensure large MCMC effective sample sizes of the posterior distribution, we chose 5,010,000 and 1,260,000 iterations, and 2,000 and 500 thinning intervals for the analyses of $V_R$ and $V_G \& M$, respectively.
This resulted in minimum samples sizes of genetic variance in the liability scale of 2,302 for $V_R$, 2,244 for $V_G$, and 2,194 for $V_M$. In all models, burn-in was set to 10,000 iterations and convergence was checked by visual inspection of the posterior distributions and by the Heidelberg and Welch (94) convergence test implemented in the coda R package (v. 0.19-4.1) (95). Furthermore, to contribute for model convergence, latent variables were truncated (trunc=TRUE) to prevent under/overflow (-/+7) of the models. The posterior distribution of random effects was saved (pr=TRUE) to obtain the best linear unbiased predictions (BLUP) of the nematode strains for each trait, using the R package postMCMCglmm (v. 0.1-2) (96). We used non-informative priors: a wide normal distribution on all fixed effects, and a chi-square distribution with one degree of freedom (V = 1, nu = 1000, alpha.mu = 0, alpha.V = 1) on all random effects due to its better inference for binary traits (97). The residual variance had to be fixed to 1 because it is not identifiable in binary responses (see MCMCglmm course notes of and QGglmm tutorial). From these models, we extracted the following QG parameters in the liability scale using their posterior distributions: trait mean (z), genetic variance ($V_g$), phenotypic variance ($V_P$) and broad-sense heritability ($H^2 = V_g / V_P$). The trait mean is the overall mean of the fixed effects of the model. The $V_g$ is half the among-line variance because it is obtained from homozygous isogenic lines that include non-additive effects such as dominance (98). The phenotypic variance is the sum of the variances from all random effects: $V_G$, block and replicate effects ($V_{Brep}$), residual variance of 1, and variance explained by fixed effects (99). These QG parameters were then computed to the observed data scale with the QGglmm package (v. 0.7.34) (92), using the binomial error distribution and integrating over the posterior distribution of each parameter. We extracted for all QG parameters the mean, median, posterior mode, and the 95% and 83% credible intervals from
their posterior distributions. Significant differences of posterior distributions can be inferred when
their 83% credible intervals do not overlap (100). Furthermore, to measure the fold-change
differences in genetic variance between P3.p and P4.p, we divided their posterior distributions and
the resulting ratio displayed in a log10 scale. We did not mean-scale the variance because binomial
traits are already in the same scale and there is an intrinsic relationship between mean and variance.

To estimate the evolutionary rate among species (V<sub>R</sub>), we additionally used the
chronogram of each nematode genus separately to run generalised linear phylogenetic mixed
models (30) following the guidelines of (101). Although this estimate only measures a certain
pattern in interspecific variation (102), the goal here is to remove the environmental variance and
the influence of possible intraspecific polymorphisms by focusing on the part of the variation that
is explained by the phylogeny. Briefly, the same logic of the QG animal models was used, but
instead of a pedigree relatedness covariance matrix, we used a phylogenetic covariance matrix.
This matrix assumes a Brownian motion model of trait evolution where there is phylogenetic
inertia and a constant increase in trait variance through evolution. In this way, phenotypic
evolutionary rates are estimated as increases of genetic variance per million years (V<sub>R</sub>) (103).
Additionally, phylogenetic heritability (or phylogenetic signal) is the proportion of phenotypic
variance explained by the phylogenetic relationships (H<sup>2</sup><sub>R</sub> = V<sub>R</sub> / V<sub>P</sub>) (29, 102) instead of pedigree-
based heritability. We first generated the inverse phylogenetic covariance matrix using all the
nodes and scaled the total branch length, from root to tips, to one (nodes="ALL", scale=TRUE
with the inverseA::MCMCglmm function). In these models we did not have fixed effects, and had
the following random effect: variance among-species dependent on the phylogeny (V<sub>R1</sub>), variance
among-species not dependent on the phylogeny (V<sub>NotPhylo</sub>), variance among species replicates
(V<sub>replicate</sub>). In this case, V<sub>P</sub> = V<sub>R</sub> + V<sub>NotPhylo</sub> + V<sub>replicate</sub> + residual variance of 1. Since the trees were
scaled in the analyses, all the random factors are in a standardise scale. To obtain the V<sub>R</sub> per
million years (MY), we divided the V<sub>R1</sub> by the total length of the chronogram: 106.30 MY for
Caenorhabditis and 81.58 MY for Oscheius.

To map the trait evolution along each phylogeny (Fig. S6), we obtained the mean trait
frequency per species and estimated ancestral states using a Brownian motion model, using the
contMap function in the R package Phytools (v. 2.3-0) (104).

The standing genetic variance (V<sub>G</sub>) is the genetic variance among wild isolates. The
mutational variance (V<sub>M</sub>) is the genetic variance among random mutant lines (RML). When
estimating these variances, we included the variance among control ancestral pseudo-lines (APL)
of the respective species to account for uncontrolled non-genetic variation (4, 83) (p. 332). For
this, we included a fixed effect named *Treatment* for control vs. mutation/wild. We ran models
only with the control APLs, the wild isolates and the RMLs, and compared them to the full V<sub>G</sub> and
V<sub>M</sub> models. We observed significant V<sub>G</sub> and V<sub>M</sub> for P3.p in all *Caenorhabditis* species and isolates,
except in the JU1200 V<sub>M</sub>; and for P4.p in all *Oscheius* species and isolates (see Figs S10, S15). In
all V<sub>M</sub> models, we ran models with isolates independently, and with the two isolates together for
each species. This species V<sub>M</sub> was estimated by including an *Ancestral* (e.g. PB306 vs. JU1200)
fixed effect that was nested within *Treatment*. In all aforementioned V<sub>G</sub> and V<sub>M</sub> models, we
included the four experimental observers (J.P.-O., N.F. and M.-A.F., and C.B.-M) as fixed effect
to take into account observational phenotyping biases. Also, in all these models we had one random
effect of block and replicate (V<sub>Brep</sub>) that depicted the experimental structure of the phenotyping.

To estimate mutational effects on the directional changes in the trait means of RMLs in
respect to their ancestors, we estimated marginal means of the fixed effect *Treatment* (control Vs
mutation) using the R package emmeans (v. 1.10.3) (105).

All plots were made using the R package ggplot2 (v. 3.5.1) (106).
**Oscheius tipulae egl-20 experiment**

The *Oti-egl-20* mutant (52) was scored in the same experimental conditions as above, in parallel
with the control reference wild-type *O. tipulae* background CEW1.

**Data availability:** All raw phenotyping data are available in the Supplementary Information.
Links to genome assemblies used in the phylogenetic analyses are provided in Dataset S2.

**Code availability:** The data analysis code and plots of the posterior distributions of the models
are available at https://github.com/Picao-Osorio/Evolution-Developmental-Bias/tree/main.

**Material availability:** Strains are available upon inquiry to MAF: felix@bio.ens.psl.eu.

**Competing interests** The authors declare no competing interest.

**Acknowledgements and funding sources**

We thank John Wang, Takao Inoue, Erik Andersen, Yen-Ping Hsueh, Walter Sudhaus, David Fitch
for strains, Karin Kiontke for help in *Oscheius* species collection and delineation, Hervé Gendrot
for plate pouring. We are grateful to Pierre de Villemereuil, Katie Pelletier, Marie Marcaillou,
Henrique Teotónio, François Mallard, Luke Noble, Archishman Raju, BingKan Xue and Stanislas
Leibler for discussions and/or providing useful comments on the manuscript. We thank the two
anonymous reviewers whose constructive and positive criticisms improved the quality and scope
of this article. Some strains were provided by the CGC, which is funded by NIH Office of Research
Infrastructure Programs (P40 OD010440). This research was funded in part, by grants from
Agence Nationale de la Recherche (France) 18-CE13-0006-01 and 22-CE13-0005-01 (MAF, CB),
Marie Skłodowska Curie Training Grants (European Union) 751530-EvoCellFate and 101110355-
EvoBias (JPO), Wellcome Trust Grant 220540/Z/20/A (MB, LS, PGR).
584 1. A. Holstad *et al.*, Evolvability predicts macroevolution under fluctuating selection. *Science* **384**, 688-693 (2024).

585 2. M. Tsuboi *et al.*, The paradox of predictability provides a bridge between micro- and macroevolution. *J Evol Biol* **37**, 1413-1432 (2024).

586 3. G. G. Simpson, *Tempo and Mode in Evolution* (Columbia University Press, New York, USA, 1944).

587 4. M. Lynch, B. Walsh, *Genetics and analysis of quantitative traits* (Sinauer Associates, Inc., Sunderland, MA, 1998), pp. 803.

588 5. J. Maynard Smith *et al.*, Developmental constraints and evolution. *Q Rev Biol* **60**, 265-287 (1985).

589 6. M. Lynch, The frailty of adaptive hypotheses for the origins of organismal complexity. *Proc Natl Acad Sci USA* **104** **Suppl 1**, 8597-8604 (2007).

590 7. A. Stoltzfus, L. Y. Yampolsky, Climbing mount probable: mutation as a cause of nonrandomness in evolution. *J Hered* **100**, 637-647 (2009).

591 8. W. Arthur, *Biased embryos and evolution* (Cambridge University Press, Cambridge, UK, 2004), pp. 248.

592 9. T. Uller, A. P. Moczek, R. A. Watson, P. M. Brakefield, K. N. Laland, Developmental bias and evolution: a regulatory network perspective. *Genetics* **209**, 949-966 (2018).

593 10. H. Cai, D. Melo, D. L. Des Marais, Disentangling variational bias: the roles of development, mutation, and selection. *Trends Genet* **41**, 23-32 (2025).

594 11. J. M. Cheverud, Quantitative genetics and developmental constraints on evolution by selection. *J Theor Biol* **110**, 155-171 (1984).

595 12. S. J. Arnold, Constraints on phenotypic evolution. *Am Nat* **140** **Suppl 1**, S85-107 (1992).

596 13. D. Schluter, Adaptive radiation along genetic lines of least resistance. *Evolution* **50**, 1766-1774 (1996).

597 14. G. P. Wagner, C.-h. Chiu, M. Laublicher, Developmental evolution as a mechanistic science: The inference from developmental mechanisms to evolutionary processes. *Amer Zool* **40**, 819-831 (2000).

598 15. T. F. Hansen, D. Houle, M. Pavlicev, C. Pélabon, *Evolvability - A unifying concept in evolutionary biology?*, Vienna Series in Theoretical Biology (MIT Press, 2023).

599 16. D. Houle, How should we explain variation in the genetic variance of traits? *Genetica* **102-103**, 241-253 (1998).

600 17. D. Houle, J. Fierst, Properties of spontaneous mutational variance and covariance for wing size and shape in *Drosophila melanogaster*. *Evolution* **67**, 1116-1130 (2013).

601 18. R. Farhadifar, J. M. Ponciano, E. C. Andersen, D. J. Needleman, C. F. Baer, Mutation is a sufficient and robust predictor of genetic variation for mitotic spindle traits in *Caenorhabditis elegans*. *Genetics* **203**, 1859-1870 (2016).

602 19. C. Braendle, C. F. Baer, M. A. Félix, Bias and evolution of the mutationally accessible phenotypic space in a developmental system. *PLoS Genet* **6**, e1000877 (2010).

603 20. P. W. Sternberg (2005) Vulval development. in *Wormbook*, ed The *C. elegans* Research Community.

604 21. J. Sulston, H. R. Horvitz, Postembryonic cell lineages of the nematode *Caenorhabditis elegans*. *Dev. Biol.* **56**, 110-156 (1977).
627 22. C. Braendle, M.-A. Félix, Plasticity and errors of a robust developmental system in different environments. *Dev Cell* **15**, 714-724 (2008).

629 23. M. Delattre, M.-A. Félix, Polymorphism and evolution of vulval precursor cell lineages within two nematode genera, *Caenorhabditis* and *Oscheius*. *Curr Biol* **11**, 631-643 (2001).

632 24. K. Kiontke *et al.*, Trends, stasis and drift in the evolution of nematode vulva development. *Curr. Biol.* **17**, 1925-1937 (2007).

634 25. J.-B. Pénigault, M.-A. Félix, Evolution of a system sensitive to stochastic noise: P3.p cell fate in *Caenorhabditis*. *Dev Biol* **357**, 419-427 (2011).

636 26. D. M. Eisenmann, J. N. Maloof, J. S. Simske, C. Kenyon, S. K. Kim, The beta-catenin homolog BAR-1 and LET-60 Ras coordinately regulate the Hox gene *lin-39* during *Caenorhabditis elegans* vulval development. *Development* **125**, 3667-3680 (1998).

639 27. S. Louvet-Vallée, I. Kolotuev, B. Podbilewicz, M.-A. Félix, Control of vulval competence and centering in the nematode *Oscheius* sp. 1 CEW1. *Genetics* **163**, 133-146 (2003).

642 28. R. J. Sommer, P. W. Sternberg, Evolution of cell lineage and pattern formation in the vulval equivalence group of rhabditid nematodes. *Dev. Biol.* **167**, 61-74 (1995).

644 29. E. A. Housworth, E. P. Martins, M. Lynch, The phylogenetic mixed model. *Am Nat* **163**, 84-96 (2004).

646 30. M. Lynch, Methods for the analysis of comparative data in evolutionary biology. *Evolution* **45**, 1065-1080 (1991).

648 31. P. D. Gingerich, *Rates of Evolution: A Quantitative Synthesis* (Cambridge University Press, Cambridge, UK, 2019).

650 32. J. C. Uyeda, T. F. Hansen, S. J. Arnold, J. Pienaar, The million-year wait for macroevolutionary bursts. *Proc Natl Acad Sci USA* **108**, 15908-15913 (2011).

652 33. J. Rolland *et al.*, Conceptual and empirical bridges between micro- and macroevolution. *Nat Ecol Evol* **7**, 1181-1193 (2023).

654 34. T. F. Hansen, Three modes of evolution? Remarks on rates of evolution and time scaling. *J. Evol. Biol.* **37**, 1523-1537 (2024).

656 35. D. Houle, G. H. Bolstad, K. van der Linde, T. F. Hansen, Mutation predicts 40 million years of fly wing evolution. *Nature* **548**, 447-450 (2017).

658 36. P. T. Rohner, D. Berger, Developmental bias predicts 60 million years of wing shape evolution. *Proc Natl Acad Sci USA* **120**, e22112-120120 (2023).

660 37. R. Farhadifar *et al.*, Scaling, selection, and evolutionary dynamics of the mitotic spindle. *Curr Biol* **25**, 732-740 (2015).

662 38. W. G. Hill, Rates of change in quantitative traits from fixation of new mutations. *Proc Natl Acad Sci USA* **79**, 142-145 (1982).

664 39. A. G. Jones, S. J. Arnold, R. Bürger, The mutation matrix and the evolution of evolvability. *Evolution* **61**, 727-745 (2007).

666 40. A. G. Jones, R. Burger, S. J. Arnold, Epistasis and natural selection shape the mutational architecture of complex traits. *Nat Commun* **5**, 3709 (2014).

668 41. A. G. Jones, S. J. Arnold, B. R., Stability of the G-matrix in a population experiencing pleiotropic mutation. *Evolution* **57**, 1747-1760 (2003).

670 42. S. J. Arnold, R. Burger, P. A. Hohenlohe, B. C. Ajie, A. G. Jones, Understanding the evolution and stability of the G-matrix. *Evolution* **62**, 2451-2461 (2008).

672 43. M.-A. Félix *et al.*, Evolution of vulva development in the Cephalobina (Nematoda). *Dev. Biol.* **221**, 68-86 (2000).
674 44. H. Zauner, R. J. Sommer, Evolution of robustness in the signaling network of *Pristionchus* vulva development. *Proc Natl Acad Sci U S A* **104**, 10086-10091 (2007).
675
676 45. H. Sawa, H. C. Korswagen (2013) Wnt signaling in *C. elegans*. in *WormBook*, ed The *C. elegans* Research Community.
677
678 46. T. R. Myers, I. Greenwald, Wnt signal from multiple tissues and lin-3/EGF signal from the gonad maintain vulval precursor cell competence in *Caenorhabditis elegans*. *Proc. Natl. Acad. Sci. USA* **104**, 20368-20373 (2007).
679
680
681 47. J.-B. Pénigault, M.-A. Félix, High sensitivity of *C. elegans* vulval precursor cells to the dose of posterior Wnts. *Dev Biol* **357**, 428-438 (2011).
682
683 48. M. Harterink *et al.*, Neuroblast migration along the anteroposterior axis of *C. elegans* is controlled by opposing gradients of Wnts and a secreted Frizzled-related protein. *Development* **138**, 2915-2924 (2011).
684
685
686 49. A. M. Pani, B. Goldstein, Direct visualization of a native Wnt in vivo reveals that a long-range Wnt gradient forms by extracellular dispersal. *eLife* **7**, e38325 (2018).
687
688 50. J. R. Kroll, J. Tsiaxiras, J. S. van Zon, Variability in beta-catenin pulse dynamics in a stochastic cell fate decision in *C. elegans*. *Dev Biol* **461**, 110-123 (2020).
689
690 51. F. Besnard, J. Picao Osorio, C. Dubois, M.-A. Félix, A broad mutational target explains a fast rate of phenotypic evolution. *eLife* **9**, e54928 (2020).
691
692 52. A. M. Vargas-Velazquez, F. Besnard, M. A. Felix, Necessity and contingency in developmental genetic screens: EGF, Wnt, and semaphorin pathways in vulval induction of the nematode *Oscheius tipulae*. *Genetics* **211**, 1315-1330 (2019).
693
694
695 53. S. Wright, Molecular and evolutionary theories of dominance. *Am. Naturalist* **63**, 24-53 (1934).
696
697 54. J. M. Rendel, The relationship between gene and phenotype. *J. Theor. Biol.* **2**, 296-308 (1962).
698
699 55. J. M. Rendel, *Canalisation and gene control* (Logos Press, 1967), pp. 166.
700 56. M. A. Félix, M. Barkoulas, Pervasive robustness in biological systems. *Nature reviews. Genetics* **16**, 483-496 (2015).
701
702 57. D. E. Cook, S. Zdraljevic, J. P. Roberts, E. C. Andersen, CeNDR, the *Caenorhabditis elegans* natural diversity resource. *Nucleic Acids Res* **45**, D650-D657 (2017).
703
704 58. A. D. Cutter, Nucleotide polymorphism and linkage disequilibrium in wild populations of the partial selfer *Caenorhabditis elegans*. *Genetics* **172**, 171-184 (2006).
705
706 59. D. Baïlle, A. Barrière, M.-A. Félix, *Oscheius tipulae*, a widespread hermaphroditic soil nematode, displays a higher genetic diversity and geographical structure than *Caenorhabditis elegans*. *Mol. Ecol.* **17**, 1523-1534 (2008).
707
708
709 60. T. Stiernagle (2006) Maintenance of *C. elegans*. in *Wormbook*, ed The *C. elegans* Research Community.
710
711 61. M. Manni, M. R. Berkeley, M. Seppey, E. M. Zdobnov, BUSCO: Assessing Genomic Data Quality and Beyond. *Curr Protoc* **1**, e323 (2021).
712
713 62. R. K. Bradley *et al.*, Fast statistical alignment. *PLoS Comput Biol* **5**, e1000392 (2009).
714 63. B. Q. Minh *et al.*, IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. *Mol Biol Evol* **37**, 1530-1534 (2020).
715
716 64. S. Kalyaanamoorthy, B. Q. Minh, T. K. F. Wong, A. von Haeseler, L. S. Jermiin, ModelFinder: fast model selection for accurate phylogenetic estimates. *Nat Methods* **14**, 587-589 (2017).
717
718
719 65. C. Zhang, M. Rabiee, E. Sayyari, S. Mirarab, ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. *BMC bioinformatics* **19**, 153 (2018).
720
721 66. S. Q. Le, N. Lartillot, O. Gascuel, Phylogenetic mixture models for proteins. *Philos Trans R Soc Lond B Biol Sci* **363**, 3965-3976 (2008).
722
723 67. L. Stevens *et al.*, Comparative genomics of ten new *Caenorhabditis* species. *Evolution Letters* **3**, 217-236 (2019).
724
725 68. S. A. Sloat *et al.*, *Caenorhabditis* nematodes colonize ephemeral resource patches in neotropical forests. *Ecology and evolution* **12**, e9124 (2022).
726
727 69. N. Kanzaki *et al.*, Biology and genome of a newly discovered sibling species of *Caenorhabditis elegans*. *Nat Commun* **9**, 3216 (2018).
728
729 70. A. Konrad, M. J. Brady, U. Bergthorsson, V. Katju, Mutational landscape of spontaneous base substitutions and small indels in experimental *Caenorhabditis elegans* populations of differing Size. *Genetics* **212**, 837-854 (2019).
730
731
732 71. A. S. Saxena, M. P. Salomon, C. Matsuba, S. D. Yeh, C. F. Baer, Evolution of the mutational process under relaxed selection in *Caenorhabditis elegans*. *Mol Biol Evol* **36**, 239-251 (2019).
733
734
735 72. D. R. Denver *et al.*, A genome-wide view of *Caenorhabditis elegans* base-substitution mutation processes. *Proc Natl Acad Sci USA* **106**, 16310-16314 (2009).
736
737 73. D. R. Denver *et al.*, Variation in base-substitution mutation in experimental and natural lineages of *Caenorhabditis* nematodes. *Genome Biol Evol* **4**, 513-522 (2012).
738
739 74. M. J. Sanderson, r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. *Bioinformatics* **19**, 301-302 (2003).
740
741 75. P. W. Sternberg *et al.*, WormBase 2024: status and transitioning to Alliance infrastructure. *Genetics* **227** (2024).
742
743 76. M.-A. Félix, A. Vierstraete, J. Vanfleteren, Three biological species related to *Rhabditis* (*Oscheius*) pseudodolichura Körner in Osche, 1952. *J. Nematol.* **33**, 104-109 (2001).
744
745 77. M.-A. Félix (2006) *Oscheius tipulae*. in *WormBook*, ed The *C. elegans* Research Community.
746
747 78. E. Andersen *et al.*, Chromosome-scale selective sweeps shape *Caenorhabditis elegans* genomic diversity. *Nature genetics* **45**, 285-290 (2012).
748
749 79. M. G. Sterken, L. B. Snoek, J. E. Kammenga, E. C. Andersen, The laboratory domestication of *Caenorhabditis elegans*. *Trends Genet* **31**, 224-231 (2015).
750
751 80. O. Thompson *et al.*, The million mutation project: a new approach to genetics in *Caenorhabditis elegans*. *Genome Res* **23**, 1749-1762 (2013).
752
753 81. S. Zhang, D. Banerjee, J. R. Kuhn, Isolation and culture of larval cells from *C. elegans*. *PLoS One* **6**, e19505 (2011).
754
755 82. D. G. Moerman, D. L. Baillie, Genetic organization in *Caenorhabditis elegans*: Fine-structure analysis of the *unc-22* gene. *Genetics* **91**, 95-104 (1979).
756
757 83. M. Lynch, Spontaneous mutations for life-history characters in an obligate parthenogen. *Evolution* **39**, 804-818 (1985).
758
759 84. W. B. Wood, *The nematode Caenorhabditis elegans* (Cold Spring Harbor Laboratory, Cold Spring Harbor, New York, 1988), pp. 667.
760
761 85. K. Rivera Gomez, M. Schvarzstein, Immobilization of nematodes for live imaging using an agarose pad produced with a Vinyl Record. *MicroPubl Biol* **2018** (2018).
762
763 86. W. S. Katz, R. J. Hill, T. R. Clandinin, P. W. Sternberg, Different levels of the *C. elegans* growth factor LIN-3 promote distinct vulval precursor fates. *Cell* **82**, 297-307 (1995).
764
765 87. M.-L. Dichtel, S. Louvet-Vallée, M. E. Viney, M.-A. Félix, P. W. Sternberg, Control of vulval cell division number in the nematode *Oscheius/Dolichorhabditis* sp. CEW1. *Genetics* **157**, 183-197 (2001).
766
767
768 88. M.-A. Félix, Cryptic quantitative evolution of the vulva intercellular signaling network in *Caenorhabditis. Curr. Biol.* **17**, 103-114 (2007).

770 89. RCoreTeam (2023) R: A language and environment for statistical computing., ed R Foundation for Statistical Computing (Vienna, Austria.).

772 90. RStudioTeam (2023) RStudio: Integrated Development for R. (RStudio, PBC Boston, MA).

774 91. J. D. Hadfield, MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R Package. *Journal of Statistical Software* **33** (2010).

776 92. P. de Villemereuil, H. Schielzeth, S. Nakagawa, M. Morrissey, General methods for evolutionary quantitative genetic inference from generalized mixed models. *Genetics* **204**, 1281-1294 (2016).

779 93. J. D. Hadfield, MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package. *Journal of Statistical Software* **33**, 1-22 (2010).

781 94. P. Heidelberger, P. D. Welch, A spectral method for confidence interval generation and run length control in simulations. *Communications of the ACM* **24**, 347-360 (1981).

782 95. M. Plummer, N. Best, K. Cowles, K. Vines, CODA: Convergence Diagnosis and Output Analysis for MCMC. *R News* **6**, 7-11 (2006).

785 96. J. F. Wiley (2013) postMCMCglmm: Average marginal predicted probabilities from Bayesian ordered probit models. ed R. p. v. 0.1-2.

787 97. P. de Villemereuil, O. Gimenez, B. Doligez, *Methods in Ecology and Evolution* **4**, 260-275 (2013).

789 98. D. S. Falconer, T. F. C. Mackay, *An introduction to quantitative genetics* (Pearson Education Ltd, Harlow, UK, 1996), vol. 4th edition, pp. 464.

791 99. P. de Villemereuil, M. B. Morrissey, S. Nakagawa, H. Schielzeth, Fixed-effect variance and the estimation of repeatabilities and heritabilities: issues and solutions. *J Evol Biol* **31**, 621-632 (2018).

794 100. P. C. Austin, J. E. Hux, A brief note on overlapping confidence intervals. *J Vasc Surg* **36**, 194-195 (2002).

796 101. P. de Villemereuil, S. Nakagawa, "General quantitative genetic methods for comparative biology" in Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology L. Z. Garamszegi, Ed. (Springer, Berlin, Heidelberg, 2014), pp. 287-303.

800 102. L. J. Revell, L. J. Harmon, D. C. Collar, Phylogenetic signal, evolutionary process, and rate. *Syst Biol* **57**, 591-601 (2008).

802 103. L. J. Revell, L. J. Harmon, *Phylogenetic comparative methods in R* (Princeton University Press, 2022).

804 104. L. Revell, phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and other things). *Peer J.* **12**, e16505 (2024).

806 105. R. Lenth *et al*. (2023) emmeans: Estimated Marginal Means, aka Least-Squares Means.

807 106. H. Wickham, *ggplot2: Elegant Graphics for Data Analysis*. (Springer-Verlag New York, 2016).

809 107. F. Besnard, G. Koutsovoulos, S. Dieudonné, M. Blaxter, M.-A. Félix, Toward universal forward genetics: Using a draft genome sequence of the nematode *Oscheius tipulae* to identify mutations affecting vulva development. *Genetics* **206**, 1747-1761 (2017).

812
813
814
815
816
817 **Author contributions** JPO, MAF, and CBr conceived the study. JPO made the mutation lines and designed the phenotyping assays. JPO, CBo, NF, and MAF phenotyped the animals. PGR, LS and MB sequenced and assembled the genome data and obtained the phylogenetic tree. JPO analyzed the data. MAF, CBr, and MB contributed funding and supervision. JPO and MAF wrote the initial draft. CBr and LS edited the manuscript.
818
819
820
821
822
823
Figure legends

<caption>**Figure 1.** Two nematode genera display divergent patterns of vulva cell fate evolution.</caption>

(A) Diagram of a L1 nematode larva with the six vulva precursor cells (VPC) highlighted in green (modified from Openworm.org). Additional Pn.p cells are shown in grey. Anterior is to the left and ventral side down.
(B) Nomarski images of a mid-L4 stage vulva displaying the VPC descendants in *C. elegans* and *O. tipulae*. In the *C. elegans* individual (top), P3.p, P4.p and P8.p all undergo a division as shown by the presence of two nuclei per mother Pn.p cell (arrowheads). In the *O. tipulae* individual (bottom), P3.p does not divide and P4.p and P8.p divide twice, producing four nuclei (arrowheads). The P3.p and P4.p lineages with no division, 2, 3 or 4 progeny cells are detailed below and color code in C. Pictures from (107). Scale bar: 10 µm.

(C) Horizontal bar plots showing division fates of the P3.p, P4.p and P5.p cells plotted along the chronogram tree of *Caenorhabditis* (red) and *Oscheius* (blue). See Fig. S5 and Dataset S3 for division patterns for the six cells. Myr: Million year (estimate). Stars mark the species used in further studies. 'S' is syncytial epidermal fate and the number of 'S' cells indicate the division pattern of the mother cell. P5.p is shown as a contrast that almost invariantly adopts a 2° fate lineage (black). As an exception, the 2° fated cells in *Oscheius* sp. L often adopt an additional division of one or more of the grand-daughters in about half of the animals (green, Dataset S3). (D) Schematic representation of the experimental design scoring vulval cell fates in two nematode genera across three evolutionary time scales: 1. among species (*V*<sub>R</sub>; top); 2. within species (*V*<sub>G</sub>; middle); and 3. upon *de novo* mutation (*V*<sub>M</sub>; bottom). Numbers of lines and individuals phenotyped in the present work are indicated.

<caption>**Figure 2.** Each nematode genus shows a specific pattern of variation for P3.p and P4.p cell fates across evolutionary scales.</caption>

(A) Division frequency of P3.p and P4.p in species of *Caenorhabditis* (salmon) and *Oscheius* (dark blue). Each dot denotes the average division frequency of a species.

(B) Division frequency of P3.p and P4.p in wild isolates of the four species as in (A). Each dot represents the Best linear unbiased prediction (BLUP) for each wild isolate.

(C) Division frequency of P3.p (x-axis) and P4.p (y-axis) upon random mutation in two isolates of *C. elegans* (red), *C. briggsae* (pink), *O. tipulae* (blue) and *O. onirici* (turquoise). Each dot represents the Best linear unbiased prediction (BLUP) in the observed data scale for each strain, corrected for block and observer structure. Darker colours are the pseudo-ancestral control lines and lighter colours the random mutant lines.

On top and on the right side of each plot are boxplots displaying the data distribution for P3.p and P4.p, respectively. Boxplots represent median values, interquartile ranges, Tukey whiskers, and outliers as individual points.

Numbers of phenotyped individuals and lines are provided on Fig. 1D and methods. As in Fig. 1C, is displayed the frequency of the most common type of cell lineage in each genus.

<caption>**Figure 3.** Each nematode genus shows consistent genetic variances for vulval traits across evolutionary scales.</caption>

(A) Rate among species genetic variance (*V*<sub>R</sub>) of the three anterior VPCs (P3.p to P5.p) estimated for the *Caenorhabditis* and *Oscheius* genera. *V*<sub>R</sub> quantifies phenotypic evolutionary rates by estimating the increase of genetic variance among species per million years assuming a Brownian motion model of trait evolution. The scale of the genetic variance (y-axis) is the same for both graphs.

(B) Standing genetic variances (*V*<sub>G</sub>) of the three anterior VPCs in *C. elegans* (red), *C. briggsae* (pink), *O. tipulae* (blue) and *O. onirici* (turquoise) estimated from the variation among wild isolates adjusted for the intrastrain variation of two APL controls on each species.

(C) Mutational variances (*V*<sub>M</sub>) of the three anterior VPCs in *C. elegans* (red), *C. briggsae* (pink), *O. tipulae* (blue) and *O. onirici* (turquoise). The mutational variance is estimated as the genetic variance among random mutation lines of two isolates corrected for their ancestral pseudo-line (APL) controls.

The scale of the genetic variance (y-axis) is customised to each species in (B) and (C).
Each estimation displays the median (black circle), 83% (coloured bar) and 95% (black bar) credible intervals of the posterior distributions of the genetic variances in the observed data scale. Each posterior distribution has a minimum effective sample size of 2000. See Fig. S17 for estimates for the six cells.

**Figure 4.** Divergence in developmental bias matches divergence in phenotypic evolutionary rates between the two nematode clades.

(A) Ten-fold change of the ratio between P3.p and P4.p genetic variances across $V_M$, $V_G$ and $V_R$ in the two genera. Bar plots display the median of the ratio of the posterior distributions of P3.p over P4.p variances, with the 83% (coloured bar) and 95% (black bar) credible intervals, plotted on a log scale.

(B) Model of the genotype-phenotype map across evolutionary time scales that schematizes the data. From bottom to top, the contrasting developmental bias in *Caenorhabditis* and *Oscheius* revealed by the orthogonal phenotypic variation upon *de novo* mutation ($V_M$) is aligned with the axes of heritable phenotypic variation within ($V_G$) and among ($V_R$) species of the respective genus. A highly schematic frequency distribution for the starting genotype (indicated in light blue and red in $V_M$ diagram) and for the population of mutated genotypes ($V_M$), intraspecific variation ($V_G$) and species ($V_R$) is indicated in the vertical axis in grey tone.

(C) Diagram of an early *C. elegans* larva with the six VPCs (green) and the posterior Wnt protein gradient (based on (49), with a simplified shape). P3.p is placed at the fading end of the gradient. This provides a simple one-dimensional underlying variable explaining the observed variational pattern for P3.p, P4.p and P5.p.

(D) Bar plots showing the dose-response effects of *egl-20/Wnt* genetic modulation on the division frequency of P3.p, P4.p and P5.p in *C. elegans* (N2 genetic background) and *O. tipulae* (CEW1 genetic background). The *C. elegans* data are replotted from (47). For simplicity, we neglected the few cases of miscentering of the 1° fate on P5.p. O.E.: overexpression.