The canonical equation of adaptive dynamics for Mendelian diploids and haplo-diploids

1. Introduction

For context, it is useful to distinguish between micro-, meso- and macro-evolution. The term micro-evolution customarily refers to changes in gene frequencies on a population dynamical time-scale. We will refer to the evolution of quantitative traits through the repeated substitution of novel mutants, including the splitting of the evolutionary path into separate evolutionary lines, as meso-evolution. The term macro-evolution then becomes restricted to larger scale changes such as anatomical innovations, where one cannot even speak in terms of a fixed set of traits.

Adaptive dynamics (AD) was devised as a mathematical framework for dealing with meso-evolution in ‘realistic’ ecological settings. It differs from more classical approaches to modelling evolutionary change, which generally assume constant fitnesses, by its focus on the population dynamical basis for those fitnesses, and hence on their inevitable change over evolutionary time. Of course, the greater realism at the ecological end is brought about by making different simplifying assumptions, this time genetically unrealistic ones. The main simplification is (i) separation of the population dynamical and mutational time scales. In order to concentrate on ecological aspects, unencumbered by the complexities of the genetic architecture and genotype to phenotype map, most AD research moreover assumes (ii) clonal inheritance.

One of the powerful tools of AD is its so-called canonical equation (CE), a differential equation describing how the prevailing trait vector changes over evolutionary time. The derivation of the CE is based on a subsequent further simplifying assumption: (iii) small mutational steps. As the speed at which a mutant substitutes is proportional to its effect, conditions (i) and (iii) will only rarely be met together. In §5, we give arguments why for small step sizes the time-scale separation need not be very strict.

The CE was first derived in [1] for well-mixed clonal populations with the simplest possible life histories. This was underpinned by a level of rigour sufficient to satisfy hard probabilists in [2] and [3]. The initial results were extended in [4] to, possibly spatially distributed, locally well-mixed populations with general life histories, with [5] giving a probabilistically rigorous underpinning for the case of simple age dependence. In this paper, we describe the extension to Mendelian diploids and haplo-diploids, once again on a physicist level of rigour. As it turns out, the CE for the simplest and most complex life histories differs only in a scalar factor summarizing how the intricacies of the life history feed through to the invasion probabilities of advantageous mutants. Mendelian diploidy brings in an additional factor 2, owing to the doubling of the number of mutant alleles per individual over a substitution.

Collet et al. [6] also derives a CE for Mendelian diploids, with full rigour. However, this in essence is a CE for a single locus trait, selected by a two-tiered ecology, within and between diploid bodies, the latter with the simplest possible life history. ([6] also lists the early applications of AD to Mendelian models.) We consider general life histories and phenotypic traits that are underlain by many loci.

In §2, we give a general heuristic derivation of the CE for clonal and for haploid and diploid Mendelian populations. In §3, we work out the details, and in §4 we consider haplo-diploids. The technicalities can be found in a suite of appendices. In the final §5, we discuss the strengths and weaknesses of the CE as a tool for evolutionary understanding.

2. Deriving the canonical equation for the textbook genetic scenarios

Mathematically, the CE is derived by taking two subsequent limits: (i) letting the system size K (and hence the average population size ) go to infinity (to make the community dynamics deterministic) and the mutation probability per birth event ɛ go to zero in such a manner that (a) ɛK ln(K) → 0 (to make the time for a substitution shorter than the time between mutations) and (b) ɛK exp(αK) → ∞ for sufficiently small positive α (to keep the population from going extinct on the time scale of the trait changes), followed by (ii) letting the mutational step sizes go to zero, all the while keeping the trait changes in view by rescaling time. Biologically, the CE is best seen as an approximation. From that perspective, it is expedient to express the result in the original time scale so that the basic biological parameters are kept in view.

Under the assumptions that K is large, ɛK ln(K) is small and ɛK exp(αK) is large for a sufficiently small positive α, that mutations are unbiased and that the environment as perceived by the individuals does not fluctuate (note that this implies a non-fluctuating resident population), the rate of change of a trait vector X can in the clonal case be expressed approximately as

2.1

The first term in square brackets after the second equals sign in (2.1) comes from multiplying the probability of a mutation per birth by the population birth rate, , a formula derived from the consistency relation (the average number of particles in a ‘reservoir’ equals the average entrance rate multiplied by the average residence time).

