How development affects evolution

Abstract

Natural selection acts on developmentally constructed phenotypes, but how does development affect evolution? This question prompts a simultaneous consideration of development and evolution. However, there has been a lack of general mathematical frameworks mechanistically integrating the two, which may have inhibited progress on the question. Here, we use a new mathematical framework that mechanistically integrates development into evolution to analyse how development affects evolution. We show that, while selection pushes genotypic and phenotypic evolution up the fitness landscape, development determines the admissible evolutionary pathway, such that evolutionary outcomes occur at path peaks rather than landscape peaks. Changes in development can generate path peaks, triggering genotypic or phenotypic diversification, even on constant, single-peak landscapes. Phenotypic plasticity, niche construction, extra-genetic inheritance, and developmental bias alter the evolutionary path and hence the outcome. Thus, extra-genetic inheritance can have permanent evolutionary effects by changing the developmental constraints, even if extra-genetically acquired elements are not transmitted to future generations. Selective development, whereby phenotype construction points in the adaptive direction, may induce adaptive or maladaptive evolution depending on the developmental constraints. Moreover, developmental propagation of phenotypic effects over age enables the evolution of negative senescence. Overall, we find that development plays a major evolutionary role.

Evo-devo, genetic constraints, G matrix, aging, evolutionary dynamics, fitness landscapes

Development may be defined as the process that constructs the phenotype across the lifespan. Natural selection must then act upon what development constructs, raising the question of how development affects evolution. Addressing this question might benefit from an integrated consideration of development and evolution, attending to both the mechanisms of evolution and the mechanisms of development. Efforts to address this question have increasingly invoked a need for a synthesis analogous to the modern synthesis between natural selection and genetics of the early twentieth century (Alberch et al., 1979; Goldschmidt, 1940; Gould, 1980; Gould & Lewontin, 1979; Laland et al., 2014, 2015; Pigliucci, 2007; Pigliucci & Müller, 2010; Waddington, 1957; West-Eberhard, 2003). Key early steps in the modern synthesis involved the formulation of general mathematical frameworks that integrated the mechanisms of natural selection and inheritance, namely, population and quantitative genetics. Similarly, general mathematical frameworks integrating the mechanisms of development into evolution might help advance our understanding of the evolutionary effects of development. Large bodies of empirical and theoretical research, as well as entire scientific journals and scientific societies have been increasingly devoted to establish the effects of development on evolution from different angles (Carroll, 2008; Müller, 2007; Schlosser & Wagner, 2004; Wagner, 2005). However, despite steady growth in research efforts on this topic, there remains a lack of general mathematical frameworks synthesizing the mechanisms of development and evolution.

One of the available frameworks, namely quantitative genetics, provides some of the still fundamental understanding of how development affects evolution. Stemming from Fisher (1918), quantitative genetics describes development by relating genotype to phenotype in terms of regression coefficients (i.e., Fisher’s additive effects of allelic substitution). From this basis, quantitative genetics established a basic principle of adaptation for phenotypic evolution as the climbing of an adaptive topography (Lande, 1979). According to this principle, under certain assumptions, evolution by natural selection can be seen as the climbing of a fitness landscape over phenotype space, where selection pushes evolution in the direction of steepest ascent in the landscape and genetic covariation may divert evolution in a less steep direction. This principle indicates that development affects evolution by shaping the genetic covariation upon which selection acts, where this covariation is quantified by the $\mathbf{G}$ matrix (Charlesworth et al., 1982; Cheverud, 1984; Klingenberg, 2010; Lande, 1979; Maynard Smith et al., 1985). This perspective gives development two fundamental possibilities. First, if development yields genetic variation in all directions of phenotype space (such that $\mathbf{G}$ is non-singular), then development may affect the direction of evolution but evolutionary outcomes still occur at peaks of the fitness landscape. Second, if development yields no genetic variation in some directions of phenotype space (such that $\mathbf{G}$ is singular), then evolution can stop away from fitness landscape peaks in which case development not only affects the direction of evolution but may also affect where the evolutionary outcomes are in the fitness landscape (Altenberg, 1995; Houle, 1991; Kirkpatrick & Lofsvold, 1992). The first possibility gives development a relatively minor role in evolution, as in such case selection alone pre-defines the evolutionary outcomes (i.e., landscape peaks) and development can at most only influence which outcome is achieved. Instead, the second possibility might give development a relatively major evolutionary role, as selection and development would jointly define the evolutionary outcomes (i.e., away from landscape peaks).

Progress on establishing which of the two possibilities is the case has been difficult. Despite substantial research effort, the answer has remained uncertain, in part due to the difficulties of empirically establishing whether $\mathbf{G}$ is singular (Kirkpatrick, 2009). More fundamentally, one reason for this uncertainty may have been that equations describing phenotypic evolution have been derived from a regression description of development, where trait values are regressed on gene content, yielding a regression-based understanding of genetic covariation (Fisher, 1918; Lande, 1979). Since regression describes relationships regardless of their underlying mechanisms, this approach to genetic covariation may have limited a mechanistic understanding of the nature of genetic covariation arising from development. In principle, a mechanistic integration of development into evolution could yield a mechanistic understanding of genetic covariation and of the $\mathbf{G}$-matrix.

The task of devising a general mathematical framework that mechanistically integrates development into evolution would involve many complexities, which has likely contributed to the persistent lack of such a framework. Calls for synthesis between developmental and evolutionary biology have asked for consideration of the mechanistic basis of phenotype construction, non-linear genotype-phenotype maps, gene-gene and gene-environment interactions, non-normal distributions, far-from-equilibrium evolutionary dynamics, dynamic fitness landscapes, evolution and the nature of the $\mathbf{G}$-matrix, evolvability and epigenetics, and a variety of other complexities (Pigliucci, 2007; Pigliucci & Schlichting, 1997). Calls for integration of development into evolution have also highlighted several developmental factors—namely phenotypic plasticity (West-Eberhard, 2003), niche construction (Odling-Smee et al., 1996), extra-genetic inheritance (Jablonka & Lamb, 2014), and developmental bias (Arthur, 2004)—as possibly having important, yet unrecognised, evolutionary consequences (Laland et al., 2014, 2015). Mathematical frameworks and specific mathematical models integrating some of these complexities have become available in recent decades (Caswell, 2001; Chantepie & Chevin, 2020; Day & Bonduriansky, 2011; Dieckmann & Law, 1996; Engen & Sæ ther, 2021; Hansen & Wagner, 2001; Lande, 2019; Mullon & Lehmann, 2018). Additionally, theoretical research has often used individual-based simulations integrating some these complexities (Jones et al., 2014a; Miloco & Salazar-Ciudad, 2022; Salazar-Ciudad & Marín-Riera, 2013; Watson et al., 2013). However, it is of particular interest to obtain general mathematical frameworks mechanistically integrating development into evolution to seek deeper insight—general in the sense of encompassing a broad class of models. Despite the progress made, a synthetic mathematical framework unifying these complexities in a general and tractable manner had remained unavailable until recently.

A new mathematical framework (González-Forero, 2021) integrates mechanistic development into evolution while incorporating the elements listed in the previous paragraph by building upon many advances in evolutionary modelling over the last decades (Caswell, 2001; Dieckmann & Law, 1996) (and references in González-Forero, 2021). This framework integrates conceptual and mathematical advances from adaptive dynamics (Dieckmann & Law, 1996), matrix population models (Caswell, 2001), and optimal control theory as used in life-history theory (Schaffer, 1983; Sydsæ ter et al., 2008). The framework yields formulas that relate mechanistic descriptions of development to genetic covariation and to plastic change separate from selection. Here we use this framework to analyze how development affects evolution. This analysis sharpens the principle of adaptation as the climbing of a fitness landscape and provides insights into a wide array of long-standing questions. We focus on developing conceptual understanding and refer the reader to González-Forero (2021) for technical details.

Methods

An overview of the framework is the following. The framework is based on standard assumptions of adaptive dynamics (Dieckmann & Law, 1996; Metz et al., 1996). It considers a resident (i.e., wild type) population where individuals can be of different ages (i.e., it is age structured), reproduction is clonal for simplicity, whereby offspring receive the same genotype of their parent, and individuals can interact socially but only with non-relatives for simplicity. The genetic architecture (e.g., number of loci, ploidy, or linkage) need not be specified given our assumption of clonal reproduction, but it may be specified in particular models. Each individual has three types of age-specific traits that we let take continuous values to take derivatives. First, individuals have genotypic traits, which we refer to as the genotype and which are directly specified by genes. For instance, genotypic traits may be a continuous representation of nucleotide sequence or may be life-history traits assumed to be under direct genetic control. Genotypic traits correspond to control variables in the terminology of optimal control theory. Second, individuals have phenotypic traits, which we refer to as the phenotype and which are constructed over life depending on the genotype, developmental history, environment, social interactions, and their interaction, and where such phenotype construction is subject to developmental constraints. For instance, phenotypic traits include morphology and behavior. Phenotypic traits correspond to state variables in the terminology of optimal control theory. Third, individuals have environmental traits, which we refer to as the environment and which describe the individual’s local environment, possibly modified by the individual, social partners, sources exogenous to the population, and their interaction, and where such environmental alteration is subject to environmental constraints. For instance, environmental traits include the individual’s ambient temperature or humidity, which the individual may adjust, such as by roosting in the shade. This terminology contrasts with standard adaptive dynamics terminology, which would call genotypic traits phenotypes while here the phenotype is only traits that are developed; also standard adaptive dynamics terminology would define the environment as global, including everything outside the individual while here the environment is only local to the individual to model niche construction. Once the resident population achieves carrying capacity, rare mutant individuals arise who have a marginally different genotype from the resident genotype, drawn from an unbiased distribution of the deviation of mutant genotypes from the resident. Thus, we assume mutation is rare, weak, and unbiased. Population dynamics is deterministic so the only source of stochasticity is mutation. We assume that if rare mutants increase in frequency, they achieve fixation, which establishes a new resident, thus yielding evolutionary change.

