Selection and the direction of phenotypic evolution
Abstract
Predicting adaptive phenotypic evolution depends on invariable selection gradients and on the stability of the genetic covariances between the component traits of the multivariate phenotype. We describe the evolution of six traits of locomotion behavior and body size in the nematode Caenorhabditis elegans for 50 generations of adaptation to a novel environment. We show that the direction of adaptive multivariate phenotypic evolution can be predicted from the ancestral selection differentials, particularly when the traits were measured in the new environment. Interestingly, the evolution of individual traits does not always occur in the direction of selection, nor are trait responses to selection always homogeneous among replicate populations. These observations are explained because the phenotypic dimension with most of the ancestral standing genetic variation only partially aligns with the phenotypic dimension under directional selection. These findings validate selection theory and suggest that the direction of multivariate adaptive phenotypic evolution is predictable for tens of generations.
Editor's evaluation
This is an important paper that takes advantage of a comprehensive evolutionary genetic dataset to tease apart the relationship between genetic variation, selection, and phenotypic divergence over 50 generations. The evidence supporting the conclusions is robust and aligns with a growing body of work that shows patterns of variation can predict divergence over long periods of time and also that evolution does not always occur in the direction of selection, particularly when selection is acting on genetically correlated traits. The questions addressed in this study will particularly appeal to evolutionary biologists and quantitative geneticists.
https://doi.org/10.7554/eLife.80993.sa0Introduction
Predicting adaptive phenotypic evolution is an important research goal in evolutionary biology, agronomy, or in conservation policy (Arnold, 2014; Nosil et al., 2020; Wortel et al., 2023). It is generally accepted that predicting adaptive phenotypic evolution should be done in the context of the whole organism because organisms are not mere collections of genetically or environmentally independent traits (Gould and Lewontin, 1979; Wagner, 2001). Many traits in natural populations are heritable across generations and under natural selection (Walsh and Lynch, 2018). Often, however, phenotypic evolution is not observed (Merilä et al., 2001; Pujol et al., 2018), or is of opposite direction than predicted, because of environmental or genetic correlations between the traits of interest with unmeasured traits (Morrissey et al., 2012; Hajduk et al., 2020).
Current approaches to predict phenotypic evolution during adaptation to novel environments, at least for infinite sexual populations and under infinitesimal assumptions of trait inheritance (Barton et al., 2017), rely on Simpson’s adaptive landscape idea and its formalization by R. Lande: (Simpson, 1944; Lande, 1979; Arnold et al., 2001). In Lande’s equation, the adaptive evolution of multiple traits’ means over one generation (the column vector ) is modeled as a function of the -matrix, which summarizes the genetic structure and heritable transmission of traits from parents to offspring, with the additive genetic variances of traits being represented in the diagonal entries and the additive genetic covariances between them in the off-diagonal entries. The direction and magnitude of phenotypic evolution depend on the size, due to genetic variances, and shape, due to genetic covariances, of the -matrix. In particular, trait combinations with more genetic variation (henceforth, called canonical traits), allow for faster and more directed responses to selection (Fisher, 1930; Lande, 1980; Schluter, 1996; Blows and McGuigan, 2015). In Lande’s equation, phenotypic evolution is also modeled as a function of the selection gradients on each trait (the vector ), gradients that describe the strength of directional selection on each trait in populations that can be under stabilizing/disruptive selection on each trait and under correlated selection between traits (Lande and Arnold, 1983). Correlated selection determines the shape of the selection surface: when traits with genetic variation are not aligned with the directional selection gradients, phenotypic evolution will be slower and distorted, resulting in less direct responses (Phillips and Arnold, 1989; Arnold et al., 2001; Svensson et al., 2021).
Lande’s equation predicts adaptive phenotypic evolution but might fail when indirect selection is important. Indirect selection results from unmeasured traits being genetically correlated with the measured traits or when there is correlated selection between measured and unmeasured traits (Lande and Arnold, 1983; Rausher, 1992). Using Robertson’s Secondary Theorem of Natural Selection (STNS) (Robertson, 1968; Walsh and Lynch, 2018), , trait changes over one episode of selection are accurately predicted because the selection differentials () equal the (additive) genetic covariance () of the trait with relative fitness () and thus the trait’s breeding value change in one generation, regardless of unmeasured traits (Morrissey et al., 2010; Morrissey and Bonnet, 2019). However, distinguishing direct from indirect selection is not possible with Robertson’s STNS, which led Stinchcombe et al., 2014 to propose combining its merits with those of Lande’s equation in a single statistical framework to predict adaptive phenotypic evolution, see also Etterson and Shaw, 2001; Morrissey et al., 2012; Hajduk et al., 2020. Using Lande’s equation retrospectively, one can estimate ‘genetic’ selection gradients as - , where the selection differentials are defined by the difference between the trait measured in an ancestral population (), and the trait of a diverging population () as predicted by the STNS.
Unfortunately, using Lande’s equation with the genetic selection gradients to predict adaptive phenotypic evolution over several generations depends not only on invariable selection gradients but also on the stability of the -matrix. The input of genetic covariances by pleiotropic mutations could in the long-term of mutation-selection balance (time being scaled by the effective population size) be aligned with correlated selection and eventually explain phenotypic differentiation among populations and species (Lande, 1980; Jones et al., 2007; Jones et al., 2014; Chebib and Guillaume, 2017; Houle et al., 2017; Farhadifar et al., 2015; Svensson et al., 2021). However, many studies find more standing genetic variation in natural populations than expected at mutation-selection balance (Walsh and Lynch, 2018; Sella and Barton, 2019). In part, this is because in the initial stages of adaptation selection might not be weak relative to recombination as required by theory (Lande, 1980; Nagylaki, 1992; Turelli and Barton, 1994), in part this is because selection is not constant or uniform in temporally changing or spatially heterogeneous environments (Gomulkiewicz and Houle, 2009; Chevin et al., 2010; de Villemereuil et al., 2020; Walter, 2023). In addition, the -matrix is bound to evolve in the short-term because of selection (Cheverud, 1996; Barton and Turelli, 1987; Turelli, 1988; Shaw et al., 1995), although there is mixed empirical evidence that the -matrix can evolve to align with the orientations of selection (Steppan et al., 2002; Arnold et al., 2008; Chenoweth et al., 2010; Ramakers et al., 2018; Johansson et al., 2021). The reduction in the size of the -matrix due to drift can be predicted because it is inversely proportional to the effective population size (Lande, 1976; Lynch and Hill, 1986; Barton et al., 2017), but the shape of the -matrix will change unpredictably because all populations are finite, and bottlenecks and founder effects are common (Phillips et al., 2001; Phillips and McGuigan, 2006). Hence, besides selection, drift might also impact ongoing phenotypic evolution long before the mutation-selection balance is reached (Whitlock et al., 2002; Mallard et al., 2022).
Here, we ask if Lande’s equation with the genetic selection gradients predicts the direction of adaptive phenotypic evolution for more than one generation. By adaptive phenotypic evolution we mean that multivariate trait responses to indirect or direct selection are correlated with adaptation to a novel and fixed target environment. We focus on the locomotion behavior and body size of the hermaphroditic nematode Caenorhabditis elegans in replicate populations gradually evolving to a high salt (NaCl) concentration in their growth-media for 35 generations and then kept in the constant high salt environment for an extra 15 generations. Our ancestral population was domesticated to constant low salt conditions for 140 generations and contains abundant but neutral standing genetic variation for locomotion behavior (Mallard et al., 2022). Although replicate populations were maintained at high population sizes, our experimental regime exacerbates drift and inbreeding because of a slow rate of environmental change until reaching the target high salt environment (Guzella et al., 2018), and because high salt favors hermaphrodite self-fertilization (Theologidis et al., 2014).
Osmotic pressure from high salt concentration in the growth media shrinks individual body size because of water cell loss (Urso and Lamitina, 2021). For the ancestral population, and in our laboratory conditions, high salt also lowers embryo to adult survival and retards growth until maturity (Theologidis et al., 2014). As hermaphrodites cannot mate with each other, delayed male development results in hermaphrodites reproducing mostly by self-fertilization. During domestication, movement was reduced from that observed among wild isolates (Mallard et al., 2022), and hermaphrodites can further reduce movement in high salt during experimental evolution as males become less frequent and cannot harass them (Barr et al., 2018; Cutter et al., 2019). We further know that the canonical trait of locomotion behavior with the most mutational variance in low salt conditions differs from that with the most standing genetic variation found after domestication (Mallard et al., 2023). Consistent with these observations, several studies have shown that C. elegans mutants insensitive to high salt have specific defects in backward or forward movement, in some of these mutants independently of body size effects (Fujiwara et al., 2002; Swierczek et al., 2011; Zhen and Samuel, 2015). On the other hand, movement can increase during experimental evolution due to more foraging and dwelling, as the bacterial food is not as dense in high salt (Gray et al., 2005). Both foraging, dwelling and mate interactions in C. elegans can be described as a complex collection of distinct behavioral states, which vary in the duration of activity and movement direction (Flavell et al., 2020). All these considerations suggest that multiple traits in locomotion behavior can respond to selection but that it is difficult to a priori define which ones are genetically or environmentally independent. For this reason, we mathematically define individual locomotion behavior in 1-dimensional space by six traits, the six transition rates between stillness, moving forward, and moving backward (Mallard et al., 2022, Mallard et al., 2023). Body size is also analyzed as a seventh trait.
In what follows, we ask whether the ancestral phenotypic plasticity between high and low salt environments is aligned with the ancestral -matrix. We use selection differentials on the seven traits in low salt and high salt environments to predict phenotypic evolution by describing the phenotypic and genetic divergence in high salt. Using Lande’s retrospective equation, we ask if the genetic selection gradients measured in the ancestral population match the phenotypic selection gradients.
Results
Experimental design and analyses
The ancestral population for experimental evolution (A6140) was ultimately derived from a hybrid population of 16 isolates and domestication to a standard laboratory environment in low salt (25 mM NaCl) growth-media conditions for 140 generations (Teotonio et al., 2012; Noble et al., 2017). GA[1,2,4] replicate populations were derived from A6140, with limited founder effects, and independently exposed for 35 generations to a gradual change in salt concentration (8 mM increase each generation) and then kept in constant high salt (305 mM NaCl) for 15 generations. During the experiment, replicate populations were maintained at high population sizes (N=104), and from generation 35 onwards, hermaphrodite self-fertilization became predominant (Theologidis et al., 2014). Using genomic data, effective population sizes have been estimated at Ne=103 in the domestication low salt environment and under partial selfing (Chelo and Teotónio, 2013).
From the ancestral population (A6140), and the three replicate populations at generation 50 (GA[1,2,4]50), inbred lines were derived by self-fertilization of hermaphrodites (Noble et al., 2017; Chelo et al., 2019). Inbred lines were measured for hermaphrodite locomotion behavior and body size at the usual reproduction time of experimental evolution in low and high salt (186, 61, 61, and 42 lines from the ancestral and evolved populations, respectively, with most lines being phenotyped twice; see Methods for details). Six traits defined locomotion behavior: the transition rates between movement states, stillness, forward, and backward (Mallard et al., 2022). For the inbred lines of the ancestral population, we also use self-fertility data at the usual reproduction time of experimental evolution in high salt, as previously reported by Chelo et al., 2019. Finally, we measure the extent of adaptation to high salt conditions using the outbred experimental populations from which the inbred lines were derived. Assays were designed so that grandmaternal and maternal environmental effects were the same for all the samples being compared.
With this data (Figure 1—source data 1, Figure 1—source data 2), we model phenotypic plasticity (mean population differences between environments) and standing genetic variation for locomotion behavior and body size in the ancestral population, the evolution of locomotion behavior and body size, and -matrix evolution in the three replicate populations at generation 50. We estimated phenotypic plasticity and phenotypic differentiation in a multivariate MANOVA model and compared it with a univariate response model similarly defined (see Methods). The MANOVA allows us to test for ancestral phenotypic plasticity and phenotypic divergence while accounting for potentially correlated trait variation. The univariate approach allows us to estimate the inbred lines trait values and to test for the phenotypic divergence of each replicate population relative to the ancestral population but does not account for correlated variation in multivariate phenotypic space. Markov chain Monte Carlo methods were used in a Bayesian framework to estimate the -matrix as half the among-line variance (see Methods) and, for the ancestral population, the -matrix together with the genetic (co)variances between traits and fitness. Table 1 defines the variables employed.
Ancestral population
Phenotypic plasticity between salt environments
Before experimental evolution to high salt, we started by characterizing phenotypic and genetic variation in low and high salt environments in the ancestral domesticated population. We find extensive phenotypic plasticity for locomotion behavior traits and body size (Figure 1, Table 2). Because we employed univariate and multivariate approaches to model phenotypic plasticity (see Methods), we compared the estimated environmental effects using both approaches. We find that univariate and multivariate modeling approaches give similar results (Figure 1—figure supplement 1, Figure 1—source data 3, Figure 1—source data 4). Most transition rates are plastic with salt, except from forward to backward movement states. As expected, body size shrinks in high salt.
-
Figure 1—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig1-data1-v2.txt
-
Figure 1—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig1-data2-v2.pdf
-
Figure 1—source data 3
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig1-data3-v2.pdf
-
Figure 1—source data 4
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig1-data4-v2.pdf
-
Figure 1—source data 5
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig1-data5-v2.pdf
The multivariate approach (MANOVA) allows us to determine the phenotypic dimension of ancestral phenotypic plasticity () that most responds to salt environmental variation (see Methods) as the first eigenvector of the MANOVA SSCP-matrix for the environment factor (Figure 1—source data 5). Transition rates from still to forward or to backward (SF or SB) have opposite loading signs in (Table 3). Body size has the same sign of the transition rates from still and from forward to backward (SB and FB) and the opposite sign relative to the other transition rates.
-
Table 3—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-table3-data1-v2.pdf
-
Table 3—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-table3-data2-v2.pdf
-
Table 3—source data 3
- https://cdn.elifesciences.org/articles/80993/elife-80993-table3-data3-v2.pdf
-
Table 3—source data 4
- https://cdn.elifesciences.org/articles/80993/elife-80993-table3-data4-v2.pdf
-
Table 3—source data 5
- https://cdn.elifesciences.org/articles/80993/elife-80993-table3-data5-v2.pdf
Standing genetic variation
We next partitioned the phenotypic (co)variances among the inbred lines of the ancestor population into the genetic (-matrix) and the residual error (co)variances of transition rates and body size (see Methods). Ancestral -matrices were estimated separately by salt environment, assuming genome-wide homozygosity and no directional non-additive genetic effects. Our estimates are robust to changes in the prior distributions (Figure 2—figure supplement 1).
We find significant genetic variance for all transition rates and body size, in high and low salt environments, except for the transition rate from backward to still in low salt (Figure 2, Figure 2—figure supplement 2). The genetic covariances between transition rates and/or body size show similar values, albeit of lower magnitude in low salt (Figure 2A). The -matrix size is similar between environments (Figure 2B). In both high and low salt environments, transition rates from still to forward or to backward (SF or SB) are negatively correlated with the other transition rates and positively correlated with each other. Body size shows positive genetic covariances with SF and SB.
-
Figure 2—source code 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig2-code1-v2.zip
-
Figure 2—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig2-data1-v2.pdf
-
Figure 2—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig2-data2-v2.pdf
Eigendecomposition of the ancestral -matrix in high or low salt reveals a similar structure between them (Figure 2C, Figure 2—source data 2). The first canonical trait (, Table 3) encompasses most genetic variation (75% in high salt, 63% in low salt). The next two canonical traits contain less genetic variation (between 6% and 20%) but are larger than null expectations (Figure 2—figure supplement 3). Only the first two canonical traits of the -matrix have a similar trait loadings between environments (Figure 2—source data 2).
Ancestral plasticity and genetic variation are not aligned
We compared the main canonical trait of phenotypic plasticity with the canonical traits of the high salt -matrix. Phenotypic plasticity is not aligned with the -matrix (Figure 3). This is because the amount of genetic variance along the dimension of phenotypic plasticity () is not different than that expected by chance (Figure 3A), and also because the angle between and is, if anything, larger than expected by chance (Figure 3B). appears to similarly summarize environmental variation as the third canonical trait from the high salt -matrix (Table 3), a trait that encompasses only 6% of standing genetic variation. As noted before, these differences stem from the association between still-to-forward and still-to-backward transition rates (SF and SB), which are genetically positive and environmentally negative, i.e., they have opposite signs in and , respectively. We suspect that positive associations between SF and SB reveal more dwelling, while negative association more individual foraging (Flavell et al., 2020).
Selection differentials are similar across environments
Selection differentials on transition rates and body size measured in the two salt environments are their genetic covariances with the self-fertility measured in the high salt environment (sk, Table 1). We used estimates of the ancestral -matrices to obtain these selection differentials (see Methods). -matrices and -matrices estimates of genetic (co)variances in locomotion behavior and body size are similar (Figure 4—figure supplement 1). We find genetic variance for self-fertility in high salt (Figure 4, Figure 4—source data 1). We also find that the transition rates from still to forward or from still to backward (SF or SB) measured in high salt have positive genetic covariances with self-fertility, and all other transition rates have negative covariances (Figure 4). Small or no selection differentials exist in low salt transition rates, and only body size shows a clear positive selection differential in both environments. These results are robust to variation in self-fertility (Figure 4—figure supplement 2).
Evolutionary divergence
Adaptation to the high salt environment
Having characterized ancestral standing variation, we describe the divergence of the three replicate populations after 50 generations of evolution. First, we measured the degree of adaptation to high salt by comparing the mean fitness of the GA[1,2,4]50 populations to the ancestral A6140 population in competition experiments against a tester GFP-strain (see Methods; Figure 5—source data 1). We find an increase in the mean fitness of all three replicate populations (Figure 5, Figure 5—source data 2). By generation 50, populations adapted to the high salt environment.
-
Figure 5—source code 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig5-code1-v2.zip
-
Figure 5—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig5-data1-v2.txt
-
Figure 5—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig5-data2-v2.txt
Locomotion traits diverged in low and high salt environments
Concomitant with adaptation, there was phenotypic divergence for the locomotion traits and body size measured in high salt (Figure 6, Table 2). Estimates of phenotypic divergence are robust to multivariate and univariate modeling (Figure 6—figure supplement 1, Figure 6—source data 1). From the univariate models, we find that for each transition rate, at least one replicate GA population differed from the ancestor and that the three replicates showed significant divergence for three transition rates (Figure 6, Figure 6—source data 2). For body size we find that only one replicate populations diverged from the ancestral population. The amount of genetic variance did not limit phenotypic divergence, as the back-to-still and forward-to-still transition rates (BS and FS) diverged while showing the lowest genetic variances in the ancestral population (Figure 2). Eigendecomposition of the MANOVA SSCP-matrix for the population factor further reveals that a single canonical trait explains most phenotypic differentiation between the four populations (; Figure 6—source data 3, Table 3).
-
Figure 6—source code 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig6-code1-v2.zip
-
Figure 6—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig6-data1-v2.pdf
-
Figure 6—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig6-data2-v2.pdf
-
Figure 6—source data 3
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig6-data3-v2.pdf
In the low salt environment, there was less phenotypic divergence than in the high salt environment, with only three out of the six transition rates having at least two replicate populations significantly different from the ancestor (Figure 6—figure supplement 2). Unlike in high salt, body size in low salt showed a marked increase after experimental evolution in all replicate populations.
Genetic variance decreased during evolution
We next characterized genetic divergence by estimating the high salt -matrices of the GA[1,2,4]50 populations and comparing them with the ancestral high salt -matrix. We did not model the evolution of the -matrix in the low salt environment. This analysis shows that the size of the high salt -matrix was reduced during experimental evolution, independently of the replicate population (Figure 7, Figure 7—source data 1). However, we continue to find that most genetic variances for the GA populations differ from null expectations (Figure 7—figure supplement 1). Eigendecomposition of the GA -matrices indicates that 3–5 canonical traits differ from null expectations (Figure 7—figure supplement 2). Furthermore, evolved populations continued to have significant genetic variances in the three canonical traits of ancestral standing genetic variation (Figure 7—figure supplement 3).
Random skewers analysis shows that when projecting the four high salt -matrices along 1000 random phenotypic directions, 250 of them showed a significant difference between the ancestral and at least one of the evolved populations (see Methods). Using these 250 vectors we build an -matrix of genetic divergence between the four populations. Eigendecomposition of the -matrix then revealed that a single canonical trait explains most divergence (, Table 1). In this canonical trait at least 2 of the 3 evolved populations showed a reduction in variance (Figure 7C). An alternative eigentensor analysis to detect genetic divergence among the 4 populations confirms the random skewers analysis (see results in our GitHub appendix).
Divergence along ‘genetic lines of least resistance’
We asked whether phenotypic divergence () and genetic divergence () occurred along the dimensions of most ancestral genetic variation (). For this analysis, we calculated the angle (, see Equation 3 in Methods) and the proportion of overlap (, see Equation 5) between these canonical traits. Table 3 summarizes the main canonical traits of ancestral standing variation, and of phenotypic and genetic divergence.
Most of the genetic variance of the ancestral high salt -matrix along is higher than expected by chance (Figure 8A). This is because and are aligned and their angle is very small when compared with other canonical traits of ancestral standing genetic variation, or with a null expectations (Figure 8B). We do not find a significant proportion of genetic variance of the high salt -matrix along (Figure 8C), but find a small angle between and (Figure 8D), which is indicative of a good alignment between these canonical traits.
Indirect selection and predicting phenotypic evolution
Expected and observed responses to selection are aligned in high salt
Using Lande’s retrospective equation, we compared the genetic selection gradients obtained with selection differentials on traits measured in high salt (; see Methods, Equation 6) to the phenotypic selection gradients obtained with the observed responses in high salt after 50 generations (; see Methods, equation 7). The ancestral population’s high salt -matrix was assumed stable during experimental evolution. Credible intervals were obtained, however, by sampling the -matrix from its posterior distribution, with fixed ancestral selection differentials (sk; Figure 4), or fixed observed phenotypic divergence for each of the three replicate populations (; Figure 6 for high salt, Figure 6—figure supplement 2 for low salt).
We find that phenotypic selection gradients for all traits are highly heterogeneous because higher phenotypic divergence has an outsize effect on the mean and error estimates. Nonetheless, the phenotypic selection gradients overlap with the corresponding genetic selection gradients for at least one replicate population (Figure 9A; Figure 9—source data 1, Figure 9—source data 2). Evidence for a lack of overlap between selection gradients in two replicate populations is found for the transition rate backward-to-forward and for body size. Only for the transition rate SB is there an overlap of selection gradients for all three replicate populations.
-
Figure 9—source code 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig9-code1-v2.zip
-
Figure 9—source data 1
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig9-data1-v2.pdf
-
Figure 9—source data 2
- https://cdn.elifesciences.org/articles/80993/elife-80993-fig9-data2-v2.pdf
Given replicate heterogeneity, testing whether selection theory predicts the direction of phenotypic evolution is possibly best estimated as the angle between expected and observed responses to selection. For traits measured in high salt, we find that the angle between selection differentials and observed phenotypic divergence is low for all replicates (Figure 9B). In contrast, when traits are measured in low salt, expected and observed responses are not aligned. Similarly, whether theory predicts the magnitude of phenotypic evolution can be estimated as the ratio between expected and observed responses. For all traits measured in high salt, and across replicate populations, observed phenotypic divergence is on average, across traits and replicates, 3.5 times the selection differentials (Figure 9C). For the traits measured in low salt, predictions about the magnitude of divergence further degenerate (not shown).
Indirect versus direct selection
There is evidence of direct selection on SB and backward-to-forward (BF), as well as on body size (Figure 9A, Figure 9—source data 1). We find a positive genetic selection gradient for SB and a negative genetic selection gradient for BF. There is also a positive selection gradient for body size, opposite in sign to its association with other traits in the canonical traits of standing genetic variation, and in the canonical traits of phenotypic and genetic divergence (Table 3). Traits with low genetic variance, FS and BS (Figure 2—figure supplement 2), do not appear to bias the genetic selection gradient estimates (Figure 9—figure supplement 1). However, none of the genetic selection gradient estimates differ from zero when sampling both the G-matrix and the selection differentials from their respective posterior distributions to obtain credible intervals (Figure 9—figure supplement 2).
Discussion
Modeling approaches for predicting adaptive phenotypic evolution before mutation-selection balance is reached are based on Lande’s equation (see Introduction). In Lande’s equation, the trait change over one generation equals the ancestral -matrix times the directional phenotypic selection gradients, but might not be accurate in the presence of indirect selection. This can be remediated by replacing the phenotypic selection gradients with the genetic selection gradients obtained after measuring selection differentials in the ancestral population. For predicting phenotypic evolution over several generations, however, one must assume invariable selection gradients and that the -matrix is stable despite selection and drift. We sought to test the selection theory by finding whether we could predict phenotypic evolution for 50 generations of experimental evolution.
We followed seven traits with different environmental and genetic dependencies, the six transition rates between movement states and body size (Table 3). Our ancestral population was adapted to the low salt conditions (Chelo and Teotónio, 2013; Theologidis et al., 2014), before challenging three replicate populations to a gradual increase in the salt concentration in the growth media for 35 generations and 15 extra generations in high salt. Fifty generations of experimental evolution led to adaptation (Figure 5) and to phenotypic and genetic divergence (Figure 6, Figure 7). Body size measured in high salt did not consistently evolve among replicate populations but individual movement increased. This is because the transition rates from the still state have increased while those to the still state have decreased. Adaptive phenotypic divergence followed the direction of the canonical trait with more ancestral standing genetic variation (Figure 8, Table 3), and therefore, we could predict phenotypic evolution, though only its direction and when the component traits of the multivariate phenotype were measured in the high salt environment (Figure 9). When considering the component traits of the multivariate phenotype individually, and due to replicate population heterogeneity, we could confidently predict the evolution of only one of the seven traits followed (SB). We could not predict the magnitude of evolution for any individual trait. These findings are relatively unique because we described -matrix evolution and measured the ancestral selection differentials to predict phenotypic evolution, but they add to the results of a growing number of experimental studies testing Lande’s equation across tens of generations. For example, a recent re-analysis of up to 60 generations in constant and homogeneous environments, for five wing traits in Drosophila melanogaster, showed an alignment between the main canonical trait of genetic variation in the evolved populations with adaptive phenotypic divergence (Yeaman et al., 2010; Walter, 2023).
As the canonical trait explaining most variation in the ancestral population, , also the second and third canonical traits differ from null expectations (, Figure 2, Figure 2—figure supplement 3), and remain so after evolution (Figure 7—figure supplement 3), despite potential variance inflation problems due to the MCMC methods we employed (Morrissey et al., 2012; Sztepanacz and Blows, 2017b). For these ancestral canonical traits, selection must have been responsible for the observed loss of genetic variance, particularly after generation 35 (Figure 7). Assuming an infinitesimal model of trait inheritance, drift is expected to lead to a loss of genetic variance by (1-1/2) at each generation (Barton et al., 2017). Even in the unrealistic situation of complete selfing during the experiment, and considering effective population sizes on the order of 1000 (Chelo and Teotónio, 2013), less than 5% was expected to be lost by genetic drift by generation 50, values that were not observed (only about half of the genetic variance was lost, Figure 7B). Supporting loss of variation by selection we before showed that allelic diversity at neutral single-nucleotide polymorphisms is reduced relative to ancestral levels by 5% by generation 22 (Theologidis et al., 2014), and to 20% only by generation 50 (Chelo et al., 2019). We did not, however, test whether the direction of phenotypic divergence occurred along the ancestral traits because they together explain 15% of the variation and there is a poor statistical power to do so. Furthermore, our previous work suggests that variation in the ancestral might have been lost by drift during the first 35 generations of the experiment such that they were perhaps of little consequence later on when populations reached the high salt environment (Guzella et al., 2018). Few studies have demonstrated that canonical traits of little genetic variance can influence selection responses (Kirkpatrick, 2009; Blows and McGuigan, 2015; Sztepanacz and Blows, 2017a). In one of the few examples, Hine et al., 2014 showed that low-variance canonical traits of eight cuticular hydrocarbons in Drosophila serrata respond to artificial selection during six generations, though inconsistently among replicate populations. In our experiment, we suspect that showing that canonical traits with a small amount of genetic variation impact adaptive phenotypic divergence will require finding the quantitative trait loci (QTL) responsible for their expression (Svensson et al., 2021; Kelly, 2009). If these low-variance canonical traits influence adaptation then allele frequency dynamics at the relevant QTL because of selection might be detected when comparing genomic data between ancestral and derived populations (Long et al., 2015; Barghi et al., 2020).
Our findings also highlight the relationship between phenotypic plasticity and adaptation to novel environments (Ghalambor et al., 2007; Pfennig et al., 2010; Teotónio et al., 2009; Draghi and Whitlock, 2012; Noble et al., 2019). While the discussion has been on showing that population persistence is more likely if plasticity is aligned with the direction of selection (Price et al., 2003; Lande, 2009; Chevin et al., 2010), our results show that plasticity only reveals the topography of the adaptive landscape. In high salt conditions, populations move away from the ancestral phenotypic optimum (Figure 1), with an associated fitness cost (Theologidis et al., 2014). Adaptation to the high salt target environment after generation 35 presumably involved recovering a phenotype similar to that of the ancestral population that alleviated this fitness cost (Table 3). In particular, high salt in the ancestral population reduces body size and SB transition rates while increasing SF transition rates. Symmetrically, selection favors increased body size, increased SB, and decreased SF. However, because all three traits show positive genetic covariances with each other (Figure 2), even if plasticity is oriented with selection (but of the opposite sign), phenotypic evolution is constrained by a lack of genetic variation in the appropriate canonical trait. The ancestral population had genetic variation in the direction of selection (the canonical trait g3, Table 3), but as argued above it was probably lost during gradual salt evolution because of drift such that when reaching high salt populations could not have further responded to selection (Matuszewski et al., 2015; Guzella et al., 2018). Future evolution in the direction of the ancestral multivariate phenotypic optimum, or close to it, should then be conditional on the appearance of de novo pleiotropic mutations. Assuming that mutational covariances do not vary with the environment, it is unclear that there can be much further phenotypic evolution as elsewhere we characterized mutational covariances and did not find any in the direction of selection (see Table 3 and Mallard et al., 2023).
Ancestral phenotypic plasticity can thus be considered ‘non-adaptive’ (Ghalambor et al., 2007). Hence, it is unsurprising that we could not predict low salt phenotypic evolution. This is not explained because of a general lack of genetic variance in locomotion traits and body size in the ancestral population (Figure 2) but because their selection differentials were small or did not differ from zero (Figure 4). An exception is body size. Body size is reduced by high salt conditions (Figure 1), and there is probably direct selection for increased body size in high salt (Figure 9). However, there was little body size evolution in the high salt environment (Figure 6). In contrast, there was evolution of increased body size when it was measured in the low salt environment, as expected from the positive selection differentials in the ancestor population (Figure 4). Body size measured in low salt could only have evolved because of indirect selection by being genetically correlated with high salt body size or some other traits expressed in high salt. Hence, the evolution of plasticity in body size might be predictable although the evolution of plasticity in the multivariate phenotype is not. The fact that body size measured in low salt evolved further reveals that the evolved populations were in a region of the adaptive landscape that was not readily accessible to the ancestral population even if it had been domesticated to low salt conditions for 140 generations. In other words, body size evolution in low salt supports the long-held hypothesis that phenotypic plasticity facilitated the evolution of novelty (Wagner and Lynch, 2010; Moczek et al., 2011; Levis et al., 2018).
Robertson’s covariance [] was originally found to describe an episode of selection (the selection differentials, ) and later to predict adaptive phenotypic evolution over one generation () as the secondary theorem of natural selection (Walsh and Lynch, 2018; Hajduk et al., 2020). Applications so far have been mostly limited to explaining the evolution of individual traits and when time-series of trait and fitness are concurrently obtained with pedigrees such that estimates of trait breeding values and genetic selection gradients can be updated every few generations (Stinchcombe et al., 2014; Morrissey et al., 2012; Hajduk et al., 2020; Hadfield, 2010). One problem that has received particular attention has been to ask whether indirect selection or phenotypic plasticity and robustness explain phenotypic stasis, despite trait heritability and a significant phenotypic selection gradient on the observed trait (Merilä et al., 2001; Kruuk et al., 2002). For example, Biquet et al., 2022 followed a blue tit population for more than 40 years, in which egg laying has changed to earlier spring dates, presumably because of climate change. Despite significant heritability of egg laying and directional selection for earlier dates, modeling the breeding values did not reveal any temporal trend, consistent with a lack of genetic covariance between egg laying date and fitness. There was thus no genetic divergence for laying date, which Biquet et al., 2022 could attribute to phenotypic robustness and the stochastic nature of individual development to maturity. Conversely, using a similar approach, Bonnet et al., 2017 found that evolution towards earlier parturition dates in a red deer population could be predicted and was consistent with the estimated change in breeding values. In a study following two traits, selection due to human harvesting of a prey species of wild salmon (for feeding domestic salmon) explained divergence in early maturity and small body sizes, despite directional selection for increased body size at maturity because of fishing (Czorlich et al., 2022). Our results suggest that when pedigrees are difficult to obtain, predicting the direction of adaptive multivariate phenotypic evolution for tens of generations may be possible without updating estimates of selection gradients and the -matrix every few generations. This is because although the environment changed for 35 generations in our experiment directional selection was maintained and only the size of the -matrix was reduced.
In sum, we have shown that using Lande’s equation with genetic selection gradients is valid to predict the direction of phenotypic evolution in a new environment, after a gradual environmental change for 35 generations and 15 generations in the new environment. However, selection theory not necessarily predict the direction or magnitude of evolutionary change in all the component traits of the multivariate phenotype, especially if the traits are not measured in the new environment. This is because there are variable and complex genetic and environmental dependencies between individual traits. There are few experimental tests of selection theory such as ours. Therefore, more will be needed to generalize our results to natural populations, particularly those challenged by changing and heterogeneous environments.
Materials and methods
Experimental populations and environmental conditions
Request a detailed protocolThe ancestral population is named A6140, where ‘A’ stands for androdioecious, ‘6’ for replicate six, and ‘140’ for the number of generations of domestication to a standard laboratory environment (Teotónio et al., 2017). A6140 resulted from the hybridization of 16 founder wild strains during 33 generations followed by 140 generations characterized by 4 day discrete and non-overlapping life-cycles at N=104census sizes and =103effective population sizes (Teotonio et al., 2012; Chelo and Teotónio, 2013; Theologidis et al., 2014). Our standard laboratory environment involves populations being maintained in 10 × 9 cm Petri dishes NGM-lite agar media containing 25 mM NaCl and a homogenous lawn of E. coli HT115 that served as food from the L1 larval stage until reproduction. Each Petri dish contains 1000 individuals which are mixed during reproduction, with embryos being collected and synchronized at the L1 larval stage to start a new generation.
We report the evolution of locomotion behavior and body size in three independent replicate populations (named GA[1,2,4] populations: ‘G’ for gradual, ‘A’ for androdioecious, ‘#’ for replicate number). They were derived from splitting into three a single pool of at least 104 individuals sampled from the A6140 population. GA populations were maintained in the same conditions as during domestication except that the NGM-lite media was supplemented with 8 mM of NaCl at each generation for 35 generations and then kept at constant 305 mM NaCl for an additional 15 generations. Details about the derivation of the GA populations can be found in Theologidis et al., 2014. We refer to the NGM-lite 305 mM NaCl environment as the ‘high salt’ target environment, while the domestication 25 mM NaCl environment as the ‘low salt’ environment.
C. elegans is an androdioecious nematode, where hermaphrodites can self but outcross only when mated with males. Natural populations are depauperate of genetic diversity and males are rare due to a long history of selfing, selective sweeps, and background selection (Andersen et al., 2012; Rockman et al., 2010). Under the domestication environment, however, outcrossing is readily maintained at frequencies between 60% and 100% (Teotonio et al., 2012, Mallard et al., 2022). In GA populations outcrossing is maintained at close to 100% for 35 generations, reduced to about 30% in GA1, 14% in GA2, and 5% in GA4, by generation 50 of the experiment (Theologidis et al., 2014).
As reported before, we derived inbred lines by selfing single hermaphrodites from the ancestor (A6140) and the three replicate populations at generation 50 (GA[1,2,4]50) for a minimum of 10 generations (Noble et al., 2017; Noble et al., 2021). Male frequency in the inbred lines is low, on the order of the mutation rate for the non-disjunction of the X-chromosome (Teotónio et al., 2006) – sex-determination is chromosomal with hermaphrodites XX and males X –.
Populations and inbred lines were cryogenically stored (Stiernagle, 1999), allowing for contemporaneous measurements of ancestral and evolved outbred populations and their inbred lines. Grandmaternal and maternal environmental effects are common to the samples being measured (Teotónio et al., 2017).
Adaptation to high salt
Request a detailed protocolWe measured the increase in mean relative fitness among the ancestral population (A6140) and evolved populations at generation 50 (GA[1,2,4]50) using pairwise competition experiments between them and a tester line (Teotónio et al., 2017). As a tester, we employed an inbred line (EEV1402) derived by selfing from the A6140 population, and that expressed a green-fluorescent-protein (GFP) morphological marker (Chelo et al., 2013). For the assays, we revived the four populations and the tester line (>1000 individuals each) and let individuals reproduce and starve for 10 days. Starved individuals were then seeded on fresh plates with food at a density of 1000 L1 larvae in low salt. We grew them for two complete generations in high salt, except the GFP tester which was only grown in high salt for one generation. At the third generation, we seeded 500 L1 larvae of the GFP tester line together with 500 L1 larvae of 1 of 4 experimental populations in high salt. For A6140, we seeded 15 plates (technical replicates), for GA150 4 plates, for GA250 four plates, and for GA450 five plates. In each of these plates, 72 hr after L1 seeding, individuals were subject to the ‘bleach/hatch-off’ protocol, the standard of our life-cycle, to recover live embryos and, 24 hr later, synchronized L1 larvae. We scored an average of 169 larvae for GFP expression in each technical replicate.
The relative proportion of non-GFP to GFP measures the relative fitness of the experimental populations to the tester after one generation of competition (Teotónio et al., 2017). To analyze this data, we used a generalized linear model in R (R Development Core Team, 2018), testing for the evolution of the ratio non-GFP/GFP, assuming a binomial error distribution (‘quasibinomial’ family option) and allowing for overdispersion of the data. Post-hoc pairwise comparisons were performed between the ancestral and the evolved populations with Tukey tests using the glht function in the multcomp package in R (Hothorn et al., 2008).
Locomotion behavior
Request a detailed protocolInbred lines were thawed from frozen stocks on 9 cm Petri plates and grown until exhaustion of food. This occurred 2–3 generations after thawing, after which individuals were washed, adults removed by centrifugation, and three plates per line seeded with 1000 larvae at mixed larval stages. Samples were then maintained in the standard domestication environment for two complete generations. At the assay generation (generation 4–6 generations post-thaw), starvation-synchronized L1 larvae were seeded in low and high salt. Adults were phenotyped for locomotion behavior 72 hr later at their usual reproduction time in one 9 cm plate (technical replicate). At the beginning of each assay we measured ambient temperature (T) and humidity (H) in the imaging room.
Given the number of lines to phenotype, we repeated the above protocol several times over several years, with each repetition defining a statistical ‘block’ on a given day. In total, we phenotyped 186 lines from the A6140 population and 61, 61, and 42 lines from each of the GA[1,2,4] populations, respectively, with most lines being phenotyped twice and always in separate blocks (average of 1.9 in low salt, and of 2 in high salt).
We imaged adult hermaphrodites using the Multi-Worm Tracker [version 1.3.0; Swierczek et al., 2011 and used the materials and protocols of Mallard et al., 2022]. Each movie contains about 1000 tracks of hermaphrodites (called objects) with a mean duration of about 1 min. Standardized to a common frame rate (4 Hz), we filtered and extracted the number and persistence of tracked objects per movie and assigned movement states across consecutive frames as forward, still or backward (assuming forward as the dominant direction of movement). Mean object density (D) per movie was also retrieved to be used as a covariate in modeling.
Locomotion behavior in 1-dimensional space is described by the transition rates between still (S), forward (F) and backward (B), plus the self-transition rates. Modeling is detailed in Mallard et al., 2022. Transition rates between movement states are assumed to follow a continuous time Markov process. The Markov process is a stochastic process modeling changes in movement state as a matrix Q. In our data, the Markovian memoryless assumption is only marginally violated (Mallard et al., 2022). The elements in Q, noted , are the transition rates from state i to state j (off-diagonal elements for , and with ). This definition constrains self-transition rates (diagonal elements) to be of the opposite sign to the sum of the two transition rates leaving the relevant movement state:
This ensures that the probability of leaving a given movement state towards any other state during a waiting time is one minus the probability of remaining in the same state (see Mallard et al., 2022 for a more detailed explanation). Therefore, only the six transition rates between movement states are mathematically independent and we thus ignore self-transition rates.
For estimation, we used log-likelihood models as defined in Mallard et al., 2022 and specified them with the msm package (Jackson, 2011) in RStan [Stan Development Team, 2018, R version 3.3.2, RStan version 2.15.1]. Because , all analyses were performed on the natural log scale to ensure normality. We used multi-log normal prior distributions with the mean transition rate and a coefficient of variation . We retained the means of the posterior distributions as the per-plate transition rates for all the subsequent analyses.
Body size
Request a detailed protocolWe included the measurements of body size obtained from the Multi-Worm Tracker movies as a seventh trait. Movie frames were sampled only for forward tracked-objects to minimize posture variation. We then extracted the per-track object mean area (Swierczek et al., 2011). These values were summarized as the per-plate median of all track mean values. We then re-scaled these measurements so that the averaged phenotypic variance in each environment is roughly similar to the average transition rates phenotypic variance. This was done by multiplying the body size by a factor of 50. We chose this procedure rather than dividing all the phenotypic values by their mean, cf. Houle et al., 2011, because our transition rates’ means are close to zero while spanning both negative and positive values. Dividing these transition rates by their means would lead to an artificial increase in phenotypic variance.
Self-fertility
Request a detailed protocolTo estimate selection differentials, we used previously-published data on hermaphrodite self-fertility of the A6140 inbred lines in high salt (Chelo et al., 2019). Self-fertility was measured under environmental conditions that closely followed those of experimental evolution. An average of 42 hermaphrodites were scored for self-fertility per inbred line (minimum 22 and maximum 85 individuals). Self-fertility includes the fecundity of hermaphrodites at the usual time of reproduction by selfing and the viability of their progeny until the L1 larval stage. The log-transformed, covariate-adjusted self-fertility values (best linear unbiased prediction estimates, BLUPs) for each inbred line were downloaded from Chelo et al., 2019, exponentiated, and divided by the mean to obtain a proxy for relative fitness (noted ; Table 1).
Phenotypic plasticity and phenotypic divergence
Request a detailed protocolWe used a multivariate analysis of variance (MANOVA) to model ancestral phenotypic plasticity and the divergence of locomotion behavior and body size. The six transition rates and body sizes were fitted as a multivariate response variable, with fixed effects of temperature and humidity at the time of movie recording and the log of object density in each Petri plate. These three environmental variables were centered and standardized before the analysis to mean=0 and sd=1. We further modeled a fixed effect of block and a fixed effect of year accounting for when the different lines were measured. The main factors of interest were the fixed effects of the salt environment and the fixed effects of evolution, together with their interaction; the last factor with four population ID levels (A6140, GA[1,2,4]50). The residual error was assumed to follow a multivariate normal distribution. We used the manova function in the stats package in R for computation (R Development Core Team, 2018), with Wilks tests are being used for the significance.
From the MANOVA results, we extracted the Sums-of-Squares and Cross-Products (SSCP) matrices for the fixed effects of environment and population and eigendecomposed these matrices to describe the orthogonal canonical traits maximizing phenotypic variation in each (Walter, 2023). For the SSCP-matrix of the environment, the first eigenvector is the dimension containing the most phenotypic plasticity in the ancestral population and is here named . For the SSCP-matrix of evolution, the first eigenvector is the dimension of divergence among the four populations in high salt and is here called are the eigenvalues measuring the variation explained by each eigenvector. Estimated mean-least square divergence per replicate population is here called , with being the environment.
Additionally, we modeled the traits individually using linear mixed-effects models to estimate the best linear unbiased predictions (BLUPs) of transition rates and body size per inbred line (used only for visualization purposes in the figures). This univariate approach allowed testing the divergence of transition rates for each replicate GA population from the ancestral A6140 population. This univariate model was similarly formulated as the MANOVA, except the block was included as having random effects. For model fitting, we employed the lme4 package in (Bates et al., 2015) in R. Post-hoc pairwise contrasts employed Tukey tests with the emmeans package (Lenth, 2021).
G-matrices and genetic divergence
Request a detailed protocolUsing the same model, we estimated the -matrices of the ancestral population A6140 and the three evolved replicate populations GA[1,2,4] separately for the traits measured in the low and high salt environments.
The six transition rates and body size were fitted as a multivariate response variable column-vector in the model:
where are the intercepts and are the environmental fixed effects of temperature (T), humidity (H), and log density (D). We denote to simplify notation of the product (×) among the environmental variables (fitting all three variables as fixed effects, the three two ways interactions, and the three-way interaction for a total of seven fixed effects). was defined as the fixed effect of year when the assays were conducted, and the random effects of line and block identity, respectively. defines the residual error.
The -matrix is half the line identity (co)variance matrix (), as we have measured homozygous diploid inbred lines and assume codominance. As estimated here, the broad-sense -matrix should be an adequate surrogate for the narrow-sense -matrix. This is because there is no inbreeding depression for self-fertility in high salt due to the self-fertilization of hermaphrodites from the experimental outbred populations (Chelo et al., 2019), and because, at least in low salt, we failed to detect average (genome-wide) directional dominance or epistasis when comparing the means of transition rates in the outbred populations with those among the inbred lines (Mallard et al., 2022).
Models were fit with the R package MCMCglmm (Hadfield, 2010). We used improper flat priors (nu=0). Model convergence was verified by visual inspection of the posterior distributions and an autocorrelation below 0.05. 100,000 burn-in iterations were done with a thinning interval of 2000 and over 2 million MCMC iterations. The A6140 -matrix in high salt was estimated using different prior distributions, chosen among the literature as the most representative including parameter-expanded priors (see Figure 2—figure supplement 1).
Because the variance estimates resulting from the MCMCglmm models are positive definite, null expectations for the -matrices were obtained by randomizing 1000 times the phenotypic data set. Randomization was done by shuffling inbred line and block identities and refitting the model at each iteration (Equation 2). We then computed the posterior mode for each of the 1000 models to construct a null distribution of genetic variances.
Eigendecomposition of each -matrix was done in R as above for phenotypic (co)variances. We define the main canonical dimension of genetic variation as the first eigenvector of the A6140 -matrix (with its eigenvalue). We calculate the angle between the two in high and low salt as the mean of the estimated posterior distribution modes (this angle is defined in the next section, see Equation 3).
We used the random skewers method described by Aguirre et al., 2014; Hine et al., 2009 to describe the genetic divergence during experimental evolution. In this method, random vectors are projected through the four -matrices to estimate the genetic variance in all phenotypic directions (Equation 5, see below). Using the -matrix posterior distributions, we tested for significant differences in genetic variance between matrices for each random vector. The vectors that showed a significant difference (i.e. no overlap between the 95% CI of the two matrices projected variance) were retained to construct an -matrix with the (co)variances of differentiation. The eigendecomposition of the -matrix then describes the canonical traits of genetic differentiation among the -matrices. The first eigenvector the -matrix is here called the vector of genetic divergence (, Table 1) because the A6140 ancestral population drives most differentiation. An alternative -matrix differentiation analysis can be done with the eigentensor approach (Aguirre et al., 2014). Eigentensor analysis of the four A6140, GA[1,2,4]50 -matrices in high salt gave similar results (see methods and results in the GitHub appendix).
Selection differentials
Request a detailed protocolFor the ancestral population A6140, we also computed the -matrices as defined in Stinchcombe et al., 2014, which is the -matrix of the 6 traits of locomotion behavior and body size expanded to include self-fertility. The last column-vector entries of the -matrix are thus the covariances between traits and relative fitness, the selection differentials (sk).
Different individuals in separate assays were measured for self-fertility and transition rates/body size. To assess for a statistical bias on selection differential estimates when using self-fertility BLUP estimates (Hadfield et al., 2010), we generated 500 ancestral -matrices in high salt with within-line self-fertility variability across the replicated measurements of transition rates. For each transition rate measurement (one per Petri dish, see above), one inbred line self-fertility value was sampled from a normal distribution using the line’s mean and, as standard deviation, the standard error of the mean multiplied by . In each line mean self-fertility was calculated from multiple individuals (at least 22 and up to 85) and on average there are two transition rate values per line. Our protocol thus mimics a random split of self-fertility into two groups of identical size. The -matrix is stable to within-line self-fertility variation and subsequent analysis was done with the initial -matrix estimates. We also ensured that the -matrix contained in the -matrix is similar to the one computed above for the ancestral population.
Phenotypic and genetic alignments
Request a detailed protocolWe used the metrics introduced by Noble et al., 2019 to compare the alignment of ancestral standing genetic variation with the first canonical dimension of phenotypic plasticity (), or with the first canonical dimensions of adaptive phenotypic () or genetic () divergence. The first metric is the angle between two vectors. The angle between the i-th eigenvector of the A6140 -matrix, gi, and is defined as:
As both gi and are eigenvectors of the -matrix, values between 90° and 180° were transformed so that always remains between 0° and 90° (=180°-, results from using instead of gi in Equation 3). Angles comparing the alignment of the ancestral gi with the axis of phenotypic and genetic divergence were calculated, by replacing in Equation 3 with and , respectively.
For each angle, we sampled the posterior distribution of the A6140 -matrix to create a credible interval. and were obtained as the first eigenvectors of the SSCP matrices from the MANOVA model, as described above. The null expectation for is calculated as the angle between 1000 pairs of random vectors sampled from a uniform distribution .
The second metric computes the proportion of ancestral genetic variance along the main canonical trait of ancestral phenotypic plasticity:
where is replaced by when computing the proportion of ancestral genetic variance in the main canonical trait of phenotypic divergence in high salt, or by when relative to the main canonical trait of genetic divergence in high salt.
is the ratio between the amount of genetic variance in that maximizes plasticity, phenotypic or genetic divergence, over the maximum possible amount of genetic variance in any phenotypic dimension (, the first eigenvalue of the -matrix):
values are comprised between 0 (no genetic variance along the plasticity/divergence canonical traits) and 1 (when the plasticity/divergence canonical traits contain all the genetic variance in ). The null distributions for were obtained by randomizing 1,000 ancestral -matrices through shuffling inbred line and assay block identities.
Selection differentials and gradients
Request a detailed protocolSelection differentials were estimated above with the ancestral matrix in high or low salt as the genetic covariance between transition rates and body size with self-fertility in high salt (sk). Comparing observed and expected responses to selection was done by estimating directional selection gradients using Lande’s retrospective equation, equation 9 in Lande, 1979. This is unlike Stinchcombe et al., 2014 or Hajduk et al., 2020, where phenotypic selection gradients were obtained by regression of fitness onto the traits, following Lande and Arnold, 1983. In our case, genetic selection gradients on each transition rate were defined as:
and the phenotypic selection gradients as:
The -matrix of the ancestral population was assumed constant during experimental evolution. The credible intervals of both selection gradients were estimated by sampling the posterior distribution of the -matrix, assuming fixed high salt sk in Equation 6 or fixed for each replicate population in equation 7. We have also obtained credible intervals for by sampling the -matrix and the posterior distribution of sk.
Whether selection theory predicts the direction of phenotypic evolution amounts to an alignment between expected and observed phenotypic divergence. We thus calculated the angle (as above, ) between the selection differentials on transition rates and body size in high or low salt (sk), with the observed phenotypic divergence in high or low salt (). The null expectations for the angle were obtained by calculating the angles between 1000 pairs of random vectors sampled from a uniform distribution . Similarly, whether selection theory predicts the magnitude of phenotypic evolution in high salt can be calculated as the ratio of observed phenotypic divergence over selection differentials (/sk). We sampled the posterior distribution of -matrix for these comparisons to obtain credible intervals.
Contrasts between posterior distributions
Request a detailed protocolThe ‘significance’ of the posterior mode estimates are based on its overlap with the posterior null distribution of the posterior modes (Walter et al., 2018). For all comparisons of posterior distributions, significance can be inferred when their 83% credible intervals do not overlap (Austin and Hux, 2002), assuming homoscedasticity.
Archiving
Request a detailed protocolSelf-fertility data has been previously published by Chelo et al., 2019, and locomotion behavior data in low salt for the ancestral population in Mallard et al., 2022. New data (adaptation, locomotion behavior and body size in high salt), R code, and modeling results are in our GitHub repository and will be archived in a public repository upon publication.
Appendix 1
Data availability
New data, R code for analysis and modeling results is freely accessible and can be found at https://github.com/ExpEvolWormLab/Mallard_Robertson (copy archived at Mallard and Teotonio, 2023).
References
-
Understanding the evolution and stability of the g-matrixEvolution; International Journal of Organic Evolution 62:2451–2461.https://doi.org/10.1111/j.1558-5646.2008.00472.x
-
Phenotypic evolution: the ongoing synthesisThe American Naturalist 183:729–746.https://doi.org/10.1086/675304
-
A brief note on overlapping confidence intervalsJournal of Vascular Surgery 36:194–195.https://doi.org/10.1067/mva.2002.125015
-
Polygenic adaptation: a unifying framework to understand positive selectionNature Reviews. Genetics 21:769–781.https://doi.org/10.1038/s41576-020-0250-z
-
The infinitesimal model: Definition, derivation, and implicationsTheoretical Population Biology 118:50–73.https://doi.org/10.1016/j.tpb.2017.06.001
-
Fitting linear mixed-effects models using lme4Journal of Statistical Software 67:1–48.https://doi.org/10.18637/jss.v067.i01
-
Phenotypic plasticity drives phenological changes in a Mediterranean blue tit populationJournal of Evolutionary Biology 35:347–359.https://doi.org/10.1111/jeb.13950
-
What affects the predictability of evolutionary constraints using a G-matrix? The relative effects of modular pleiotropy and mutational correlationEvolution; International Journal of Organic Evolution 71:2298–2312.https://doi.org/10.1111/evo.13320
-
An experimental test on the probability of extinction of new genetic variantsNature Communications 4:2417.https://doi.org/10.1038/ncomms3417
-
The opportunity for balancing selection in experimental populations of Caenorhabditis elegansEvolution; International Journal of Organic Evolution 67:142–156.https://doi.org/10.1111/j.1558-5646.2012.01744.x
-
Partial selfing can reduce genetic loads while maintaining diversity during experimental evolutionG3: Genes, Genomes, Genetics 9:2811–2821.https://doi.org/10.1534/g3.119.400239
-
The contribution of selection and genetic constraints to phenotypic divergenceThe American Naturalist 175:186–196.https://doi.org/10.1086/649594
-
Developmental integration and the evolution of pleiotropyAmerican Zoologist 36:44–50.https://doi.org/10.1093/icb/36.1.44
-
BookThe Genetical Theory of Natural SelectionOxford University Press.https://doi.org/10.5962/bhl.title.27468
-
Demographic and genetic constraints on evolutionThe American Naturalist 174:E218–E229.https://doi.org/10.1086/645086
-
The spandrels of San Marco and the Panglossian paradigm: a critique of the adaptationist programmeProceedings of the Royal Society of London. Series B, Biological Sciences 205:581–598.https://doi.org/10.1098/rspb.1979.0086
-
MCMC methods for multi-response generalized linear mixed models: The mcmcglmm r packageJournal of Statistical Software 33:1–22.https://doi.org/10.18637/jss.v033.i02
-
The misuse of BLUP in ecology and evolutionThe American Naturalist 175:116–125.https://doi.org/10.1086/648604
-
The “algebra of evolution”: the Robertson-Price identity and viability selection for body mass in a wild bird populationPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190359.https://doi.org/10.1098/rstb.2019.0359
-
Characterizing the evolution of genetic variance using genetic covariance tensorsPhilosophical Transactions of the Royal Society B 364:1567–1578.https://doi.org/10.1098/rstb.2008.0313
-
Evolutionary constraints in high-dimensional trait setsThe American Naturalist 184:119–131.https://doi.org/10.1086/676504
-
Simultaneous inference in general parametric modelsBiometrical Journal. Biometrische Zeitschrift 50:346–363.https://doi.org/10.1002/bimj.200810425
-
Measurement and meaning in biologyThe Quarterly Review of Biology 86:3–34.https://doi.org/10.1086/658408
-
Multi-state models for panel data: the msm package for RJournal of Statistical Software 38:1–28.https://doi.org/10.18637/jss.v038.i08
-
Natural selection mediated by seasonal time constraints increases the alignment between evolvability and developmental plasticityEvolution; International Journal of Organic Evolution 75:464–475.https://doi.org/10.1111/evo.14147
-
The mutation matrix and the evolution of evolvabilityEvolution; International Journal of Organic Evolution 61:727–745.https://doi.org/10.1111/j.1558-5646.2007.00071.x
-
Connecting QTLs to the g-matrix of evolutionary quantitative geneticsEvolution; International Journal of Organic Evolution 63:813–825.https://doi.org/10.1111/j.1558-5646.2008.00590.x
-
Antler size in red deer: heritability and selection but no evolutionEvolution; International Journal of Organic Evolution 56:1683–1695.https://doi.org/10.1111/j.0014-3820.2002.tb01480.x
-
Natural selection and random genetic drift in phenotypic evolutionEvolution; International Journal of Organic Evolution 30:314–334.https://doi.org/10.1111/j.1558-5646.1976.tb00911.x
-
The measurement of selection on correlated charactersEvolution; International Journal of Organic Evolution 37:1210–1226.https://doi.org/10.1111/j.1558-5646.1983.tb00236.x
-
Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilationJournal of Evolutionary Biology 22:1435–1446.https://doi.org/10.1111/j.1420-9101.2009.01754.x
-
Softwareemmeans: estimated marginal means, Aka least-squares meansR Package Version 1.7.1-1.
-
Morphological novelty emerges from pre-existing phenotypic plasticityNature Ecology & Evolution 2:1289–1297.https://doi.org/10.1038/s41559-018-0601-8
-
Elucidating the molecular architecture of adaptation via evolve and resequence experimentsNature Reviews. Genetics 16:567–582.https://doi.org/10.1038/nrg3937
-
Phenotypic evolution by neutral mutationEvolution; International Journal of Organic Evolution 40:915–935.https://doi.org/10.1111/j.1558-5646.1986.tb00561.x
-
SoftwareMallard_Robertson, version swh:1:rev:dbb9f78b983cf20dd1acdf0ae3a92e02d24a6782Software Heritage.
-
Explaining stasis: microevolutionary studies in natural populationsGenetica 112–113:199–222.
-
The role of developmental plasticity in evolutionary innovationProceedings. Biological Sciences 278:2705–2713.https://doi.org/10.1098/rspb.2011.0971
-
The danger of applying the breeder’s equation in observational studies of natural populationsJournal of Evolutionary Biology 23:2277–2288.https://doi.org/10.1111/j.1420-9101.2010.02084.x
-
The prediction of adaptive evolution: empirical application of the secondary theorem of selection and comparison to the breeder’s equationEvolution; International Journal of Organic Evolution 66:2399–2410.https://doi.org/10.1111/j.1558-5646.2012.01632.x
-
Increasing our ability to predict contemporary evolutionNature Communications 11:5592.https://doi.org/10.1038/s41467-020-19437-x
-
Phenotypic plasticity’s impacts on diversification and speciationTrends in Ecology & Evolution 25:459–467.https://doi.org/10.1016/j.tree.2010.05.006
-
Visualizing multivariate selectionEvolution; International Journal of Organic Evolution 43:1209–1222.https://doi.org/10.1111/j.1558-5646.1989.tb02569.x
-
Evolutionary Genetics: Concepts and Case StudiesEvolution of genetic variance-covariance structure, Evolutionary Genetics: Concepts and Case Studies, Oxford University Press.
-
The role of phenotypic plasticity in driving genetic evolutionProceedings. Biological Sciences 270:1433–1440.https://doi.org/10.1098/rspb.2003.2372
-
The missing response to selection in the wildTrends in Ecology & Evolution 33:337–346.https://doi.org/10.1016/j.tree.2018.02.007
-
The measurement of selection on quantitative traits: biases due to environmental covariances between traits and fitnessEvolution; International Journal of Organic Evolution 46:616–626.https://doi.org/10.1111/j.1558-5646.1992.tb02070.x
-
SoftwareR: A language and environment for statistical computingR Foundation for Statistical Computing, Vienna, Austria.
-
Adaptive radiation along genetic lines of least resistanceEvolution; International Journal of Organic Evolution 50:1766–1774.https://doi.org/10.1111/j.1558-5646.1996.tb03563.x
-
Thinking about the evolution of complex traits in the era of Genome-Wide Association studiesAnnual Review of Genomics and Human Genetics 20:461–493.https://doi.org/10.1146/annurev-genom-083115-022316
-
SoftwareRstan: the R interface to StanR Package Version 2.18.2.
-
Comparative quantitative genetics: evolution of the G matrixTrends in Ecology & Evolution 17:320–327.https://doi.org/10.1016/S0169-5347(02)02505-3
-
Estimating uncertainty in multivariate responses to selectionEvolution; International Journal of Organic Evolution 68:1188–1196.https://doi.org/10.1111/evo.12321
-
Correlational selection in the age of genomicsNature Ecology & Evolution 5:562–573.https://doi.org/10.1038/s41559-021-01413-3
-
High-throughput behavioral analysis in C. elegansNature Methods 8:592–598.https://doi.org/10.1038/nmeth.1625
-
Artificial selection to increase the phenotypic variance in gmax FailsThe American Naturalist 190:707–723.https://doi.org/10.1086/693959
-
Genetic variation for outcrossing among Caenorhabditis elegans isolatesEvolution; International Journal of Organic Evolution 60:1300–1305.
-
Phenotypic Plasticity of InsectsPhenotypic plasticity and Evolvability: an empirical test with experimental evolution, Phenotypic Plasticity of Insects, Science Publishers, Inc, 10.1201/b10201-22.
-
Phenotypic evolution, constant covariances, and the maintenance of additive varianceEvolution; International Journal of Organic Evolution 42:1342–1347.https://doi.org/10.1111/j.1558-5646.1988.tb04193.x
-
The C. elegans hypertonic stress response: big insights from shrinking wormsCellular Physiology and Biochemistry 55:89–105.https://doi.org/10.33594/000000332
-
Evolution of genetic variance during adaptive radiationThe American Naturalist 191:E108–E128.https://doi.org/10.1086/696123
-
Persistence of changes in the genetic covariance matrix after a bottleneckEvolution; International Journal of Organic Evolution 56:1968–1975.https://doi.org/10.1111/j.0014-3820.2002.tb00122.x
-
Towards evolutionary predictions: current promises and challengesEvolutionary Applications 16:3–21.https://doi.org/10.1111/eva.13513
-
C. elegans locomotion: small circuits, complex functionsCurrent Opinion in Neurobiology 33:117–126.https://doi.org/10.1016/j.conb.2015.03.009
Article and author information
Author details
Funding
European Research Council (ERC-St-243285)
- Henrique Teotónio
Agence Nationale pour la Recherche (ANR-14-ACHN-0032-01)
- Henrique Teotónio
Agence Nationale pour la Recherche (ANR-17-CE02-0017-01)
- Henrique Teotónio
National Science Foundation (PHY-1748958)
- Henrique Teotónio
Gordon and Betty Moore Foundation (2919.02)
- Henrique Teotónio
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank I Chelo, H Gendrot, T Guzella, and L Noble for help with worm handling, data acquisition or data analysis; C Baer, C Dillmann, T Long, L Noble, J Pemberton, P Phillips, S Proulx, J Sztepanacz and P de Villemereuil for discussion; J Sztepanacz, G Walter, and reviewers for suggestions on improving the presentation of this study. Funding: This work was supported by the European Research Council (ERC-St-243285), the Agence Nationale pour la Recherche (ANR-14-ACHN-0032–01, ANR-17-CE02-0017-01), and the KITP Quantitative Biology program (National Science Foundation, PHY-1748958; Gordon and Betty Moore Foundation, 2919.02).
Copyright
© 2023, Mallard et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,268
- views
-
- 249
- downloads
-
- 5
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Evolutionary Biology
Signs of ageing become apparent only late in life, after organismal development is finalized. Ageing, most notably, decreases an individual’s fitness. As such, it is most commonly perceived as a non-adaptive force of evolution and considered a by-product of natural selection. Building upon the evolutionarily conserved age-related Smurf phenotype, we propose a simple mathematical life-history trait model in which an organism is characterized by two core abilities: reproduction and homeostasis. Through the simulation of this model, we observe (1) the convergence of fertility’s end with the onset of senescence, (2) the relative success of ageing populations, as compared to non-ageing populations, and (3) the enhanced evolvability (i.e. the generation of genetic variability) of ageing populations. In addition, we formally demonstrate the mathematical convergence observed in (1). We thus theorize that mechanisms that link the timing of fertility and ageing have been selected and fixed over evolutionary history, which, in turn, explains why ageing populations are more evolvable and therefore more successful. Broadly speaking, our work suggests that ageing is an adaptive force of evolution.
-
- Evolutionary Biology
- Genetics and Genomics
The evolutionary origins of Bilateria remain enigmatic. One of the more enduring proposals highlights similarities between a cnidarian-like planula larva and simple acoel-like flatworms. This idea is based in part on the view of the Xenacoelomorpha as an outgroup to all other bilaterians which are themselves designated the Nephrozoa (protostomes and deuterostomes). Genome data can provide important comparative data and help understand the evolution and biology of enigmatic species better. Here, we assemble and analyze the genome of the simple, marine xenacoelomorph Xenoturbella bocki, a key species for our understanding of early bilaterian evolution. Our highly contiguous genome assembly of X. bocki has a size of ~111 Mbp in 18 chromosome-like scaffolds, with repeat content and intron, exon, and intergenic space comparable to other bilaterian invertebrates. We find X. bocki to have a similar number of genes to other bilaterians and to have retained ancestral metazoan synteny. Key bilaterian signaling pathways are also largely complete and most bilaterian miRNAs are present. Overall, we conclude that X. bocki has a complex genome typical of bilaterians, which does not reflect the apparent simplicity of its body plan that has been so important to proposals that the Xenacoelomorpha are the simple sister group of the rest of the Bilateria.