The second term comes from combining three approximation formulae to calculate the probability p that a mutant invades into the resident population, all coming from a branching process approximation for the invasion phase of the mutant dynamics, followed by averaging the product of the resulting expression and the mutational effect over the distribution g of Z.

2.2

2.3

2.4

2.5

2.6

In the derivation of (2.6), we have used the very strong interpretation of ‘unbiased’ that Z not only has mean zero but also is distributed symmetrically around that mean. Relaxing these assumptions gives an expression with considerably less appeal [1–3]. As relaxing them would make the following arguments less easy to follow while their essence stays the same, we have chosen to stick to the time-honoured simplification.

Moreover, in writing down (2.1), we have tacitly assumed that an invading mutant that makes it through the stochastic boundary layer, where its population dynamics can be approximated by a branching process, also makes it to a full substitution. This ‘invasion implies substitution’ rule presumably holds good for small mutational steps, away from population dynamical bifurcation points and evolutionarily singular strategies (characterized by ∂₁s(X|X) = 0). A hard proof is available for the case where the community dynamics allows a finite dimensional representation [21–23]. However, it looks as if with the right mathematical expertise the rule should be extendable to the required generality by combining the approaches in [24] and [25,26].

The two terms in square brackets in (2.1) connect to the simplifying assumptions in the following manner. The first term is contingent on the time-scale separation assumption. Only when substitutions do not interfere does the speed of trait movement become proportional to the number of invasion attempts per time unit. The derivation of the second term is contingent on the smallness of the mutational steps and also on the time-scale separation to produce the constant environment underlying the branching process calculations (which assume fixed probability distributions for the lifetime offspring numbers). More in particular, time-scale separation means that in between the negligibly short substitution events the environment is stationary and resident populations are genetically homogeneous. (More about the latter, seemingly unrealistic, consequence in §5.)

We now consider the Mendelian case. For chromosomal sex determination, we focus on the autosomes, deferring the allosomal contributions to §4.

Thanks to the genetic homogeneity of the residents, the argument in (2.1) extends seamlessly to haploids. The only difference from the clonal case is that thinking genetics points one to the concept of genotype to phenotype map and mutations that occur on multiple loci. However, when substitutions occur singly the latter multiplicity becomes phenotypically inconsequential.

Moving on to diploids, we first take a closer look at genotype to phenotype maps. The prevalent view in Evo-Devo nowadays is that the trait changes that AD attempts to model are mostly caused not by changes in the coding regions of genes but at their regulatory regions (e.g. [27,28]). The latter determine the activity of the genes in different parts of the body, at different times during development and under different micro-environmental conditions. This scenario allows us to look at the genotype as a sequence of vectors of expression levels, with a (a_i) a placeholder for the name of (an allele on) a generic locus. The genotype to phenotype map Φ transforms this sequence into phenotypic traits. It is from this perspective that one should judge the assumption of smallness of mutational steps: the influence of any specific regulatory site among its many colleagues tends to be relatively minor. And it is this perspective that allows us to assume that (iv) genotype to phenotype maps are smooth.

Lemma

([29], A. Pugliese 1996, personal communication). When there are no parental effects, smooth genotype to phenotype maps are locally additive, i.e. if at a the expression vector U_a mutates to U_A, then

2.7

Proof.

Without parental effects

2.8

Hence,

2.9

▪

In diploids, an invading mutant allele A practically always shares a body with a resident allele a and this aA reproduces through backcrossing with a resident aa. Hence, the allele population initially grows as clonally reproducing aAs (producing aas on the side), and its invasion fitness corresponds to the asymptotic average per capita growth of that clonal population in an environmental background provided exclusively by, also seemingly clonally reproducing (for homogeneous), residents. The invasion implies that substitution theorem also applies unchanged (thanks to the local additivity). However, after substitution the population consists of mutant homozygotes, making the resulting step in phenotype space twice as large as in the clonal and haploid cases. We thus conclude that the CE for Mendelian diploids reads

2.10