The framework uses the following notation. The genotype across all genotypic traits for a mutant individual of age $a\in \{1,\dots ,N_{\mathrm{a}}\}$ is described by the vector ${\mathbf{y}}_a=(y_{1a},\dots ,y_{N_{\mathrm{g}}a}{)}^{⊺}$ of $N_{\mathrm{g}}$ genotypic traits (control variables; e.g., y_ia may be a continuous variable representing—via a differentiable approximation of the Heaviside step function—whether nucleotide $I=\mathrm{c}\mathrm{e}\mathrm{i}\mathrm{l}(i/n)$ is present at locus $J=i-\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{r}((i-1)/n)n$ at age a for n loci). While the genotype of an individual may often be constant with age, we let it depend on age as genotypic traits vary with age in life-history models. The phenotype of a mutant individual of age a is described by the vector ${\mathbf{x}}_a=(x_{1a},\dots ,x_{N_{\mathrm{p}}a}{)}^{⊺}$ of $N_{\mathrm{p}}$ phenotypic traits (state variables; e.g., x_ia may be the size of trait i at age a). The environment of a mutant individual of age a is described by the vector ${\mathit{\epsilon }}_a=({\epsilon }_{1a},\dots ,{\epsilon }_{N_{\mathrm{e}}a}{)}^{⊺}$ of $N_{\mathrm{e}}$ mutually independent environmental traits (e.g., ${\epsilon }_{ia}$ may be the temperature experienced at age a). A mutant individual’s phenotype, genotype, and environment across all ages are described respectively by the block-column vectors $\mathbf{x}=({\mathbf{x}}_1;\dots ;{\mathbf{x}}_{N_{\mathrm{a}}})$⁠, $\mathbf{y}=({\mathbf{y}}_1;\dots ;{\mathbf{y}}_{N_{\mathrm{a}}})$⁠, and $\mathit{\epsilon }=({\mathit{\epsilon }}_1;\dots ;{\mathit{\epsilon }}_{N_{\mathrm{a}}})$ (the semicolon indicates a linebreak). The framework finds that two aggregates of traits have special importance. First, we aggregate genotype and phenotype into a single vector called the “geno-phenotype” (Feldman & Zhivotovsky, 1992 use the term phenogenotype). Thus, the geno-phenotype of a mutant individual of age a is described by ${\mathbf{z}}_a=({\mathbf{x}}_a;{\mathbf{y}}_a)$⁠. The mutant geno-phenotype across all ages is $\mathbf{z}=(\mathbf{x};\mathbf{y})$⁠. Second, we similarly aggregate genotype, phenotype, and environment into a single vector called the “geno-envo-phenotype”. The geno-envo-phenotype of a mutant individual of age a is described by the block-column vector ${\mathbf{m}}_a=({\mathbf{z}}_a;{\mathit{\epsilon }}_a)$⁠. The mutant geno-envo-phenotype across all ages is $\mathbf{m}=(\mathbf{z};\mathit{\epsilon })$⁠. Resident variables are denoted analogously with an overbar (e.g., $\overline{\mathbf{z}}$⁠). We assume that the mutant genotype $\mathbf{y}$ is a random variable, whose deviation from the resident is given by the mutational distribution $M(\mathbf{y}-\overline{\mathbf{y}})$ with mean zero and vanishingly small and unbiased variance (i.e., $0<\mathrm{t}\mathrm{r}(\mathrm{c}\mathrm{o}\mathrm{v}[\mathbf{y},\mathbf{y}])\ll 1$ and M is even).

To mechanistically incorporate development, we describe an individual’s phenotype at a given age as a function of her genotype, phenotype, and environment at the immediately preceding age and of the social interactions experienced at that age. Thus, the phenotype of a mutant individual at age a + 1 satisfies the developmental constraint

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{x}}_{a+1}={\mathbf{g}}_a({\mathbf{m}}_a,{\displaystyle \overline{\mathbf{z}}}),}\end{array}}$$

(1)

where the developmental map ${\mathbf{g}}_a({\mathbf{m}}_a,\overline{\mathbf{z}})=(g_{1a}({\mathbf{m}}_a,\overline{\mathbf{z}}),\dots ,g_{N_{\mathrm{p}}a}({\mathbf{m}}_a,\overline{\mathbf{z}}{))}^{⊺}$ depends on the mutant’s genotype, phenotype, and environment at age a and on the genotype and phenotype of social partners of any age. For simplicity, we assume that the developmentally initial phenotype ${\mathbf{x}}_1$ is constant. We also assume that the genotype $\mathbf{y}$ is developmentally independent, that is, each y_ia is entirely determined by the relevant genotypic trait (e.g., nucleotide presence at the relevant locus) but is not influenced by other genotypic traits (e.g., nucleotide presence at other loci), the phenotype, or environment (i.e., in the terminology of optimal control theory, we assume open-loop control; Sydsæ ter et al., 2008). Equation (1) is a constraint in that the phenotype $\mathbf{x}$ cannot take arbitrary values but only those satisfying that equation. The developmental map can evolve and take any (differentiable) form as the genotype, phenotype, and environment evolve. The developmental map can also change over development (i.e., the functions ${\mathbf{g}}_a$ may be different to ${\mathbf{g}}_j$ for $j\ne a$⁠), for example, with metamorphosis.

To mechanistically incorporate niche construction and plasticity, we describe an individual’s environment as a function of the genotype and phenotype of the individual or social partners, and of processes exogenous to the population. The dependence of the environment on the genotype and phenotype of herself and social partners can describe niche construction by her or her social partners. The dependence of the environment on exogenous processes can describe, for instance, eutrophication or climate change caused by members of other species. For simplicity, we assume that exogenous environmental change is slow to allow the population to achieve the carrying capacity. Thus, a mutant’s environment at age a satisfies the environmental constraint

$${\displaystyle \begin{array}{l}{\displaystyle {\mathit{\epsilon }}_a={\mathbf{h}}_a({\mathbf{z}}_a,{\displaystyle \overline{\mathbf{z}}},\tau ),}\end{array}}$$

(2)

where the environmental map ${\mathbf{h}}_a({\mathbf{z}}_a,\overline{\mathbf{z}},\tau )=(h_{1a}({\mathbf{z}}_a,\overline{\mathbf{z}},\tau ),\dots ,h_{N_{\mathrm{e}}a}({\mathbf{z}}_a,\overline{\mathbf{z}},\tau {))}^{⊺}$ depends on the mutant’s genotype and phenotype at that age (e.g., from behavioural choice or environmental modification), on her social partners’ genotypes and phenotypes at any age (e.g., from social niche construction), and on the evolutionary time $\tau$ due to exogenous environmental change. The environment can then change over development (i.e., ${\mathit{\epsilon }}_a$ may be different from ${\mathit{\epsilon }}_j$ for $a\ne j$⁠) and evolution (either as the population evolves or due to exogenous causes). Equation (2) is also a constraint in that the environment $\mathit{\epsilon }$ can only take values allowed by that equation. Although social interactions are part of the environment (⁠${\mathit{\epsilon }}_a$ depends on $\overline{\mathbf{z}}$⁠), the developmental map (1) may also directly depend on social interactions to allow for modelling social interactions without having to also consider environmental traits. A mutant’s fertility $f_a({\mathbf{m}}_a,\overline{\mathbf{m}})$ and survival probability $p_a({\mathbf{m}}_a,\overline{\mathbf{m}})$ at age a depend on her genotype, phenotype, and environment at that age, and on the genotype, phenotype, and environment of her social partners of any age.

The developmental constraint (1) incorporates: developmental bias, in that the phenotype $\mathbf{x}$ may be predisposed to develop in certain ways; phenotypic plasticity, in that the same genotype $\mathbf{y}$ can generate a different phenotype $\mathbf{x}$ under a different environment $\mathit{\epsilon }$⁠; adaptive phenotypic plasticity, for example, via somatic selection or reinforcement learning (e.g., if g_ia is proportional to the gradient $\mathrm{\partial }{g}_{ja}/\mathrm{\partial }{x}_{ia}$ of a material payoff g_ja earned at age a that increases fitness); niche construction, in that environmental traits depend on developmentally constructed traits of the individual and social partners; and extra-genetic inheritance, in that an individual’s phenotype at age a + 1 can be an identical or modified copy of non-relatives’ phenotypes—this includes Jablonka & Lamb’s (2014) notion that extra-genetic inheritance involves the copying or reconstruction of others’ phenotype rather than their genotype. Thus, the developmental constraint (1) formalizes Maynard Smith et al.’s (1985) verbal definition of developmental constraints, according to which such constraints are “biases on the production of variant phenotypes or limitations on phenotypic variability caused by the structure, character, composition, or dynamics of the developmental system”. Moreover, the developmental constraint (1) can evolve as explained above, be non-linear, and mechanistically describe development, gene-gene interaction, and gene-environment interaction. To explore suggestions of adaptation without selection (Laland et al., 2015), we allow for the possibility that the developmental map ${\mathbf{g}}_a$ depends on the selection gradient $\mathrm{\partial }w/\mathrm{\partial }{\mathbf{z}}_a$ (defined below), in which case we say that development is selective whereby development has information about the adaptive direction. Allowing for selective development does not violate the assumptions of the developmental constraint (1), even though such perfectly idealised selective development may be biologically infeasible. Mechanistically, selective development might occur to some extent with individual or social learning, or with somatic selection (West-Eberhard, 2003).

The framework describes the evolutionary developmental (evo-devo) dynamics as follows. At each evolutionary time $\tau$⁠, the developmental dynamics of the resident phenotype are given by the developmental constraint (1) evaluated at the resident genotype. In turn, the evolutionary dynamics of the resident genotype are given by the canonical equation of adaptive dynamics (Dieckmann & Law, 1996; Dieckmann et al., 2006; Durinx et al., 2008)

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\Delta }{\displaystyle \overline{\mathbf{y}}}}{\mathrm{\Delta }\tau }=\iota {\mathbf{H}}_{\mathbf{y}}\frac{\mathrm{d}w}{\mathrm{d}\mathbf{y}}.}\end{array}}$$

(3)

Throughout, derivatives are evaluated at $\mathbf{y}=\overline{\mathbf{y}}$ and are in matrix calculus notation (Supplementary Material section S1; we use an equal sign in Equation (3) and evolutionarily dynamic equations below but they are strictly first-order approximations). The evo-devo dynamics are thus given by the evolutionary dynamics of the genotype (Equation 3) and the concomitant developmental dynamics of the phenotype (Equation 1) (a juxtaposition already used by Parvinen et al., 2013 and Metz et al., 2016). The right-hand side of Equation (3) has three components. First, $\iota$ is a non-negative scalar measuring mutational input, which is proportional to the mutation rate and the carrying capacity. Second, ${\mathbf{H}}_{\mathbf{y}}=\mathrm{c}\mathrm{o}\mathrm{v}[\mathbf{y},\mathbf{y}]$ is the mutational covariance matrix (H for heredity). Third, $\mathrm{d}w/\mathrm{d}\mathbf{y}$ is the total selection gradient of the genotype, which measures total genotypic selection, that is, total directional selection on genes considering the ability of genes to affect the phenotype (from the chain rule, the total derivative of a function $f(x,y)$ with respect to y, where $x=g(y)$⁠, is $\mathrm{d}f/\mathrm{d}y=(\mathrm{\partial }f/\mathrm{\partial }x)(\mathrm{d}x/\mathrm{d}y)+\mathrm{\partial }f/\mathrm{\partial }y$⁠). In contrast, $\mathrm{\partial }w/\mathrm{\partial }\mathbf{y}$ measures direct directional selection on genes without considering such ability (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathbf{y}$ is traditionally assumed zero, but is not zero in standard life-history models where the genotype modulates resource allocation to fertility and is thus under direct selection; González-Forero, 2021). Thus, $\mathrm{\partial }w/\mathrm{\partial }\mathbf{y}$ is a vector that points in the direction of steepest increase in fitness in genotype space without constraints, whereas $\mathrm{d}w/\mathrm{d}\mathbf{y}$ is a vector that points in the direction of steepest increase in fitness in genotype space subject to the developmental (1) and environmental (2) constraints. González-Forero (2021) finds closed-form formulas for the total selection gradient of the genotype. Because of age-structure, a mutant’s relative fitness is

$${\displaystyle \begin{array}{l}{\displaystyle w={\displaystyle \sum _{a=1}^{N_{\mathrm{a}}}}{w}_a=\frac{1}{T}{\displaystyle \sum _{a=1}^{N_{\mathrm{a}}}}({\varphi }_a{f}_a+{\pi }_a{p}_a),}\end{array}}$$

(4)

where w_a is the mutant’s relative fitness at age a, T is generation time, and ${\varphi }_a$ and ${\pi }_a$ are the forces (Hamilton, 1966) of selection on fertility and survival at age a (given in Equations S5 in the Supplementary Material). The quantities T, ${\varphi }_a$⁠, and ${\pi }_a$ depend on resident but not mutant values. The total selection gradient of the genotype is a form of Caswell’s (1982, 2001) total derivative of fitness, Charlesworth’s (1994) total differential of the population’s growth rate, van Tienderen’s (1995) integrated sensitivity of the population’s growth rate, and Morrissey’s (2014, 2015) extended selection gradient.

Results

Recovery of classic but insufficient results

We investigate how development as described by Equation (1) affects evolution by analysing dynamically sufficient equations describing the long-term evolution of a developed phenotype $\overline{\mathbf{x}}$ as the climbing of a fitness landscape. Before arriving at such equations, we first describe a mechanistic version of the Lande (1979) equation and explain why it is generally insufficient to describe long-term evolution of developed phenotypes.

Assume for now that there is (i) no niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}=\mathbf{0}$⁠), (ii) no social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}=\mathbf{0}$⁠), (iii) no exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau =\mathbf{0}$⁠), and (iv) no direct genotypic selection (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathbf{y}=\mathbf{0}$⁠). Then, González-Forero (2021) shows that the (expected) evolutionary dynamics of the resident phenotype (as $\mathrm{\Delta }\tau \to 0$⁠) satisfy

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{x}}}}{\mathrm{d}\tau }=\iota {\mathbf{H}}_{\mathbf{x}}\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{x}},}\end{array}}$$

(5)

where ${\mathbf{H}}_{\mathbf{x}}$ guarantees that the developmental constraint (1) is satisfied at all times. The right-hand side of Equation (5) has two different components relative to those of Equation (3). First, $\mathrm{\partial }w/\mathrm{\partial }\mathbf{x}$ corresponds to Lande’s selection gradient, which measures direct directional selection on the phenotype: it is a vector that points in the direction of steepest increase in fitness in phenotype space without constraints. Second, Equation (5) depends on the mechanistic additive genetic covariance matrix of the phenotype

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{H}}_{\mathbf{x}}=\mathrm{c}\mathrm{o}\mathrm{v}[{\mathbf{b}}_{\mathbf{x}},{\mathbf{b}}_{\mathbf{x}}]=\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}{\mathbf{y}}^{⊺}}{\mathbf{H}}_{\mathbf{y}}\frac{\mathrm{d}{\mathbf{x}}^{⊺}}{\mathrm{d}\mathbf{y}},}\end{array}}$$

(6)

where the mechanistic breeding value of $\mathbf{x}$ is defined as

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{b}}_{\mathbf{x}}={\displaystyle \overline{\mathbf{x}}}+\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}{\mathbf{y}}^{⊺}}(\mathbf{y}-{\displaystyle \overline{\mathbf{y}}}).}\end{array}}$$

(7)

The key difference of Equation (7) with breeding value is that the latter is defined using a matrix of regression coefficients in place of $\mathrm{d}\mathbf{x}/\mathrm{d}{\mathbf{y}}^{⊺}$ (Lynch & Walsh, 1998; Walsh & Lynch, 2018). Such regression coefficients are Fisher’s (1918) additive effects of allelic substitution and the matrix formed by them is Wagner’s (1984) developmental matrix. Whereas such regression coefficients are obtained via least squares, the matrix $\mathrm{d}\mathbf{x}/\mathrm{d}{\mathbf{y}}^{⊺}$ is obtained from the mechanistic description of development in Equations (1) and (2). Thus, the matrix $\mathrm{d}\mathbf{x}/\mathrm{d}{\mathbf{y}}^{⊺}$ is a mechanistic counterpart of both Fisher’s (1918) additive effects of allelic substitution and Wagner’s (1984) developmental matrix. Consequently, ${\mathbf{H}}_{\mathbf{x}}$ is a mechanistic counterpart of Lande’s (1979),$\mathbf{G}$-matrix. González-Forero (2021) obtained closed-form formulas for $\mathrm{d}\mathbf{x}/\mathrm{d}{\mathbf{y}}^{⊺}$ from the developmental (1) and environmental (2) constraints, which allow for a mechanistic understanding of genetic covariation. The fact that mechanistic breeding value is obtained from $\mathrm{d}\mathbf{x}/\mathrm{d}{\mathbf{y}}^{⊺}$ rather than regression coefficients means that breeding value and mechanistic breeding value have different properties. In particular, while breeding value is uncorrelated with the residual of predicting the phenotype because of least squares, mechanistic breeding value can be correlated with the residual (i.e., for a single developed trait $x=b_x+e_x$ with residual e_x, but possibly $\mathrm{c}\mathrm{o}\mathrm{v}[x,e_x]\ne 0$⁠). Consequently, the classic partition of phenotypic variance into additive genetic and “environmental” variances does not hold with mechanistic breeding value. Hence, narrow-sense mechanistic heritability of a single developed trait x defined with the traditional formula $\mathrm{v}\mathrm{a}\mathrm{r}[b_x]/\mathrm{v}\mathrm{a}\mathrm{r}[x]$ can be greater than one.

Equation (5) may then be understood as a mechanistic Lande (1979) equation, apart from some differences due to our adaptive dynamics assumptions but these differences do not affect the points we make. Indeed, the expression for relative fitness (4) has the form of that in the Lande (1982) equation for quantitative genetics under age structure. The expression for the mechanistic additive genetic covariance matrix of the phenotype (6) has the same form of previous expressions obtained under quantitative genetics assumptions, which in place of our closed-form derivatives had regression coefficients or derivatives of unknown form (see Equation II of Fisher, 1918, Equation + of Wagner, 1984, Equation 3.5b of Barton & Turelli, 1987, and Equation 4.23b of Lynch & Walsh, 1998; see also Equation 22a of Lande, 1980, Equation 3 of Wagner, 1989, and Equation 9 of Charlesworth, 1990). While in quantitative genetics the matrix corresponding to ${\mathbf{H}}_{\mathbf{y}}$ would describe the realized standing variation in gene content, here ${\mathbf{H}}_{\mathbf{y}}$ describes the expected variation in genotype due only to mutation in the current evolutionary time and so we call it the mutational covariance matrix. This different meaning of ${\mathbf{H}}_{\mathbf{y}}$ in the two frameworks does not mean that ${\mathbf{H}}_{\mathbf{y}}$ in our framework induces absolute constraints that are absent under quantitative genetics assumptions: our definition of ${\mathbf{H}}_{\mathbf{y}}$ allows for mutational variation in all directions of genotype space, or only in some directions (i.e., ${\mathbf{H}}_{\mathbf{y}}$ may be full rank, or less than full rank), as in quantitative genetics. While Lande’s (1979),$\mathbf{G}$ matrix is the covariance matrix of breeding values of the phenotypes, ${\mathbf{H}}_{\mathbf{x}}$ is the covariance matrix of the mechanistic breeding values of the phenotypes. Mechanistic breeding value (7) is the first-order estimate of the phenotype with respect to genotype as predictor, which corresponds to the quantitative genetics notion of breeding value as the best linear estimate of the phenotype with respect to gene content as predictor (Falconer & Mackay, 1996; Lynch & Walsh, 1998). Hence, genotypic traits play here an analogous role to that of gene content in quantitative genetics in that they emerge as the relevant first-order predictor of the phenotype to describe inheritance.

In contrast to the Lande (1979) equation, the mechanistic Lande Equation (5) has been derived from a mechanistic account of development so we have formulas to relate development (1) to mechanistic breeding value and so to mechanistic genetic covariation ${\mathbf{H}}_{\mathbf{x}}$⁠. These formulas guarantee that the developmental (1) and environmental (2) constraints are satisfied at all times. The formulas show that ${\mathbf{H}}_{\mathbf{x}}={\mathbf{H}}_{\mathbf{x}}(\overline{\mathbf{x}},\overline{\mathbf{y}})$ is generally a function of the resident genotype $\overline{\mathbf{y}}$⁠, in particular, if the development map is non-linear, in which case the developmental matrix $\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\mathbf{y}$ depends on $\overline{\mathbf{y}}$⁠. Specifically, gene-gene interaction (products between y’s in the developmental map), gene-phenotype interaction (products between y’s and x’s), and gene-environment interaction (products between y’s and $\epsilon$’s) cause the resident genotype $\overline{\mathbf{y}}$ to remain in $\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\mathbf{y}$⁠. Thus, in general, the right-hand side of the mechanistic Lande Equation (5) depends on the resident genotype $\overline{\mathbf{y}}$ if the mutation rate, carrying capacity, mutational covariation ${\mathbf{H}}_{\mathbf{y}}$⁠, the developmental matrix $\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\mathbf{y}$⁠, or direct phenotypic selection $\mathrm{\partial }w/\mathrm{\partial }\mathbf{x}$ depend on the resident genotype $\overline{\mathbf{y}}$⁠. Consequently, Equation (5) is under-determined as it generally depends on the resident genotype but does not describe the evolutionary dynamics of the resident genotype. Hence, Equation (5) cannot be used alone to describe the evolutionary dynamics of the resident phenotype in general.

This problem equally applies to the Lande equation. Since breeding value depends on the regression coefficients of phenotype on gene content at the current time, such regression coefficients depend on the current state of the population, including allele frequency. Thus, the Lande equation depends on allele frequency but it does not describe allele frequency change. The standard approach to address this problem is to assume Fisher’s (1918) infinitesimal model, whereby each phenotype is controlled by an arbitrarily large number of loci so in the short term allele frequency change per locus is negligible (Barton et al., 2017; Hill, 2017; Walsh & Lynch, 2018). Then, $\mathbf{G}$ is assumed constant or its evolution is described by the Bulmer equation, which considers change in $\mathbf{G}$ due to change in linkage disequilibrium while still assuming negligible allele frequency change (Barton et al., 2017; Lande & Arnold, 1983; Turelli, 1988). Such an approach allows for the Lande equation to describe evolution only in the short term, where allele frequency change per locus remains negligible (Barton et al., 2017). However, long-term evolution involves non-negligible allele frequency change, thus limiting the ability of the Lande equation to describe long-term phenotypic evolution.

In turn, the canonical equation of adaptive dynamics (Dieckmann & Law, 1996) describes long-term phenotype evolution, but where the phenotype does not have developmental constraints. In our terminology, the canonical equation of adaptive dynamics describes the long-term evolution of genotypic traits, which do not have developmental constraints (Equation 3; Dieckmann & Law, 1996). Thus, the canonical equation does not describe the long-term evolution of developed traits as an adaptive topography. The canonical equation for function-valued traits (Dieckmann et al., 2006) describes the evolution of genotypic traits that can affect the construction of phenotypic traits (Parvinen et al., 2013), but this equation still describes evolution in gradient form only for genotypic traits, not phenotypic traits (i.e., not for state variables; indeed, given the age-block-structure of $\mathbf{y}$⁠, Equation (3) may be understood as the canonical equation for function-valued traits in discrete age). Equations describing the long-term evolution in gradient form of phenotypic traits (i.e., with explicit developmental constraints, that is, of state variables) such as Equation (5), where formulas for ${\mathbf{H}}_{\mathbf{x}}$ guarantee that the developmental constraints are met, do not seem to have been available until recently (González-Forero, 2021). Yet, the dynamic insufficiency of Equation (5) limits its ability to describe the evolution of developed traits in a fitness landscape. (For a different view of how the canonical and Lande equations are related see pp. 1084-1086 of Geritz et al., 2016).