On the right-hand side, we have put first the ecologically determined number of resident haplotypes (sets of chromosomes as present in gametes), followed by the life-history statistics controlling the initial demographic stochasticity of allelic substitutions, followed by an expression quantifying the per birth mutational variability generated by the genetic architecture and genotype to phenotype map, to conclude with the selection gradient, summarizing the ecology's current tendency for filtering novel genetic variation. In Metz & Jansen [30], it is argued that times the second group of quantities precisely equals the effective population size from population genetics.

3. Filling in the details

It may seem that with (2.10) we are done. However, the devil is in the detail, to wit the calculation of s, T_s, T_r, and C. For background material on the ecological models covered, see [20,31] for the resident and [32,33] for the invader dynamics.

In principle also presents a problem, although only when the resident population fluctuates. This is even so in the clonal case as is not just a time average but a peculiarly weighted one (e.g. [34,35]). To keep things simple, we have confined the argument to non-fluctuating residents.

We first consider R₀, although this quantity does not appear in our (2.1)-derived list. The reason that it does not do so is that we wanted to write the CE in the form customary in the literature. However, for most structured populations R₀ is far easier to calculate than s. Therefore we might just as well in (2.1) and (2.10) drop T_r and replace s with R₀.

R₀ equals the dominant eigenvalue of the next-generation matrix L(Y|X) (or operator in the case of infinitely many birth states; e.g. [9]). For the clonal case, L is constructed by calculating from a model for the behaviour of individuals how many offspring in different birth states they produce on average, dependent on their own birth state.

Given the next-generation matrix, we can introduce two further quantities for later use: the stationary birth state distribution, i.e. normalized positive right eigenvector of L(X|X) going with R₀ = 1, , 1^TU = 1 with , and the corresponding reproductive values, i.e. co-normalized left eigenvector, , VU = 1. v_i equals the expected contribution of a newborn of type i to future birth rates.

How can we define R₀ for evolving sexual diploids? The first answer is that in hermaphrodites one can just add the numbers of offspring that individuals father and mother over their life (i.e. produce through the micro- and macro-gametic routes) and divide by 2. The factor 1/2 comes from the fact that in the Mendelian process each allele is only transferred with that probability. When hermaphrodites are born stochastically equal (i.e. their birth states have the same probability distribution),

3.1

3.2

Equations (3.1) and (3.2) also hold good for dioecious organisms, but with a twist. For later use, we note that we then can rewrite (3.1) and (3.2) by letting p_f and p_m denote the fractions of newly produced females and males and f₊ and m₊ the average lifetime numbers of offspring begotten by a female or male. Then f = p_ff₊ and m = p_mm₊, so that

3.3

3.4

3.5

Appendix C indicates how the preceding considerations extend in the presence of additional birth state distinctions.

Next we consider T_s, the quantity that had to be combined with to get the population birth rate.

3.6

3.7

3.8

, appearing as a coefficient in the approximation formula for the invasion probability p, is the most complicated beast. We first give its formula for the clonal case

3.9

To calculate for the Mendelian diploid case, we apply (3.9) to the alleles. We only consider the case with all offspring born stochastically equal.

3.10

For dioecious organisms, c = −1 (as Ek_fk_m = 0 and Ek_f = Ek_m = 1), and , , θ_f+ = σ_f+/f₊ and θ_m+ = σ_m+/m₊. The latter result is obtained from the following lemma together with the observation that E ₊ k_f = f₊ and E ₊ k_m = m₊, with E ₊ the expectation conditional on the parental sex.

Lemma.

Let the discrete non-negative random variable k be constructed from a random variable k₊ by setting k = k₊ with probability p and k = 0 with probability 1 − p, then Ek = pEk₊ and var(k) = Ek² − (Ek)² = p[var(k₊) + (Ek₊)²] − (pEk₊)². When moreover Ek = 1 so that Ek₊ = p⁻¹,

3.11

The mutational covariance matrix C is primarily a phenomenological quantity, although in principle it can be expressed in terms of per locus statistics and the genotype to phenotype map,

3.12