Development blocks evolutionary change

Even though the mechanistic Lande Equation (5) is insufficient to describe the long-term evolution of developed phenotypes particularly if development is non-linear, such equation can be extended to allow for this.

To see this, remove the assumption (iv) above, and still assume that there is (i) no niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}=\mathbf{0}$⁠), (ii) no social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}=\mathbf{0}$⁠), and (iii) no exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau =\mathbf{0}$⁠). Then, González-Forero (2021) shows that the evolutionary dynamics of the resident phenotype and genotype are given by

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{z}}}}{\mathrm{d}\tau }=\iota {\mathbf{H}}_{\mathbf{z}}\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{z}},}\end{array}}$$

(8)

where ${\mathbf{H}}_{\mathbf{z}}$ guarantees that the developmental constraint (1) is satisfied at all times. Equation (8) now depends on the mechanistic additive genetic covariance matrix of the geno-phenotype

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{H}}_{\mathbf{z}}=\mathrm{c}\mathrm{o}\mathrm{v}[{\mathbf{b}}_{\mathbf{z}},{\mathbf{b}}_{\mathbf{z}}]=\frac{\mathrm{d}\mathbf{z}}{\mathrm{d}{\mathbf{y}}^{⊺}}{\mathbf{H}}_{\mathbf{y}}\frac{\mathrm{d}{\mathbf{z}}^{⊺}}{\mathrm{d}\mathbf{y}},}\end{array}}$$

(9)

where the mechanistic breeding value ${\mathbf{b}}_{\mathbf{z}}$ of the geno-phenotype is given by Equation (7) replacing $\mathbf{x}$ with $\mathbf{z}$⁠.

Equation (8) is an extended mechanistic Lande equation that, in contrast to the Lande equation or its mechanistic counterpart (5), is dynamically sufficient and so describes the long-term evolution of developed traits as an adaptive topography. Equation (8) is dynamically sufficient because it describes the dynamics of all the variables involved by simultaneously describing genotypic and phenotypic evolution since $\mathbf{z}=(\mathbf{x};\mathbf{y})$ includes both genotype and phenotype. In particular, the matrix ${\mathbf{H}}_{\mathbf{z}}={\mathbf{H}}_{\mathbf{z}}(\overline{\mathbf{x}},\overline{\mathbf{y}})$ is still a function of the resident phenotype and genotype, but now the evolutionary dynamics of both are described by Equation (8). The right-hand side of Equation (8) may also be a function of the resident environmental traits, but because of assumptions (i) and (iii), niche construction and exogenous environmental change are absent so the environment remains constant.

As Equation (8) describes the long-term evolution of developed traits as an adaptive topography, we can now use it to analyse the effect of development on genetic covariation and evolution. A biologically crucial property of Equation (8) is that it implies that there are always absolute genetic constraints to adaptation within the assumptions made (absolute constraints mean that the constraining matrix is singular; Houle, 2001; Klingenberg, 2005, 2010). Indeed, Equation (8) simultaneously describes genotypic and phenotypic evolution, so the climbing of the fitness landscape is in geno-phenotype space rather than only phenotype space as in the Lande equation. Yet, because the geno-phenotype $\mathbf{z}$ contains the genotype $\mathbf{y}$⁠, genetic covariation in geno-phenotype is absolutely constrained by the genotypic space (i.e., ${\mathbf{H}}_{\mathbf{z}}$ is always singular because $\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}\mathbf{y}$ has fewer rows than columns). That is, to achieve dynamic sufficiency, one cannot generally consider any traits for evolutionary analysis, but must consider both the phenotype and genotype, and since the phenotype and genotype are related by development, genetic covariation in geno-phenotype space is absolutely constrained to the region where the relationship between the two is met. Thus, there cannot be genetic variation in as many directions in geno-phenotype space as there are phenotypes across life (i.e., ${\mathbf{H}}_{\mathbf{z}}$ has at least $N_{\mathrm{a}}{N}_{\mathrm{p}}$ eigenvalues that are exactly zero). Along such directions evolutionary change is blocked (i.e., along the directions of the eigenvectors corresponding to the zero eigenvalues of ${\mathbf{H}}_{\mathbf{z}}$⁠). In this sense, development can be seen as blocking evolutionary change in some directions. Therefore, evolutionary stasis (⁠$\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau =\mathbf{0}$⁠) can occur away from landscape peaks, where direct directional selection persists (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}\ne \mathbf{0}$⁠; see also Houle, 1991; Kirkpatrick & Lofsvold, 1992, and Altenberg, 1995). Thus, consideration of direct directional selection alone (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}$⁠) is not sufficient for predicting possible evolutionary outcomes, which depend also on ${\mathbf{H}}_{\mathbf{z}}$⁠. Predicting evolutionary outcomes still depends on ${\mathbf{H}}_{\mathbf{z}}$ even if there is no genotypic selection (i.e., if $\mathrm{\partial }w/\mathrm{\partial }\mathbf{y}=\mathbf{0}$⁠).

The singularity of ${\mathbf{H}}_{\mathbf{z}}$ and generalizations thereof has important implications, which we detail below. One immediate implication is that, without absolute mutational constraints (i.e., if ${\mathbf{H}}_{\mathbf{y}}$ is non-singular), evolutionary outcomes are jointly defined by direct selection and development as described by $\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\mathbf{y}$⁠, which is influenced by gene-gene, gene-phenotype, and gene-environment interactions. Another immediate implication is that this singularity may help explain common empirical observations of evolutionary stasis in wild populations despite directional selection and genetic variation, termed the “paradox of stasis” (Kingsolver & Diamond, 2011; Kirkpatrick, 2009; Merilä et al., 2001). Because ${\mathbf{H}}_{\mathbf{z}}$ is singular, evolutionary stasis is expected to generally occur with persistent directional selection and genetic variation, although there are other explanations for the paradox of stasis, including measurement error (Estes & Arnold, 2007; Haller & Hendry, 2013; Merilä et al., 2001; Morrissey, 2015).

Even though plasticity may be present in the extended mechanistic Lande Equation (8), plasticity has no evolutionary effect under the assumptions made in that equation. Indeed, the formulas for ${\mathbf{H}}_{\mathbf{z}}$ show that, under the framework’s assumptions, plasticity (i.e., $\mathrm{\partial }{\mathbf{z}}^{⊺}/\mathrm{\partial }\mathit{\epsilon }$⁠) only affects the evolutionary dynamics by interacting with niche construction or exogenous environmental change (see below). Thus, although in Equation (8) there may be plasticity, it has no effect since by assumptions (i) and (iii) niche construction and exogenous environmental change are absent so the environment remains constant, but plasticity requires environmental change.

Development determines the path

While development can be seen as blocking evolutionary change, development can be more specifically seen as determining the admissible evolutionary path. This is because one must follow the evolution of the genotype and phenotype for a dynamically sufficient phenotypic adaptive topography, and the phenotype is related to the genotype by the developmental constraint (1) along which evolution is allowed. Consequently, if there are no absolute mutational constraints and no exogenous environmental change, development and selection jointly define the evolutionary equilibria and development determines which of these equilibria are admissible.

To see this, continue to assume that there is (i) no niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}=\mathbf{0}$⁠), (ii) no social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}=\mathbf{0}$⁠), and (iii) no exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau =\mathbf{0}$⁠). Then, the evolutionary dynamics of the resident phenotype and genotype are equivalently given by

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{z}}}}{\mathrm{d}\tau }=\iota {\mathbf{H}}_{\mathbf{z}\mathbf{y}}\frac{\mathrm{d}w}{\mathrm{d}\mathbf{y}},}\end{array}}$$

(10)

where ${\mathbf{H}}_{\mathbf{z}\mathbf{y}}$ guarantees that the developmental constraint (1) is satisfied at all times. This equation is closely related to Morrissey’s (2014) Equation 4 but they differ as Morrissey’s equation is for the evolutionary change of the phenotype with linear development rather than geno-phenotype with possibly non-linear development, considers the total selection gradient of the phenotype rather than the genotype, considers regression-based rather than mechanistic genetic covariation (if the classic partition of phenotypic variance is to hold), and is dynamically insufficient as it is a rearrangement of the Lande (1979) equation. Equation (10) depends on the total selection gradient of the genotype, which under assumptions (i–iii), is

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}w}{\mathrm{d}\mathbf{y}}=\frac{\mathrm{d}{\mathbf{z}}^{⊺}}{\mathrm{d}\mathbf{y}}\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{z}}=\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{y}}+\frac{\mathrm{d}{\mathbf{x}}^{⊺}}{\mathrm{d}\mathbf{y}}\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{x}}.}\end{array}}$$

(11)

Thus, total genotypic selection depends on the developmental matrix $\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\mathbf{y}$ (see Figure S1 for a description of all the factors affecting total genotypic selection removing assumptions i-iii). Equation (10) also depends on the mechanistic additive genetic cross-covariance matrix between the geno-phenotype and genotype

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{H}}_{\mathbf{z}\mathbf{y}}=\mathrm{c}\mathrm{o}\mathrm{v}[{\mathbf{b}}_{\mathbf{z}},{\mathbf{b}}_{\mathbf{y}}]=\frac{\mathrm{d}\mathbf{z}}{\mathrm{d}{\mathbf{y}}^{⊺}}{\mathbf{H}}_{\mathbf{y}}.}\end{array}}$$

(12)

Crucially, the transpose matrix of the total effects of genotype on geno-phenotype, $\mathrm{d}\mathbf{z}/\mathrm{d}{\mathbf{y}}^{⊺}$⁠, is non-singular because the genotype is developmentally independent by assumption. Thus, Equation (10) implies that evolutionary stasis occurs when total genotypic selection vanishes (⁠$\mathrm{d}w/\mathrm{d}\mathbf{y}=\mathbf{0}$⁠) if there are no absolute mutational constraints. That is, total genotypic selection can identify evolutionary equilibria, even though direct directional selection typically cannot. Because total genotypic selection depends on the developmental matrix (Equation 11), this is a different way of showing that evolutionary equilibria are generally defined jointly by direct selection and development if there are no absolute mutational constraints. Moreover, there are always an infinite number of evolutionary equilibria, because $\mathrm{d}w/\mathrm{d}\mathbf{y}=\mathbf{0}$ provides fewer equations than those in $\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau =\mathbf{0}$⁠. Since the developmental constraint (1) must be satisfied at all times, this constraint provides the remaining equations and the admissible evolutionary trajectory. That is, the developmental constraints not only influence the evolutionary equilibria, but also provide the admissible evolutionary trajectory and so the admissible equilibria.

To gain intuition, consider the following simple example. Let there be one phenotype x and one genotypic trait y, and let age structure be collapsed such that development occurs instantaneously. Collapsing age structure is a heuristic simplification, but doing so makes $\mathbf{x}$ and $\mathbf{y}$ scalars, which can then be easily visualised. Also, let there be no social interactions, no niche construction, no exogenous environmental change, and no density dependence. By removing social interactions and density dependence, the latter of which is also a heuristic simplification, evolutionary change can be not only abstractly but also visually described as the climbing of a fitness landscape. The extended mechanistic Lande Equation (8) states that selection can be seen as “pushing” the geno-phenotype uphill on the fitness landscape in geno-phenotype space in the direction of steepest ascent, whereas genetic covariation diverts evolutionary change in a possibly less steep direction. If genetic constraints on the traits considered were not absolute (here, if ${\mathbf{H}}_{\mathbf{z}}$ were non-singular), the population would eventually reach a landscape peak, which is commonly implicitly or explicitly assumed. Yet, as the traits considered must be both the genotype and the phenotype arising from the developmental constraint (1), genetic constraints are necessarily absolute (i.e., ${\mathbf{H}}_{\mathbf{z}}$ is singular). Hence, evolution is restricted to an admissible path on the landscape where the developmental constraint is met (Figure 1a–c; the computer code used to generate all the figures is in the Supplementary Data). Adaptive evolution may thus be understood as the climbing of the fitness landscape along an admissible path determined by development. The evolutionary process eventually reaches a path peak if there are no absolute mutational constraints (Figure 1a–c). Formally, a path peak is a point in geno-phenotype space where the developmental constraints are met and fitness is locally and totally maximized with respect to genotype: here “totally” means after substituting in the fitness function for both the developmental and environmental constraints. Selection response vanishes at path peaks which are not necessarily landscape peaks (Figure 1a–c). The admissible path yields an elevation profile of fitness, namely, the total fitness landscape of the genotype, on which adaptation is constrained to occur, and can have peaks and valleys created by the developmental map (Figure 1d–f). Hence, selection pushes the genotype uphill on the total fitness landscape of the genotype, but in this profiled landscape evolutionary change is not blocked provided that there are no absolute mutational constraints. Consequently, without absolute mutational constraints, evolutionary stasis generally occurs at a peak on the total fitness landscape of the genotype, even though this does not generally correspond to a peak on the fitness landscape of the phenotype and genotype. The ${\mathbf{H}}_{\mathbf{z}}$-matrix evolves as an emergent property as the resident phenotype and genotype evolve, and genetic variances and covariances may substantially increase or decrease if development is non-linear (consistent with previous individual-based simulations; Miloco & Salazar-Ciudad, 2022), rather than being approximately constant as is often considered under short-term evolution with an infinite number of loci (Barton et al., 2017; Fisher, 1918) which we do not assume (Figure 1g–i). Ultimately, development constrains the path of adaptation and thus defines its outcome jointly with selection (Figure 1j–l). Overall, development has a permanent dual role in adaptation by influencing evolutionary equilibria (blue lines in Figure 1j–l) and determining the admissible evolutionary path (red line in Figure 1j–l).

Figure 1.

Development determines the evolutionary path. a–c, Fitness landscape for a phenotype x and a genotype y, with three developmental maps. The developmental map determines the path on the landscape as x can only take values along the path. d–f, The total fitness landscape of the genotype gives the elevation profile considering developmental and environmental constraints. g–i, Mechanistic additive genetic covariances depend on the developmental map and evolve as the genotype evolves. j–l, The evolutionary dynamics occur along the developmental constraint (red; gray arrows are parallel or antiparallel to leading eigenvectors of ${\mathbf{H}}_{\mathbf{z}}$⁠, called “genetic lines of least resistance”; Schluter, 1996). Evolutionary outcomes (black dots) occur at path peaks and thus depend on development. Evolutionary equilibria (blue) are infinite in number because ${\mathbf{H}}_{\mathbf{z}}$ is singular and occur when total genotypic selection vanishes despite persistent direct selection. The intersection of the developmental map and evolutionary equilibria yields the admissible evolutionary equilibria. Different developmental maps yield different evolutionary trajectories and outcomes with the same single-peak landscape. Developmental bias is quantified by the slope of the red line (⁠$\mathrm{\partial }x/\mathrm{\partial }y$⁠). Evolutionary change satisfies $\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau ={\mathbf{H}}_{\mathbf{z}}\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}$⁠, where ${\mathbf{H}}_{\mathbf{z}}=(\mathrm{d}\mathbf{z}/\mathrm{d}y)H_y(\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}y)=(\mathrm{d}x/\mathrm{d}y{,1)}^{⊺}{H}_y(\mathrm{d}x/\mathrm{d}y,1)=H_y((\mathrm{d}x/\mathrm{d}y{)}^2,\mathrm{d}x/\mathrm{d}y;\mathrm{d}x/\mathrm{d}y,1)=(H_x,H_{xy};H_{xy},H_y)$ which is singular. Hence, the mechanistic additive genetic variance of x is $H_y{(\mathrm{d}x/\mathrm{d}y)}^2$ . The total effects of a mutant’s genotype on her geno-phenotype are $\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}y=(\mathrm{d}x/\mathrm{d}y,\mathrm{d}y/\mathrm{d}y)=(\mathrm{d}x/\mathrm{d}y,1)$ . Fitness is $w(x,y)=\exp (-x^2-y^2)$ and mutational variance is H_y = 0.01 (so that $0<\mathrm{t}\mathrm{r}({\mathbf{H}}_{\mathbf{y}})\ll 1$⁠). As we use our heuristic assumptions of collapsed age structure and no density dependence, we use here a simpler fitness function than that in Equation (4). Functional forms are chosen for illustration.

That evolution stops at path peaks rather than landscape peaks may help explain abundant empirical data. Empirical estimation of selection has found that direct stabilizing selection is rare relative to direct directional selection across thousands of estimates in more than 50 species (Kingsolver & Diamond, 2011; Kingsolver et al., 2001). The rarity of stabilizing selection has been puzzling under the common assumption that evolutionary outcomes occur at landscape peaks, where stabilizing selection is prevalent (Charlesworth et al., 1982; Walsh & Lynch, 2018). In contrast, the rarity of stabilizing selection is consistent with evolutionary outcomes necessarily occurring at path peaks. Indeed, if path peaks occur outside landscape peaks, evolutionary outcomes occur with persistent direct directional selection (Altenberg, 1995; Houle, 1991; Kirkpatrick & Lofsvold, 1992) and relatively weak direct stabilizing selection (because direct stabilizing selection is relatively strong only near landscape peaks). Thus, evolution generally stopping at path peaks outside landscape peaks may help explain the otherwise puzzling observation of rare stabilizing selection. There are other explanations for the rarity of stabilizing selection (Haller & Hendry, 2013; Morrissey, 2015), and Morrissey’s is closely related to ours.

Development can drive diversification

The evolution of phenotypic and genetic diversity is typically explained in terms of selection or drift (Doebeli, 2011; Gavrilets & Losos, 2009), but the relevance of development has been less clear. We find that development can drive the evolution of phenotypic and genotypic diversity, even in a constant, single-peak fitness landscape, where selection provides an evolutionary force while development translates it into diversification. This can happen in two ways. First, if the developmental map changes inducing a shifted path peak, the population can evolve different phenotypes and genotypes (Figure 2a–f). Second, the developmental map may change, generating new path peaks, if it is or becomes non-linear, such that the population may evolve different phenotypes and genotypes depending on the initial conditions (Figure 2g–l). Diversification may then occur if a population subdivides and one of the descendant populations undergoes changes in their developmental map due to changes in its genotype, environment, or how the phenotype is constructed. Mathematically, multiple path peaks can arise from non-linear development because it can generate bifurcations altering the number of admissible evolutionary equilibria even though the number of evolutionary equilibria remains the same (i.e., infinite). This is a form of evolutionary branching (Geritz et al., 1998), which may be called evo-devo branching, where the admissible path creates a valley in the total fitness landscape of the genotype so there is disruptive total genotypic selection (Figure 1e; see also Morrissey, 2015 where stabilizing total selection is analysed). There may not be a valley in the fitness landscape so no disruptive direct selection on the phenotype or genotype (Figure 1b) as evo-devo branching only requires disruptive total genotypic selection possibly created by development rather than direct disruptive selection. Thus, evo-devo branching may occur in situations where evolutionary branching driven by direct disruptive selection would not occur. Hence, development can lead to phenotypic and genotypic diversification, even with a constant single-peak fitness landscape in geno-phenotype or phenotype space. Consequently, phenotypic and genotypic diversification might arise from the evolution of development.

Figure 2.

Development-driven diversification. Development may drive phenotypic and genotypic diversification in two ways. First, by shifting path peaks (a–f). Second, by creating path peaks (g–l). The evolutionary dynamics are given as in the legend of Figure 1. As c changes, development g changes, but total developmental bias $\mathrm{d}x/\mathrm{d}y$ does not change, so the evolutionary equilibria remain the same (blue); however, the admissible evolutionary equilibria do change (open and solid circles) and so does the outcome even though the fitness landscape is constant and single peaked.

Our analysis partly substantiates a classic explanation of punctuated equilibria in terms of developmental constraints, an explanation that has been dismissed in the light of previous theory. Punctuated equilibria refer to long periods of morphological stasis alternated by periods of rapid change as observed in the fossil record (Eldredge & Gould, 1972; Hunt et al., 2015). The dominant explanation for punctuated equilibria is that fitness landscapes remain relatively invariant for long periods so there is stabilizing selection for long periods at an unconstrained optimum, and fitness landscapes subsequently undergo substantial change particularly during adaptive radiations triggering directional selection (Charlesworth et al., 1982; Walsh & Lynch, 2018). A classic alternative explanation is that (a) developmental constraints prevent change for long periods and (b) revolutions in the developmental program bring sudden change (Gould, 1980). The ability of constraints to prevent change has been thought to be refuted by the common observation of selection response, sometimes under seeming developmental constraints (Beldade et al., 2002), while revolutions in the developmental program have been refuted on other grounds (Charlesworth et al., 1982; Eldredge et al., 2005; Klingenberg, 2010). Our analysis substantiates the view that developmental constraints may prevent change for long periods, despite possible selection response: indeed, even if developmental constraints halt a population at a path peak, selection response is possible since path peaks may change, either by change in the path (e.g., by change in the genotype, environment, or how the phenotype is constructed) or the landscape (e.g., by artificial selection as in Beldade et al., 2002). Thus, available evidence does not necessarily rule out that developmental constraints can prevent change for long periods in some directions. Moreover, punctuated equilibria could arise by change in the developmental map, rather than only by change of the fitness landscape. In particular, gradual change in the developmental map as the genotype, phenotype, or environment gradually evolve may yield sudden path-peak creation events (sudden since such events are bifurcations, which are sudden by definition; Figure 2g–l; see also Figure 3.4 of Metz, 2011). Path-peak creation events yielding sudden phenotypic diversification might then generate punctuated equilibria, arising by gradual evolution of the developmental map, without revolutions in the developmental program. Furthermore, by shifting the admissible path, the evolution of development may similarly lead to the crossing of valleys in a multi-peak fitness landscape in geno-phenotype space if a path peak moves from a landscape peak to another.

Niche construction reshapes the path

Niche construction—whereby organisms affect their environment—has been suggested to be a developmental factor having major evolutionary effects (Laland et al., 2014, 2015; Odling-Smee et al., 1996). How does niche construction affect evolution? We find that niche construction has a dual effect on adaptation. First, niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}$⁠) reshapes the fitness landscape as previously known (Laland et al., 2014, 2015; Odling-Smee et al., 1996). Second, niche construction reshapes the path on the landscape by altering absolute genetic constraints.