As a final point, we need to say something about the traits and the concept of phenotype. In general, the phenotypes of AD should be seen as reaction norms, i.e. maps from micro-environmental conditions to characteristics of individuals (another term is conditional strategies). Only in the simplest cases a reaction norm is degenerate, taking only a single value. The dioecious case is similar in that the development of the sexes need not be, and in fact rarely is, equal. Hence, trait vectors will generally consist of two components, traits of the male, X_m, and of the female, X_f. In general, the traits of the two sexes evolve dependently as they are coupled by their mutational covariances. In the extreme case that the mutational covariances between male and female traits are all zero the female and male coevolve as if they are separate species. At the opposite end of the spectrum X_m = G_m(X_f), or X_f = G_f(X_m), and the mutational variation in the male and female traits is fully correlated. The upshot is that, except in the fully correlated case, we cannot work with a monolithic trait vector influencing both macro- and micro-gametic reproduction, as in hermaphrodites. Instead, we should take into account the fact that when the sexes are separate they can evolve in their own ways.

4. Haplo-diploids: a not uncommon, but often-neglected, reproductive mode

In addition to the haploid and diploid ones, there exist all sorts of other lifestyles. One common type is where haploid and diploid phases alternate (as in, for example, ferns, mosses and a great variety of algae). It is then necessary, as in dioecious diploids, to introduce trait vectors for each separate phase. We can then consider a diploid plus its haploid offspring as a single generalized individual ([36] gives a relaxed introduction to this useful concept) and apply the theory of the previous two sections, where for the diploid phase traits we again have to put in a factor 2.

Still another lifestyle is the so-called haplo-diploid one, where one sex is diploid and the other haploid (the supplementary material to [37] lists the many known haplo-diploid taxa). Although this case can also be treated through the mental construction of appropriate generalized individuals, we follow the strategy of §3, and treat the two sexes as different birth states, as this is simpler and allows us to illustrate further tricks of the trade.

Although the opposite also occurs, we shall for definiteness take the hymenopteran situation as reference and assume that females are diploid. In that case only the female reproductive output needs to be discounted by 1/2. We allow any physiological structure and only assume that within each sex all offspring are born equal. The next-generation matrix then becomes

4.1

For the resident f_f = 1 (as the density of females is constant) and hence f_m = r, the sex ratio at birth (the number of males born relative to females), and m = r⁻¹ (since all females have a father and one female is born for r males). Hence, the resident's next-generation matrix is

4.2

4.3

Lemma.

Let and , then

4.4

So far, we have listed the additional ideas. The rest is work, following the route outlined in §§2 and 3.

We begin with the selection gradient

4.5

4.6

Our arguments for diploids only considered autosomes. To obtain the full CE in the case of chromosomal sex determination, it suffices to add the allosomal contributions to the autosomal one. The contribution of the chromosome occurring only in the heterogametic sex follows the rules for clonal reproduction, that of the chromosome of the homogametic sex follows those for haplo-diploids.

5. Discussion

Away from evolutionarily singular strategies (where ∂₁s(X|X) = 0), the time-scale separation assumption combined with the assumption of small mutational steps guarantees practically permanent genetic homogeneity of the resident population. More in particular, a population with sufficiently restricted variability will become homogeneous on the time scale of the substitutions, thanks to the close to additive genetics, and will effectively stay so until the next substitution. The exceptions to this homogeneity occur when an evolutionary trajectory comes close to an evolutionarily singular strategy. The reason is that, where everywhere else we have an approximately linear selection regime, near the singular strategies the quadratic terms in the local expansion of the fitness landscape start to dominate, creating epistasis (non-additivity) at the level of the fitnesses.

Genetic homogeneity of the resident population lies at the base of the approximations made to reach the CE and of its easy extension to a Mendelian world. At first sight this seems a strong argument against the CE's applicability, as genetic variability appears to be rampant. There are two arguments for yet thinking that the CE may often be a fair description of reality. The first one is that we need mutation limitation only for genes affecting the trait. The thrust of our theoretical deductions is not affected by selectively protected variability that is developmentally and selectively unconnected to the focal traits or variability on neutral loci subjected to mutation, drift and draft (hitchhiking with loci under selection). The only effect that selectively kept variability at unconnected loci may have is that it makes the lifetime offspring numbers of the substituting alleles more variable, enlarging . However, as we treated the components of as empirical quantities, any genetic causes of this variability are automatically included.