To see this, we relax the assumption that niche construction is absent (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}\ne \mathbf{0}$⁠), while maintaining our assumptions of (ii) no social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}=\mathbf{0}$⁠), and (iii) no exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau =\mathbf{0}$⁠). The evolutionary dynamics of the resident phenotype and genotype are then given by

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{z}}}}{\mathrm{d}\tau }=\iota {\mathbf{H}}_{\mathbf{z}}\frac{\delta w}{\delta \mathbf{z}},}\end{array}}$$

(13)

where ${\mathbf{H}}_{\mathbf{z}}$ guarantees that the developmental constraint (1) is met at all times. Equation (13) depends on the total immediate selection gradient of the phenotype and genotype

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\delta w}{\delta \mathbf{z}}=\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathbf{z}}+\frac{\mathrm{\partial }{\mathit{\epsilon }}^{⊺}}{\mathrm{\partial }\mathbf{z}}\frac{\mathrm{\partial }w}{\mathrm{\partial }\mathit{\epsilon }},}\end{array}}$$

(14)

which measures total immediate selection considering environmental constraints but not developmental constraints. So, selection pushes uphill, no longer on the fitness landscape that ignores environmental constraints, but on the reshaped total immediate fitness landscape that considers such constraints (Figure 3a–e). Equation (14) shows that such reshaping of the fitness landscape is done by the interaction of niche construction and environmental sensitivity of selection (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathit{\epsilon }$⁠) (a term coined by Chevin et al., 2010 for a more specific notion). The total immediate fitness landscape can have new peaks and valleys relative to the fitness landscape. Additionally, Equation (13) depends on ${\mathbf{H}}_{\mathbf{z}}$⁠, which still has the form in Equation (9) so it is still singular but now depends on the interaction of niche construction and plasticity (⁠$\mathrm{\partial }{\mathbf{x}}^{⊺}/\mathrm{\partial }\mathit{\epsilon }$⁠). Thus, development determines the path on the total immediate fitness landscape, but now this path is reshaped by the interaction of niche construction and plasticity (Figure 3f–j). Consequently, niche construction reshapes the fitness landscape by interacting with environmental sensitivity of selection, and reshapes the admissible path on the landscape by interacting with phenotypic plasticity. Thus, niche construction affects both selection and development, where selection and development jointly determine the evolutionary outcomes if there are no absolute mutational constraints.

Figure 3.

Niche construction reshapes the fitness landscape and the evolutionary path. a–e, There is niche construction by the phenotype x which non-linearly affects the environmental trait $\epsilon$⁠, which in turn increases fitness. a, Fitness landscape vs. phenotype, genotype, and environment. Slices for constant environment (yellow) and for the environmental constraint (green) are, respectively, shown in b and c. The interaction of niche construction (⁠$\mathrm{\partial }\epsilon /\mathrm{\partial }x$⁠) and environmental sensitivity of selection (⁠$\mathrm{\partial }w/\mathrm{\partial }\epsilon$⁠) reshapes the landscape (from b to c) by affecting total immediate selection $\delta w/\delta \mathbf{z}$ but it does not affect genetic covariation ${\mathbf{H}}_{\mathbf{z}}$⁠. f–j, There is niche construction by the genotype y which linearly affects the environmental trait $\epsilon$⁠, which in turn increases the phenotype. The interaction of niche construction (⁠$\mathrm{\partial }\epsilon /\mathrm{\partial }y$⁠) and plasticity (⁠$\mathrm{\partial }x/\mathrm{\partial }\epsilon$⁠) reshapes the path (from g to h). This interaction affects genetic covariation ${\mathbf{H}}_{\mathbf{z}}$ but does not affect total immediate selection $\delta w/\delta \mathbf{z}$⁠. In all panels, evolutionary change satisfies $\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau ={\mathbf{H}}_{\mathbf{z}}\delta w/\delta \mathbf{z}$⁠, where $\delta w/\delta \mathbf{z}=\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}+(\mathrm{\partial }\epsilon /\mathrm{\partial }\mathbf{z})(\mathrm{\partial }w/\mathrm{\partial }\epsilon )$ and ${\mathbf{H}}_{\mathbf{z}}$ is given in Figure 1. In a–e,$\mathrm{d}x/\mathrm{d}y=\mathrm{\partial }x/\mathrm{\partial }y$⁠, whereas in f–j,$\mathrm{d}x/\mathrm{d}y=\mathrm{\partial }x/\mathrm{\partial }y+(\mathrm{\partial }\epsilon /\mathrm{\partial }y)(\mathrm{\partial }x/\mathrm{\partial }\epsilon )$⁠.

Equation (13) could also be seen as corresponding to an extended Lande equation describing genotypic and phenotypic evolution where environmental traits are not explicitly considered in a quantitative genetics analysis. Now, Equation (13) generally depends on the resident environment $\overline{\mathit{\epsilon }}$⁠, which evolves due to niche construction. Provided that the environmental map is known or assumed, the environment is then known so Equation (13) still yields a dynamically sufficient description of long-term phenotypic evolution despite environmental evolution. Yet, if one is interested in understanding how the evolution of the environment proceeds on the fitness landscape, an equation describing the evolution of the environment in gradient form would be needed. González-Forero (2021) derives such an equation, showing it has the form of a further extended mechanistic Lande equation describing the evolutionary dynamics of the geno-envo-phenotype with an associated mechanistic additive genetic covariance matrix that is also necessarily singular for the same reasons given above.

Social development alters the outcome

Extra-genetic inheritance, including social learning and epigenetic inheritance, has also been suggested to be a developmental factor having major evolutionary effects (Baldwin, 1896; Laland et al., 2014, 2015). Current understanding is that extra-genetic inheritance has only transient evolutionary effects if elements so inherited are transmitted over only a few generations (Beltran et al., 2020; Charlesworth et al., 2017; Jablonka & Lamb, 2014; Walsh & Lynch, 2018). For instance, the permanent evolutionary effects discussed by Day & Bonduriansky (2011) and Bonduriansky & Day (2018) involve the co-evolution of genetic and extra-genetic elements, which requires transmission of extra-genetic elements over multiple generations. In contrast, we find that extra-genetic inheritance can have permanent evolutionary effects even without transmission of extra-genetic elements over multiple generations, because it changes development and so the absolute genetic constraints and the admissible evolutionary path.

Extra-genetic inheritance can be seen as a particular form of social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}\ne \mathbf{0}$⁠), where the developed phenotype depends on the phenotype and genotype of social partners. For instance, acquiring a skill or epimark from a role model or a parent, respectively, involves an interaction between the recipient and donor, and these extra-genetic elements affect the subsequent phenotype of the individual (although recall we assume that social interactions are only with non-relatives for simplicity). Social development can also mechanistically describe both indirect genetic effects, where the genes of social partners affect an individual’s phenotype (Moore et al., 1997), and social phenotypic effects that need not stem from genes (e.g., if the individual’s developed traits depend on social partners’ environmentally induced developed traits). How does social development, including extra-genetic inheritance, affect evolution? We find that, since it alters the developmental map, social development alters the evolutionary outcome (i.e., the admissible evolutionary equilibria) by altering the admissible evolutionary path without necessarily altering the evolutionary equilibria themselves. Thus, extra-genetically inherited elements can have permanent evolutionary effects without being transmitted over many generations. What is needed is that the ability to acquire the (possibly evolving) extra-genetic element persists over many generations so that individuals have a consistently socially altered development over many generations. To see this, consider the following.

Social development introduces a complication to evolutionary analysis as follows: With social development, the developed phenotype depends on the frequency of other phenotypes and genotypes and even if a phenotype is fixed in the population, offspring may develop a different phenotype if the fixed parental phenotype developed in a different social context (i.e., the phenotype is frequency dependent and may not breed true despite clonal reproduction). We handle this complication with the following notions. We say a genotype, phenotype, and environment ${\mathbf{m}}^{\ast \ast }=({\mathbf{x}}^{\ast \ast };{\mathbf{y}}^{\ast \ast };{\mathit{\epsilon }}^{\ast \ast })$ is a socio-devo equilibrium if and only if an individual with such genotype and environment developing in the context of individuals having such genotype, phenotype, and environment develops their same phenotype. A socio-devo equilibrium is socio-devo stable if and only if socio-devo dynamics are locally stable to small perturbations in the phenotype (González-Forero, 2021). A phenotype in socio-devo stable equilibrium breeds true when fixed. Although social development can yield socio-devo unstable phenotypes, our analyses apply only to socio-devo stable phenotypes (which are guaranteed to be so if all the eigenvalues of $\mathrm{d}\mathbf{x}/\mathrm{d}{\overline{\mathbf{x}}}^{⊺}$ have absolute value strictly less than one, which we assume).

Thus, for socio-devo stable phenotypes, their evolutionary dynamics are described by the following equations. Allowing for niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}\ne \mathbf{0}$⁠) and social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}\ne \mathbf{0}$⁠), but assuming (iii) that there is no exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau =\mathbf{0}$⁠), then the evolutionary dynamics of the resident phenotype and genotype are given by

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{z}}}}{\mathrm{d}\tau }=\iota {\mathbf{L}}_{\mathbf{z}}\frac{\delta w}{\delta \mathbf{z}}=\iota {\mathbf{L}}_{\mathbf{z}\mathbf{y}}\frac{\mathrm{d}w}{\mathrm{d}\mathbf{y}},}\end{array}}$$

(15)

where both ${\mathbf{L}}_{\mathbf{z}}$ and ${\mathbf{L}}_{\mathbf{z}y}$ guarantee that the developmental constraint (1) is met at all times. This equation now depends on the mechanistic additive socio-genetic cross-covariance matrix of the geno-phenotype (L for legacy)

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{L}}_{\mathbf{z}}=\mathrm{c}\mathrm{o}\mathrm{v}[{\mathbf{b}}_{\mathbf{z}}^{\mathrm{s}},{\mathbf{b}}_{\mathbf{z}}]=\frac{\mathrm{s}\mathbf{z}}{\mathrm{s}{\mathbf{y}}^{⊺}}{\mathbf{H}}_{\mathbf{y}}\frac{\mathrm{d}{\mathbf{z}}^{⊺}}{\mathrm{d}\mathbf{y}},}\end{array}}$$

(16)

where the stabilized mechanistic breeding value of $\mathbf{z}$ is defined as

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{b}}_{\mathbf{z}}^{\mathrm{s}}={\displaystyle \overline{\mathbf{z}}}+\frac{\mathrm{s}\mathbf{z}}{\mathrm{s}{\mathbf{y}}^{⊺}}(\mathbf{y}-{\displaystyle \overline{\mathbf{y}}}).}\end{array}}$$

(17)