The second argument has a bearing on variability owing to a lack of strict time-scale separation. As long as that variability stays sufficiently limited, it should in Mendelian populations have little effect on the quality of the CE approximation. The reason is that the effect of such variability on the invasion fitness is of higher order in the mutational step size than the terms retained in the derivation of the CE. (The argument for the latter statement may be found in [41]. This paper admittedly only considers unstructured populations. However, the nature of the argument suggests that with the right mathematical expertise it should be extendable to structured ones.) In clonal populations, a lack of mutation limitation yet presents an obstacle to the quantitative applicability of the CE owing to clonal interference: a substituting clone may be ousted en route by a better mutant so that only a fraction of the mutants that invade contribute to longer run evolutionary change. We may thus expect ‘reality’ to be slower than the CE. Mendelian populations do not suffer this slowing down as mutants on different loci effectively substitute in parallel without interfering thanks to recombination (and the approximate additivity coming from the smallness of mutational effects).

The picture of multiple substituting mutants with small additive effect may seem close to the one considered by Lande [42,43]. However, there is a difference as in our picture only a small and variable number of loci are simultaneously substituting, with new mutants continually coming and old ones going. By contrast, the quantitative genetics picture underlying Lande's work has its basis in Fisher's picture of evolutionary change coming from changes in the allele frequencies on a large number of loci with all alleles present from the start [44]. (An argument against that picture is that the genetic differences between, say, a choanoflagellate and a human are so many that the attendant polymorphism would contain more genotypes than the number of atoms in the Earth, invalidating the classical population genetic calculations.) Of course, Lande did not subscribe to this simplified picture (e.g. [45,46]), and neither did Fisher, at least not wholeheartedly [47, chapter IV]. However, the mathematical details of the connection between the approximations introduced in Lande's various papers remains to be worked out.

One of the consequences of having only a few loci substituting at any one time is that the within-population variance of the trait fluctuates stochastically and appears to be highly variable (J. Ripa 2005, personal communication). This notwithstanding a negative feedback: when the standing variability is high, mutant alleles will find it more difficult to pass through the stochastic boundary layer because of the larger fluctuations in their offspring numbers. Unfortunately, this picture does not easily lead to explicit expressions, as the fluctuations in variability occur on a time scale similar to that of the substitutions so that we cannot use our earlier branching process calculation. However, we surmise that in reality the allelic reproductive variability coming from co-substituting alleles can usually be neglected relative to the variability in offspring numbers caused by alleles not affecting the traits and by micro-environmental variability such as some individuals running into a predator and others not.

Although the previous paragraphs may give the impression that we are quite optimistic about the applicability of the CE, there is one problem that in applications is usually ignored: there is no need for C to stay constant. It may perhaps seem, when the mutations with small and hence additive effect all occur on different loci, that also their cumulative effect will be additive. However, this is not the case. This can be seen by looking at the developmental process as a sequence of maps, each of which transforms the genetic information plus environmental inputs further towards the final phenotype. Even if linearity were to prevail early on, the accumulated change will be appreciable, and when the output from those early stages is fed into a nonlinear map to get to the next stage, approximate linearity is lost and the derivatives of Φ that appear in (3.12) change over evolutionary time.

In principle, it is possible to write down a CE for the extended trait vector (X, C) to obtain the change in C as correlated response to selection on X. However, going through the calculations in the manner of (3.12) shows that the expression for the covariances between the mutational steps in X and C generically involves higher moments. We should thus see the CE as the first step in a truncated moment expansion.

All this does not mean that using the CE with fixed C is never more than a heuristic tool, with no strong connection to reality. There are scenarios where one may with impunity assume C to be constant, in particular when investigating the behaviour of evolutionary trajectories close to evolutionarily singular strategies. In a linearized stability analysis as in [48] or in the analysis of scenarios for initial divergence close to an adaptive branching point (as defined in [49–51]) as in [4,52–54], only those situations are considered where the change in X, and hence the associated change in C, is small. (Such arguments involve two trait scales, a small one for linearizing the CE and a smaller one for the mutational steps.)

Funding statement

C.G.F.d.K. was supported by the Netherlands Organization of Scientific Research (NWO), grant 813.04.001. J.A.J.M. benefited from the support of the Chair Modélisation Mathématique et Biodiversité VEOLIA-École Polytechnique-MNHN-F.X.c.