The matrix $\mathrm{s}\mathbf{z}/\mathrm{s}{\mathbf{y}}^{⊺}$ is the stabilized total derivative of $\mathbf{z}$ with respect to ${\mathbf{y}}^{⊺}$⁠, which measures the total effects of ${\mathbf{y}}^{⊺}$ on $\mathbf{z}$ after the effects of social development have stabilized in the population; $\mathrm{s}\mathbf{z}/\mathrm{s}{\mathbf{y}}^{⊺}$ reduces to $\mathrm{d}\mathbf{z}/\mathrm{d}{\mathbf{y}}^{⊺}$ if development is not social. The ${\mathbf{L}}_{\mathbf{z}}$-matrix generalizes ${\mathbf{H}}_{\mathbf{z}}$ by including the effects of social development, so if development is non-social ${\mathbf{L}}_{\mathbf{z}}$ reduces to ${\mathbf{H}}_{\mathbf{z}}$ (⁠${\mathbf{L}}_{\mathbf{z}}$ is analogous to another generalization of the $\mathbf{G}$-matrix used in the indirect genetic effects literature; Moore et al., 1997). As for ${\mathbf{H}}_{\mathbf{z}}$⁠, the matrix ${\mathbf{L}}_{\mathbf{z}}$ is always singular because the geno-phenotype contains the genotype (i.e., because $\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}\mathbf{y}$ has fewer rows than columns). Hence, also under social development, evolutionary stasis can generally occur with persistent total immediate selection on the phenotype and genotype. In turn, the second equality in Equation (15) depends on the additive socio-genetic cross-covariance matrix between the geno-phenotype and genotype

$${\displaystyle \begin{array}{l}{\displaystyle {\mathbf{L}}_{\mathbf{z}\mathbf{y}}=\mathrm{c}\mathrm{o}\mathrm{v}[{\mathbf{b}}_{\mathbf{z}}^{\mathrm{s}},{\mathbf{b}}_{\mathbf{y}}]=\frac{\mathrm{s}\mathbf{z}}{\mathrm{s}{\mathbf{y}}^{⊺}}{\mathbf{H}}_{\mathbf{y}},}\end{array}}$$

(18)

which generalizes the corresponding matrix ${\mathbf{H}}_{\mathbf{z}\mathbf{y}}$ we had in Equation (12) to consider social development. Similarly to ${\mathbf{H}}_{\mathbf{z}\mathbf{y}}$⁠, this matrix ${\mathbf{L}}_{\mathbf{z}\mathbf{y}}$ is non-singular if ${\mathbf{H}}_{\mathbf{y}}$ is non-singular because the genotype is developmentally independent by assumption. Thus, evolutionary equilibria still occur when total genotypic selection vanishes (⁠$\mathrm{d}w/\mathrm{d}\mathbf{y}=\mathbf{0}$⁠) if there are no absolute mutational constraints.

The singularity of ${\mathbf{L}}_{\mathbf{z}}$ entails that social development, including extra-genetic inheritance, may alter both the evolutionary equilibria (by altering $\mathrm{d}w/\mathrm{d}\mathbf{y}$⁠) and the admissible path (by altering ${\mathbf{L}}_{\mathbf{z}\mathbf{y}}$⁠; Figure 4). Note that $\overline{y}$ in Figure 4a,b can be interpreted as being an extra-genetically acquired element that modifies development, but such acquired effect is not transmitted to future generations as the individual transmits the y it has. Yet, such extra-genetically acquired element alters the developmental constraint and evolutionary equilibria, thus altering the evolutionary outcome relative to that without extra-genetic inheritance (Figure 1a,d,g,j). Alternatively, social development may not affect the evolutionary equilibria if total genotypic selection $\mathrm{d}w/\mathrm{d}\mathbf{y}$ is independent of social partners’ geno-phenotype $\overline{\mathbf{z}}$ , but social development may still affect the admissible evolutionary path given by the developmental constraint in Equation (1) (by altering only ${\mathbf{L}}_{\mathbf{z}\mathbf{y}}$⁠; Figure 5). Similarly, $\overline{y}$ in Figure 5a,b can also be interpreted as being an extra-genetically acquired element that modifies development, but it is not transmitted to future generations and still alters the evolutionary outcome relative to that without extra-genetic inheritance (Figure 1b,e,h,k). Indeed, in both cases, social development (e.g., extra-genetic inheritance) affects the evolutionary path and so the admissible evolutionary equilibria and hence the evolutionary outcome. Mathematically, this results because ${\mathbf{L}}_{\mathbf{z}}$ is singular, so changing the path generally changes the outcome. Thus, because ${\mathbf{L}}_{\mathbf{z}}$ is singular, extra-genetic inheritance can have permanent evolutionary effects even if extra-genetically acquired elements are not transmitted to other individuals.

Figure 4.

Social development alters the evolutionary outcome. a, b, There is social development, particularly social developmental bias from the genotype (⁠$\mathrm{\partial }x/\mathrm{\partial }\overline{y}\ne 0$⁠). Here social development introduces a non-linearity relative to the developmental map of Figure 1a,d,g,j, which changes the evolutionary equilibria (blue), admissible path (red), and outcome (solid dot). c, d, As in a, b, but there is social developmental bias from the phenotype (⁠$\mathrm{\partial }x/\mathrm{\partial }\overline{x}\ne 0$⁠). The developmental map is that of a, b replacing $\overline{y}$ with $\overline{x}$⁠. This also changes the evolutionary equilibria, admissible path, and outcome. In all panels, evolutionary change satisfies $\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau ={\mathbf{L}}_{\mathbf{z}}\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}$ , where $w(x,y)=\exp (-x^2-y^2)$ and ${\mathbf{L}}_{\mathbf{z}}=(\mathrm{s}\mathbf{z}/\mathrm{s}y)H_y(\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}y)=(L_x,L_{xy};L_{yx},L_y)$⁠, with $\mathrm{s}\mathbf{z}/\mathrm{s}y=(\mathrm{s}x/\mathrm{s}y;1)$⁠. For a, b, $\mathrm{s}x/\mathrm{s}y=\mathrm{d}x/\mathrm{d}y+\mathrm{d}x/\mathrm{d}\overline{y}$ (from Layer 5, Equation 2a in González-Forero, 2021). Similarly, for c–d, $\mathrm{s}x/\mathrm{s}y=(1-\mathrm{\partial }x/\mathrm{\partial }\overline{x}{)}^{-1}\mathrm{d}x/\mathrm{d}y$⁠, with $\mathrm{d}{\mathbf{z}}^{⊺}/\mathrm{d}y$ and $\mathrm{s}x/\mathrm{s}y$ evaluated at the socio-devo equilibrium x^** that solves $x^{\ast \ast }=g(y,x^{\ast \ast })$⁠. Since development is social, a fitness landscape visualization as in Figure 1a is not possible. In c, the stream plot (gray arrows) is only locally accurate around the developmental constraint where the socio-devo equilibrium holds. The 1/2 in the developmental constraint is used to meet the assumption $|\mathrm{d}x/\mathrm{d}\overline{x}|<1$ but does not affect the points made.

Figure 5.

Social development can alter the evolutionary outcome without altering evolutionary equilibria. The plots are as in Figure 4, except that now the developmental maps are a linear modification of those in Figure 1b,e,h,k. Relative to that figure, social development does not alter total developmental bias $\mathrm{d}x/\mathrm{d}y$ so it does not affect evolutionary equilibria (blue). However, social development affects the admissible path (red) and so the evolutionary outcome (solid dot). Note L_x can be negative in b as it gives the covariance between the stabilized mechanistic breeding value of x and the mechanistic breeding value of x (in contrast, H_x gives the variance of the mechanistic breeding value of x, which is non-negative). The 4/5 in the developmental constraint is used to meet the assumption $|\mathrm{d}x/\mathrm{d}\overline{x}|<1$ but does not affect the points made.

When social development alters only the admissible path, it does so because of social developmental bias (⁠$\mathrm{\partial }{\mathbf{x}}^{⊺}/\mathrm{\partial }\overline{\mathbf{z}}$⁠; e.g., extra-genetic inheritance and indirect-genetic effects), and/or the interaction of social niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\overline{\mathbf{z}}$⁠) and plasticity (⁠$\mathrm{\partial }{\mathbf{x}}^{⊺}/\mathrm{\partial }\mathit{\epsilon }$⁠) (from Layer 3, Equation 3 of González-Forero, 2021).

Plastic response alters the path

Research on the evolutionary effects of plasticity has intensified, with a particular focus on whether adaptive plastic change precedes genetic evolution (Waddington, 1961; West-Eberhard, 2003). How does plasticity affect evolution? We find that plasticity also has a dual evolutionary effect. First, we saw that plasticity (⁠$\mathrm{\partial }{\mathbf{x}}^{⊺}/\mathrm{\partial }\mathit{\epsilon }$⁠) alters the evolutionary path by interacting with niche construction (i.e., endogenous environmental change; $\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}$⁠). This is because niche construction alters the environment, which through plasticity alters the developed phenotype and so the evolutionary path (by altering ${\mathbf{L}}_{\mathbf{z}}$⁠). Second, plasticity alters the evolutionary path by interacting with exogenous environmental change (e.g., eutrophication or climate change caused by other species; $\mathrm{\partial }\epsilon /\mathrm{\partial }\tau$⁠). This is similarly because exogenous environmental change alters the environment, which through plasticity alters the developed phenotype and so the evolutionary path (but not by directly altering ${\mathbf{L}}_{\mathbf{z}}$⁠). With exogenous environmental change, the evolutionary dynamics comprise selection response and exogenous plastic response (see also Chevin et al., 2010), the latter of which modifies the developmental map and thus the evolutionary path as evolutionary time advances. Through exogenous plastic response, the evolutionary effects of development go beyond modulating genetic covariation. Exogenous plastic response has some ability to increase adaptation separate from selection response via selective development, although this ability is limited in various ways.

To see this, assume now that there is niche construction (⁠$\mathrm{\partial }{\mathit{\epsilon }}^{⊺}/\mathrm{\partial }\mathbf{z}\ne \mathbf{0}$⁠), social development (⁠$\mathrm{d}{\mathbf{x}}^{⊺}/\mathrm{d}\overline{\mathbf{z}}\ne \mathbf{0}$⁠), and exogenous environmental change (⁠$\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau \ne \mathbf{0}$⁠). Then, the evolutionary dynamics of the resident phenotype and genotype are given by

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{\displaystyle \overline{\mathbf{z}}}}{\mathrm{d}\tau }=\iota {\mathbf{L}}_{\mathbf{z}}\frac{\delta w}{\delta \mathbf{z}}+\frac{\mathrm{s}\mathbf{z}}{\mathrm{s}{\mathit{\epsilon }}^{⊺}}\frac{\mathrm{\partial }\mathit{\epsilon }}{\mathrm{\partial }\tau }=\iota {\mathbf{L}}_{\mathbf{z}\mathbf{y}}\frac{\mathrm{d}w}{\mathrm{d}\mathbf{y}}+\frac{\mathrm{s}\mathbf{z}}{\mathrm{s}{\mathit{\epsilon }}^{⊺}}\frac{\mathrm{\partial }\mathit{\epsilon }}{\mathrm{\partial }\tau }.}\end{array}}$$

(19)

The matrix $\mathrm{s}\mathbf{z}/\mathrm{s}{\mathit{\epsilon }}^{⊺}$ describes what we call stabilized plasticity of the phenotype and genotype. Stabilized plasticity includes the effects of socio-devo stabilization and reduces to total plasticity $\mathrm{d}\mathbf{z}/\mathrm{d}{\mathit{\epsilon }}^{⊺}$ if development is non-social. The first term in Equation (19) is selection response (⁠$\iota {\mathbf{L}}_{\mathbf{z}}\delta w/\delta \mathbf{z}$⁠) whereas the second term is exogenous plastic response (⁠$(\mathrm{s}\mathbf{z}/\mathrm{s}{\mathit{\epsilon }}^{⊺})(\mathrm{\partial }\mathit{\epsilon }/\mathrm{\partial }\tau )$⁠), which makes the developmental map evolve as the environment exogenously changes (Figure 6). Equation (19) shows that exogenous plastic response can induce selection response, which allows for plastic change to preceed genetic change (Waddington, 1961; West-Eberhard, 2003). For example, if the population is at a path peak so there is no selection response (⁠${\mathbf{L}}_{\mathbf{z}}\delta w/\delta \mathbf{z}=\mathbf{0}$⁠) and there is exogenous plastic response, this generally changes the path thus inducing future selection response (e.g., in Figure 6a, selection response is absent at the two path peaks at x = 1/2 at $\tau =0$⁠, but environmental change and plasticity change the path thereby creating selection response at these points and the path peaks move approximately to x = 0 and x = 0.18 at $\tau =1000$ in Figure 6c). In the idealized case where development is selective in the sense that the developmental map is a function of the selection gradient, stabilized plasticity can depend on the selection gradient. With selective development, exogenous plastic response can induce evolution either towards or away landscape peaks depending on the developmental constraints and initial conditions, as the population is trapped in local peaks of the total fitness landscape of the genotype (Figure 6). Exogenous plastic response with selective development has a limited ability to induce evolution toward a landscape peak since plasticity must depend on developmental-past environments, which induces a developmental lag in the plastic response (Supplementary Material section S2). Additionally, exogenous plastic response with selective development may induce evolution away a landscape peak when the rate of environmental change switches sign (Supplementary Material section S2).

Figure 6.

Exogenous plastic response changes the evolutionary path. We let there be exogenous environmental change and plasticity, but no social development, so ${\mathbf{L}}_{\mathbf{z}}$ reduces to ${\mathbf{H}}_{\mathbf{z}}$ and no niche construction, so plasticity does not affect ${\mathbf{H}}_{\mathbf{z}}$ but only the exogenous plastic response. For illustration, development is selective by letting plasticity equal the selection gradient at the optimal phenotype, $\mathrm{\partial }x/\mathrm{\partial }\epsilon =\mathrm{\partial }w(0,y)/\mathrm{\partial }y$⁠. a–l, Exogenous environmental change induces exogenous plastic response which raises part of the path on the landscape. m, Resulting evolutionary dynamics showing that selective development induces some evolutionary trajectories to converge to the landscape peak, whereas others do not by getting trapped in another path peak due to their initial condition. Evolutionary change satisfies $\mathrm{d}\overline{\mathbf{z}}/\mathrm{d}\tau ={\mathbf{H}}_{\mathbf{z}}\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}+(\mathrm{d}\mathbf{z}/\mathrm{d}\epsilon )(\mathrm{\partial }\epsilon /\mathrm{\partial }\tau )$⁠, where $\mathrm{d}\mathbf{z}/\mathrm{d}\epsilon =(\mathrm{d}x/\mathrm{d}\epsilon ;\mathrm{d}y/\mathrm{d}\epsilon )=(\mathrm{\partial }x/\mathrm{\partial }\epsilon ;0)$ and ${\mathbf{H}}_{\mathbf{z}}$ is given in Figure 1. Fitness is $w(x,y)=\exp (-x^2-y^2)$⁠.

Development enables negative senescence

Thus far, we have considered evolutionary effects of development that occur even without considering that development takes time. Indeed, in our illustrative examples we have let age structure be collapsed so development happens instantaneously. We now show that explicit consideration of age progression in development may help explain an otherwise puzzling ageing pattern.

Individuals may show senescence as they age, that is, decreasing survival or fertility after the onset of reproduction. Leading hypotheses for the evolution of senescence are based on the fact that the forces of selection decline with age (⁠${\varphi }_a$ and ${\pi }_a$ in Equation (4), and Equations S8 in the Supplementary Material; Hamilton, 1966). For example, if a mutation has a beneficial effect on survival or fertility early in life and a pleiotropic, deleterious effect of similar magnitude late in life, the mutation would be favoured (Medawar, 1952; Williams, 1957). The universality of declining forces of selection has suggested that senescence should be universal (Hamilton, 1966), at least for organisms with a bodily reset at conception (Lehtonen, 2020). However, empirical research has increasingly reported organisms with negative senescence (Jones et al., 2014b), that is, where survival or fertility increase with age after the onset of reproduction, and models have suggested that this may stem from indeterminate growth (Lehtonen, 2020; Vaupel et al., 2004).

Consistently with the latter, we find that developmental propagation—whereby an early-acting mutation has a large phenotypic effect late in life—can drive the evolution of negative senescence. For example, if a mutation in genotypic trait i has a deleterious effect on survival or fertility at an early age a, so that $\delta w_a/\delta y_{ia}<0$⁠, and a pleiotropic, beneficial effect on survival or fertility of similar magnitude at a later age j>a, so that $\delta w_j/\delta x_{kj}>0$ for phenotype k, then total immediate selection is weaker at the later age because of declining selection forces (i.e., $|\delta w_a/\delta y_{ia}|\ge |\delta w_j/\delta x_{kj}|$⁠). Yet, crucially, such late-beneficial mutation is still favoured if total genotypic selection on it is positive, that is, if

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}w}{\mathrm{d}{y}_{ia}}=\frac{\delta w_a}{\delta y_{ia}}+\frac{\mathrm{d}{x}_{kj}}{\mathrm{d}{y}_{ia}}\frac{\delta w_j}{\delta x_{kj}}>0}\end{array}}$$

(20)

(Supplementary Material section S3), where the effect of the early mutation on the late phenotype is

$${\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_{kj}}{\mathrm{d}{y}_{ia}}=\frac{\delta {\mathbf{x}}_{a+1}^{⊺}}{\delta y_{ia}}\left({\displaystyle \prod _{i=a+1}^{j-2}}\frac{\delta {\mathbf{x}}_{i+1}^{⊺}}{\delta {\mathbf{x}}_i}\right)\frac{\delta x_{kj}}{\delta {\mathbf{x}}_{j-1}}}\end{array}}$$

(21)

(from Equation C17 of González-Forero, 2021). This effect of the early mutation on the late phenotype may increase (or decrease) with the phenotype age j because of developmental propagation of phenotypic effects as age advances, a propagation that is described by the term in parentheses in Equation (21): in particular, indeterminate growth allows for indefinite developmental propagation. So total genotypic selection may increase with the age j at which the pleiotropic effect is expressed. Consequently, a mutation with late beneficial effects may be favoured despite early deleterious effects. Thus, developmental propagation may enable selection to shape traits at advanced age, allowing for adaptive evolution of increasing survival or fertility as ageing proceeds despite decreasing selection forces after the onset of reproduction. This suggests that allowing for developmental propagation in models with explicit development in which negative senescence fails to evolve (Nelson & Masel, 2017) might still enable negative senescence (perhaps via cell maintenance).

Discussion

Evolution by natural selection can be understood as the climbing of a fitness landscape under some assumptions. This was first established for genetic evolution (Wright, 1937) and was later established for phenotypic evolution under the tacit assumption that there is negligible genetic evolution (Barton et al., 2017; Hill, 2017; Lande, 1976, 1979; Walsh & Lynch, 2018). However, long-term phenotypic evolution involves non-negligible genetic evolution. We have shown that understanding long-term phenotypic evolution as the climbing of a fitness landscape requires that genotypic and phenotypic evolution are considered simultaneously. Doing so necessarily makes the associated constraining matrix (⁠$\mathbf{H}$⁠, $\mathbf{L}$⁠, or $\mathbf{G}$⁠) singular because the phenotype is related to the genotype via development. Thus, long-term phenotypic evolution necessarily involves absolute genetic constraints, where the evolutionary path is given by the developmental process that associates genes to phenotype.

Consequently, we have shown that development affects evolution by determining the evolutionary path. This sharpens the basic principle of adaptation as the climbing of the fitness landscape (Lande, 1979) (sensu Dieckmann & Law, 1996), such that evolutionary outcomes occur at best at peaks on the path rather than on the landscape. Developmental constraints thus stop evolution at path peaks that can happen mid-way on a fitness slope. Hence, selection and development jointly define the evolutionary outcomes if there are no absolute mutational constraints. Selective development, if feasible, can bring some of the path closer to the landscape peak by means of exogenous plastic response. Our findings substantiate several although not all intuitions previously given in relation to the role of development in evolution (Gould & Lewontin, 1979; Jablonka & Lamb, 2014; Laland et al., 2014, 2015; Pigliucci, 2007; Waddington, 1959; West-Eberhard, 2003). In particular, seemingly disparate development-related factors previously suggested to be evolutionarily important (Laland et al., 2014, 2015)—namely phenotypic plasticity (West-Eberhard, 2003), niche construction (Odling-Smee et al., 1996), extra-genetic inheritance (Jablonka & Lamb, 2014), and developmental bias (Arthur, 2004)—variously alter development, and so the absolute genetic constraints and consequently the evolutionary outcome. Our analysis offers answers to various major questions, namely the origin of diversity, the punctuated fossil record, the paradox of stasis, the rarity of stabilizing selection, and the origin of negative senescence. Our observations entail that an understanding of development is instrumental for evolutionary understanding: change in development alone changes evolutionary outcomes by changing absolute socio-genetic constraints (⁠${\mathbf{L}}_{\mathbf{z}}$⁠), even if direct (⁠$\mathrm{\partial }w/\mathrm{\partial }\mathbf{z}$⁠) or total (⁠$\mathrm{d}w/\mathrm{d}\mathbf{z}$ or $\mathrm{d}w/\mathrm{d}\mathbf{y}$⁠) selection remain constant (Figure 5). These observations call for empirical estimation of developmental maps, for instance, via the rapidly developing methods to estimate dynamic equations from data (Brunton et al., 2016; Ghadami & Epureanu, 2022, and papers in the special issue, Schmidt & Lipson, 2009). Overall, our analysis finds that development has major evolutionary consequences.

Data availability

The Mathematica notebook used to generate the figures is accessible from the data repository https://doi.org/10.5281/zenodo.7180175. All other data relevant to this publication is in the main paper or supplementary information.

Author contributions

M.G.F. conceived and performed the research and wrote the paper.

Funding statement

This work was funded by an European Research Council Consolidator Grant to A. Gardner (grant no. 771387), by the School of Biology of the University of St Andrews, and by a John Templeton Foundation grant to K.N. Laland and T. Uller (grant ID 60501).

Conflict of interest: The authors declare no conflict of interest.

Acknowledgments

I thank A. Gardner for extensive support throughout this project, by discussing, reading, and commenting on the many drafts, and offering interpretation, advice, and funding. I thank N.W. Bailey, M.W. Feldman, K.N. Laland, R. Lande, L.C. Mikula, A.J. Moore, and M.B. Morrissey for comments on previous versions of the manuscript, and J.A.J. Metz, I. Salazar-Ciudad, and D.M. Shuker for discussion. I thank J. Masel, M. Pavlicev, and two anonymous reviewers for detailed criticism.