## Abstract

Natural selection favors alleles that increase the number of offspring produced by their carriers. But in a world that is inherently uncertain within generations, selection also favors alleles that reduce the variance in the number of offspring produced. If previous studies have established this principle, they have largely ignored fundamental aspects of sexual reproduction and therefore how selection on sex-specific reproductive variance operates. To study the evolution and consequences of sex-specific reproductive variance, we present a population-genetic model of phenotypic evolution in a dioecious population that incorporates previously neglected components of reproductive variance. First, we derive the probability of fixation for mutations that affect male and/or female reproductive phenotypes under sex-specific selection. We find that even in the simplest scenarios, the direction of selection is altered when reproductive variance is taken into account. In particular, previously unaccounted for covariances between the reproductive outputs of different individuals are expected to play a significant role in determining the direction of selection. Then, the probability of fixation is used to develop a stochastic model of joint male and female phenotypic evolution. We find that sex-specific reproductive variance can be responsible for changes in the course of long-term evolution. Finally, the model is applied to an example of parental-care evolution. Overall, our model allows for the evolutionary analysis of social traits in finite and dioecious populations, where interactions can occur within and between sexes under a realistic scenario of reproduction.

- Within-generation variance
- sex-specific evolution
- dioecious population, probability of fixation
- convergent stable strategies

IN the absence of mutation, the change in allele frequency is the result of natural selection and genetic drift. Natural selection favors alleles that maximize their representation within the gene pool, and a large body of work has investigated how alleles achieve this by increasing the expected number of offspring produced by their carriers. However, genetic drift, which arises from randomness in reproduction, reduces the efficacy of natural selection, thereby slowing down or even preventing adaptation altogether.

While many studies have investigated how natural selection affects the expected number of offspring produced by an individual, less attention has been given to the degree to which selection acts on the variance in offspring number, or reproductive variance. Gillespie (1974, 1975, 1977) investigated how natural selection dampens randomness in within-generation fertility in a haploid population. He demonstrated that between two alleles that on average produce the same number of offspring, natural selection favors the allele that produces offspring with lesser variance when the breeding adults of the next generation are sampled from a finite pool of offspring.

Reproductive variance also correlates with the intensity of genetic drift. By decreasing effective population size, reproductive variance mitigates the effect of selection (Wright 1931). Gillespie (1974)’s haploid model also revealed that the level of genetic drift affecting the segregation of two alleles increases with the reproductive variance they code for. As a consequence, fixation of the allele coding for lower fertility variance reduces the intensity of genetic drift and facilitates adaptive evolution in future generations.

The variance in fertility considered in Gillespie’s (1974, 1975, 1977) seminal articles had arbitrary causes and could have resulted from randomness at any stage of an individual’s life history, such as its development, its fertility, or the survival of its offspring. Extensions of Gillespie’s models have since investigated the effect of selection against reproductive variance in the context of more specific life histories. Shpak (2007) investigated the evolution of this variance in age-structured populations and showed that selection favors alleles that code for lower stochasticity in age-specific survival and fertility. Selection against reproductive variance has also been demonstrated to affect the evolution of traits as diverse as sex allocation in hermaphrodites (Proulx 2000), dispersal in structured populations (Shpak 2005; Shpak and Proulx 2007; Lehmann and Balloux 2007), and helping behaviors in social animals (Lehmann and Balloux 2007; Beckerman *et al.* 2011).

These models have highlighted that selection against reproductive variance may be a subtle yet significant force in the evolution of many different traits in natural populations (Rice 2008, for a general discussion). However, it remains unclear how selection on reproductive variance, and its feedback with genetic drift, affect the reproductive biology and life history of sexual organisms. Most models so far have omitted sex altogether. The only study that has taken into account selection on reproductive variance in a dioecious (two-sexes) population approximated reproduction as the random union of gametes and assumed that gamete production of different individuals was uncorrelated (Taylor 2009). These assumptions miss fundamental aspects of the reproductive biology of the vast majority of organisms, where it is individuals, rather than gametes, that unite to mate. But considering a realistic mating system would have significant consequences for the variances and covariances in offspring number among individuals of one sex and how these (co)variances differ across the sexes (Bateman 1948; Wade 1979, for examples). If the effect of the mating system on genetic drift has successfully been captured by calculations of the effective population sizes (Nunney 1993; Nomura 2002), a more general approach is needed to make evolutionary predictions that incorporate selection acting on the traits that generate this variance.

In this article, we present a population-genetic model of male and female phenotypic evolution that makes it possible to predict the evolution of sex-specific reproductive traits under the influence of selection and genetic drift. Using an individual-based approach, the model incorporates a description of the mating system based on the first and second moments (means and (co)variances) of the distribution of individual offspring number. The article is divided into two broad sections. In the first part of the article, we present our model. We start by deriving the probability of fixation of a mutant allele that affects male and/or female reproduction in a finite dioecious population. Our derivation accounts for sex-specific levels of reproductive variance, as well as for covariances between members of the same sex. We then extend the analysis of short-term evolutionary change by deriving predictions of long-term phenotypic evolution, which in turn makes it possible to calculate the equilibrium for sex-specific phenotypes based on the probability of fixation and properties of the mutational input. In the second part of the article, we provide an illustration of how our model can be applied and study the effects of sex-specific reproductive variance on the evolution of various parental-care strategies. These simple examples allow us to demonstrate that sex-specific reproductive variance can lead to differences between the probability of fixation for mutants affecting female traits and those affecting male traits and between the phenotypic equilibria of these traits.

## Model and Results

### Probability of fixation in dioecious populations

We derive the probability of fixation of a mutant allele (*A*) introduced as a single copy into a population fixed for a resident allele (*a*). The population is dioecious and of any, but constant, size, with *N*_{m} adult males and *N*_{f} females. Generations are nonoverlapping. The mutant *A* allele alters the expression of a continuously varying phenotypic trait, which may affect one or more aspects of reproductive biology, such as mating success, fecundity, or offspring survival. The trait can have different values in males and females and we denote by *z*_{m} the phenotypic value of a male homozygous for the resident allele (genotype *aa*). The phenotype of a heterozygous male (*Aa*) is denoted by *z*_{m} + *hδ*_{m}, where *h* is the dominance coefficient of *A*. A male homozygous for the mutant allele (*AA*) has phenotype *z*_{m} + *δ*_{m}. Similarly, the phenotypes of the three genotypes in a female are *z*_{f} (*aa*), *z*_{f} + *hδ*_{f} (*Aa*), and *z*_{f} + *δ*_{f} (*AA*).

#### Weak selection approximation:

The fixation probability of the mutant is derived using an individual-based approach that builds on previous works (Rousset 2003; Roze and Rousset 2004; Lessard and Ladret 2007; Lehmann and Rousset 2009; and supporting information, File S1). Under weak selection (small mutant deviations *δ*_{m} and *δ*_{f}), and if the mutation rate is the same in males and females, the fixation probability of a single mutant copy arising at random in a monomorphic population with phenotypes *z*_{m} in males and *z*_{f} in females is (1)where *N* = *N*_{m} + *N*_{f} is the total number of adults, and *δ* is such that *δ*_{m} ∼ *O*(*δ*), *δ*_{f} ∼ *O*(*δ*) (File S1, Equations SI.1–SI.39). The functions *G*_{m}(*z*_{m}, *z*_{f}) and *G*_{f}(*z*_{m}, *z*_{f}) are fitness gradients: they measure the effect of the mutant on male and female fitness and are further explained below (*Fitness gradients* section). The functions *K*_{m}(*z*_{m}, *z*_{f}, *h*) and *K*_{f}(*z*_{m}, *z*_{f}, *h*) are measures of the variance of mutant frequencies in males and females over the segregation process, from the appearance until the eventual fixation of the mutant. They are inversely proportional to the intensity of genetic drift and capture the efficacy with which selection acts on the mutant (see *Efficacy of selection* section). Whether a mutant is under positive selection [*π* > 1/(2*N*)], evolves neutrally [*π* = 1/(2*N*)], or is counterselected (*π* < 1/(2*N*)), therefore depends on the balance between male and female mutant effects *δ*_{m} and *δ*_{f}, the fitness gradients *G*_{m} and *G*_{f}, and the measures of genetic variances *K*_{m} and *K*_{f}. In the following sections, we lay out how these quantities depend on the reproductive system of a population.

#### Fitness gradients:

The fitness gradients express how the expected number of adult offspring of a focal individual changes in response to small alterations of its phenotype. They are given by (2a) (2b)where is the expected number of adult offspring of sex *u* of a focal individual of sex *v* ∈ {m, f}. This fitness function depends on the phenotype *z _{vi}* of the focal individual of sex

*v*., on the average phenotype among sex

*v*, but excluding the focal, and on the average phenotype in the population of the opposite sex ( for males and for females). The model can thus easily accommodate for sex-specific interactions based on games, like the classic battle of the sexes. The derivatives of focal fitness are evaluated at the resident phenotypes , , so that

*G*

_{m}and

*G*

_{f}measure the effects of phenotypic changes on male and female fitness with respect to the resident population.

The fitness gradients in Equation 2a indicate the direction of phenotypic evolution in each sex that is favored by selection. If *G*_{m}(*z*_{m}, *z*_{f}) and *G*_{f}(*z*_{m}, *z*_{f}) are positive, then selection favors an increase in the trait in males and females; if they are negative, then selection favors a decrease. Although the gradients are derivatives of the expected average number of adult offspring produced ( and ), the direction of selection depends on how the phenotype affects the average number of juvenile offspring produced as well as reproductive variance. To demonstrate why this is the case, we derive the fitness of a focal individual *i* of sex *v* in terms of the distribution of juveniles in the population. Fitness is the expected number of *i*’s offspring that become part of the adult breeding population of the next generation. We separate fitness gained through male and female offspring. We write for the number of juvenile offspring of sex *u* born to *i*, itself of sex *v*. In each generation, the set of reproductive individuals is established by independently sampling *N*_{m} males from a pool of surviving male offspring and *N*_{f} females from a pool of surviving female offspring. The conditional fitness of individual *i* of sex *v* gained through offspring of sex *u*, *i.e.*, the expected number of its adult offspring of sex *u*, can then be calculated in terms of *i*’s reproductive success, relative to that of the total population (3)where is the realized offspring production of all parents in the population. Note that the total number of juveniles of either sex must be the same when counted as the offspring of males or females (*i.e.*, ). We assume nonextinction of the population . To describe the fitness of individual *i*, we need to calculate the expectation of Equation 3 over the distribution of . Following the approach of previous work (Gillespie 1975; Proulx 2000; Shpak and Proulx 2007; Lehmann and Balloux 2007; Rice 2008), is approximated using the delta method (Oehlert 1992), so that the expected fitness becomes (4)where is the expected number of juveniles of sex *u* produced by individual *i*, is the expected total number of juveniles of sex *u* produced in the population, is the variance of the number of offspring of individual *i* , and is the covariance between the number of offspring of sex *u* of individuals *i* and *k* of sex *v* . The fitness of an individual therefore takes into account all first and second moments of the probability distribution that describes individual reproduction (see Figure 1 for a depiction of those moments for a focal male).

The remainder *R* in Equation 4 is composed of central cross moments of of order 3 and higher. These terms may be significant in certain scenarios (Rice 2008), but we omit them in this analysis by assuming that the distribution of is well behaved as the population size *N* increases. Previous models used the central limit theorem to justify this assumption (Shpak and Proulx 2007; Lehmann and Balloux 2007, Equation A6). This is not strictly valid here because the number of offspring produced by different individuals is not necessarily independent. However, the remainder terms can be ignored if we assume that offspring numbers are close to independence and that the “total” covariance between a given set of individuals decreases as the number of individuals in that set increases (see *Appendix*, Equation A1). In this case, the remainder terms in Equation 4 are of order (1/*N*^{2}). Therefore, while the expression for the fixation probability (Equation 1) holds for any population size, the approximation for fitness (Equation 4) takes into account only the first-order effects of finite population size on offspring means and variances. If condition (A1) also holds for the first and second moments of the distribution of reproduction, then the effects of (co)variances on individual fitness vanish as *N* → ∞, in agreement with previous studies (Gillespie 1974).

Equation 4 shows that individual fitness depends on four terms. The first is the relative expected number of offspring produced , which has a positive effect on fitness. The remaining three terms capture the effects of reproductive variance. Fitness decreases with the variance in offspring number , increases with the variance in offspring number produced by the remaining individuals in the population ( and ), and decreases with the covariance between the number of offspring produced and that of the remaining individuals in the population .

The fitness effects of the variance terms stem from the nonlinear relationship between fitness and the offspring production of the focal (; see Figure 2 and Equation 3) and the rest of the population (; see Figure 2B and Equation 3). For a given number of offspring produced by the rest of the population, the fitness returns of a focal individual diminishes with its production of more offspring as a consequence of the increased competition between related juveniles for access to breeding. This results in a net negative effect of variance in the focal individual’s reproductive output on its fitness ( in Equation 4 and Figure 2A). Conversely, for a given offspring production by the focal, the advantage of competing within a population of individuals that are on average less fecund is expected to be greater than the disadvantage of competing in a more productive population. This leads to a net positive effect of population variance on the focal individual’s fitness ( and in Equation 4 and Figure 2B). Finally, using arguments similar to those presented in Figure 2, one can see that the benefit of overperforming in a less competitive population is on average greater than the cost of underperforming in a more competitive population. As a consequence, the covariance between the numbers of offspring produced by the focal individual and the rest of the population has a negative impact on focal fitness ( in Equation 4).

Selection on a phenotype (Equation 2) then reflects the balance between the impact of the trait on the different terms of Equation 4. Since the difference between two phenotypes are small [of order *O*(*δ*)], we can describe the dependence between the moments and phenotypes without explicitly characterizing the interactions between every individual in the population. Rather, we average the sums in Equation 2 over mean population phenotypes (File S1, Equation SI.4). Then, if the trait of interest affects all first and second moments of individual reproduction, the fitness function of a focal individual of sex *v* can be written as (5)where denotes the average phenotype of the sex opposite to that of the focal, *v*, (*e.g.*, ). The function is the expected number of juveniles of sex *u* produced by a focal of sex *v*, and is the average expected number of juveniles of sex *u* produced by sex *v* individuals in the population; similar interpretations are given to the variance and covariance functions ( and ). Therefore, calculating the individual fitness functions that go into the fitness gradients (Equation 2) requires characterizing only the individual mean, variance, and covariance functions, (, , and ), and these depend only on the phenotype of the focal individual (*z _{vi}*), the average phenotype in the opposite sex , and the average among other individuals of the same sex (, as the average is written as ). Examples of such calculations are given in the

*Example*section.

The first line in Equation 5 reflects the fact that an individual who produces on average a greater number of offspring than the average individual in the population has higher fitness. The second and third lines reflect the fact that an individual with a lower variance in progeny number than the average individual has higher fitness, as originally described by Gillespie (1974). Finally, the last two lines of Equation 5 reflect the fact that an individual whose offspring production covaries with that of another individual to a lesser degree than the average individual also has higher fitness. In addition, we see that the effect of covariance on fitness [of order (*N _{v}* − 1)/

*N*] is potentially greater than that of the variance (of order 1/

_{v}*N*).

_{v}#### Efficacy of selection:

In addition to the fitness gradients, the fixation probability (Equation 1) also depends on *K*_{m}(*z*_{m}, *z*_{f}, *h*) and *K*_{m}(*z*_{f}, *z*_{f}, *h*), which weigh on the fitness gradients and measure the sex-specific efficacy of selection in males and females, respectively. The weight *K*_{m} (*K*_{f}) is the expected covariance between between genotype and phenotype in males (females), cumulated over the neutral segregation of the mutant. Mathematically, this is (6)where *C*[*z _{ui}*,

*p*]

_{ui}*is the covariance between the phenotype*

_{t}*z*of an individual of sex

_{ui}*u*and the frequency

*p*∈ {0, 1/2, 1} of the mutant it carries at generation

_{ui}*t*, and

*E*[⋅] denotes expectation under neutral evolution,

*i.e.*, where only genetic drift affects fluctuations of genotypic frequency

*p*from one generation

_{ui}*t*to the next (File S1, Equations SI.20, SI.21, and SI.39). The factor 1/(2

*N*) in Equation 6 is the initial frequency of the mutant, while the factor 1/4 is the product of the frequency of transmission of a gene by a parent to an offspring (

*i.e.*, 1/2; File S1, Equation SI.2) and the reproductive value of that class of offspring, which is here 1/2 for both males and females. If mutant effects are additive (

*h*= 1/2), then

*C*[

*z*,

_{ui}*p*]

_{ui}*=*

_{t}*V*[

*p*]

_{ui}*and*

_{t}*K*reduces to the cumulative genetic variance in sex

_{u}*u*, highlighting that the larger genetic variance is, the more efficient selection can be. In the simple case of asexual haploids, this variance reduces to the familiar

*p*(1 −

_{t}*p*), where

_{t}*p*is the average frequency of the mutant at generation

_{t}*t*. More generally,

*K*

_{m}(

*K*

_{f}) depends on dominance and captures the association between phenotypic and genetic variance in males (females) on which selection is then able to act.

Since the calculations for *K*_{m} and *K*_{f} are made over the segregation of a neutral mutant (Equation 6), they can be expressed in terms of coalescence times of neutral genes (File S1, Equations SI.37–SI.38), which themselves can be expressed in terms of how genes coalesce within individuals of different sexes. Doing so links *K*_{m} and *K*_{f} back to the mating system and hence to the evolving reproductive traits *z*_{m} and *z*_{f} that are under study. The general expressions for *K*_{m} and *K*_{f} in terms of the coalescence process depend on the probabilities that individuals share the same parent in the absence of selection, referred to here as “probabilities of sibship.” We find that the coalescence process can be described by 14 probabilities of sibship. Six of these describe the probability that a pair of individuals share the same parent; they are written as , where *u* ∈ {♂, ♀} indicates the sex of the parent, and *v* ∈ {m, f, *c*} indicates whether a pair of individuals consist of two males, two females or a male and a female (Figure 3). The remaining eight probabilities of sibship describe the probability of three individuals (three males, two males and a female, two females and a male, or three females) sharing the same parent (male or female). Providing a general characterization of the neutral coalescence process is complicated, but the system can be simplified by taking into account only the first-order effects of finite population size [*O*(1/*N*)]. In this case, we find that the eight three-way probabilities of sibship may be expressed in terms of the pairwise probabilities ’s (File S1, Equations SI.40–SI.42), which are then sufficient to describe the entire coalescence process.

The pairwise probabilities of sibship capture different aspects of reproductive variance. The probabilities that two males have the same father , that two females have the same father , and that a male and female have the same father measure the level of reproductive variance of adult males and depend on the moments of the distribution of male reproduction (Table 1). In a situation where all males contribute equally to the next generation of individuals of either sex, the paternal probabilities of sibship are all (*x* ∈ {m, f, *c*}). With skewed paternity, the variance in the number of offspring a male has increases and so do the paternal probabilities of sibship (Table 1). The paternal probabilities of sibship , , and differ from one another if the variance in the number of sons produced differs from the variance in the number of daughters produced, *i.e.*, if gene transmission is more variable through offspring of one sex than through offspring of the other sex.

The difference between maternal and paternal probabilities of sibship reflect the difference between the reproductive variances of males and females. In a monogamous population, with equal number of males and females, males and females have the same reproductive variance and . In contrast, a polygynous population has greater reproductive variance in males than females , while a polyandrous population exhibits the opposite pattern .

Using the probabilities of sibship, we can express the weights *K _{u}* in terms of the coalescence process of a neutral gene. In the simplest case, where a mutant is additive and there is no difference between the probabilities of a gene being transmitted through a son or a daughter ( and ),

*K*(Equation A2) simplifies to (7)This expression decreases hyperbolically with the probabilities of sibship. Thus, as reproductive variance increases in males and females, the transmission of the gene from one generation to the next becomes more stochastic and genetic variance is reduced. As a consequence, the efficacy of selection

_{u}*K*decreases, reflecting the smaller impact of selection on the probability of fixation in the face of increased drift (Equation 1).

_{u}Similar patterns are observed when reproductive variance varies with the sex of the parent and the sex of the offspring. Analytical results for additive mutants (Equation A2, Figure 4, A–C) and numerical results for nonadditive mutants show that *K*_{m} and *K*_{f} both decrease hyperbolically with all six probabilities of sibship. Numerical results also show that *K*_{m} and *K*_{f} increase linearly with dominance (Figure 4D). This stems from the fact that dominance increases the phenotype–genotype covariance at lower allele frequencies (*p _{ui}* < 1/2) and that the frequency of neutral mutants remains on average low.

Because *K*_{m} and *K*_{f} weight male (*G*_{m}) and female (*G*_{f}) fitness gradients independently in the probability of fixation (Equation 1), reproductive variance may affect male and female evolution differently. For example, selection on females has a greater impact on the probability of fixation than selection on males when *K*_{m} < *K*_{f}. In this case, a female-limited mutation (*δ*_{m} = 0 or *G*_{m} = 0) would have a greater chance of reaching fixation than a male-limited mutation (*δ*_{f} = 0 or *G*_{f} = 0), even if both improve fitness by the same amount. In the longer term, we would then observe a faster rate of adaptation in females than males. The reverse patterns are predicted when *K*_{m} > *K*_{f}.

Differences between *K*_{m} and *K*_{f}, and hence differences between the efficacy of selection in males and females, occur whenever genetic variance is lower in one sex than in the other. For additive mutants (Equation A2 and Figure 4, A and B), *K*_{m} < *K*_{f} requires that (8)*i.e.*, the probability that two males have at least one parent in common exceeds the probability that two females have at least one parent in common. This inequality (Equation 8) reflects that if at each generation male offspring are more related than female offspring, then genetic variance in males is lower than in females, and as a consequence, selection is less efficient in males.

#### Calculating the probability of fixation:

Based on the above derivations, the probability of fixation can be explicitly calculated, taking into account the fitness change caused by the mutant, through its effect on first and second moments of the distribution of offspring production and the impact of segregation in the two sexes on the efficacy of selection. To calculate the fixation probability, the probabilities of sibship (Table 1) are substituted into *K*_{m} and *K*_{f} (Equation A2 if the mutant is additive and File S1, SI.38 otherwise), and the expressions for focal fitness (Equation 4) are substituted into the fitness gradients *G*_{m} and *G*_{f} (Equation 2). Finally, *K*_{m}, *K*_{f}, *G*_{m}, and *G*_{f} are substituted into the fixation probability (Equation 1).

### Long-term phenotypic evolution in dioecious populations

The fixation probability of a mutant is useful for predicting short-term evolution and to understanding how the interplay between selection and genetic drift affects the fate of a new mutation. However, it is often desirable to predict the long-term evolution of phenotypes as a result of selection and drift acting on an influx of new mutations. In this section, we use the probability of fixation (Equation 1) to determine the phenotypes most likely to be observed in males and females at a selection–mutation–drift balance. To that aim, we assume that the autosomal locus mutates at a constant rate *ξ*. This rate is sufficiently weak compared to the rate of fixation to ensure that there are only ever two alleles segregating, thus complying with the weak-mutation strong-selection limit of population genetics (e.g., Gillespie 1994; Sella and Hirsh 2005) and/or the trait substitution sequence limit of evolutionary game and inclusive fitness theory (e.g., Metz *et al.* 1995; Champagnat and Lambert 2007; Lehmann 2012). The effects (*δ*_{m}, *δ*_{f}, *h*) of a mutation are drawn from a distribution *u*(*δ*_{m}, *δ*_{f}, *h*), which is such that the dominance (*h*) of a mutant is independent from its homozygotic effects (*δ*_{m}, *δ*_{f}), and mutants have on average no phenotypic effect *E*[*δ*_{m}] = *E*[*δ*_{f}] = 0. The rate at which a population monomorphic for (*z*_{m}, *z*_{f}) is substituted by a population with traits (*z*_{m} + *δ*_{m}, *z*_{f} + *δ*_{f}) can then be written as (9)where 2*N* is the number of gene copies in the population, ξ*u*(*δ*_{m}, *δ*_{f}, *h*) is the probability that a single copy produces a mutation of type (*δ*_{m}, *δ*_{f}, *h*), and the term within brackets is the probability that this mutation fixes in the population, which is given by Equation 1.

The substitution rate *k* determines a jump process (Gardiner 2009), which describes the stochastic evolution of male and female traits as jumps between monomorphic states in phenotypic space. Ignoring terms of order *O*(*δ*^{2}) in Equation 9, the jump process can be described in continuous time by a diffusion process that eventually reaches a stationary distribution *ψ*(*z*_{m},*z*_{f}) (*Appendix*, *Diffusion equation for phenotypic evolution in dioecious populations*). This long-run stationary state reflects a balance between the forces of mutation, selection, and genetic drift, and the maxima of *ψ*(*z*_{m}, *z*_{f}) correspond to phenotypes around which the populations spend the greatest amount of time. These maxima are the most likely outcomes of phenotypic evolution, and in single phenotype models, they are the “convergence stable” states of the system (Lehmann 2012). A phenotype is convergent stable if populations sitting close to this phenotype are attracted toward it. Convergence stability is an important concept of attainability of equilibrium points, common to evolutionary game and inclusive fitness theory (Rousset and Billiard 2000; Leimar 2009). Under the trait-substitution sequence limit, a phenotype that is convergent stable is also evolutionary stable (Wakano and Lehmann 2012).

Phenotypes that are multidimensional convergence stable can be found by considering the attractor points of the system of differential equations(10)which describes the deterministic trajectory associated to the underlying diffusion process. Here, is the average dominance of the mutation distribution and *b _{uv}* =

*C*[

*δ*

_{u},

*δ*

_{v}] is the covariance between mutation effect in sex

*u*and

*v*. If none of the attractor points of system Equation 10 lie on the boundary of the phenotypic space, then large deviation theory shows that as the population size grows, the stationary distribution

*ψ*(

*z*

_{m},

*z*

_{f}) becomes peaked around these attractor points (use Theorem 4.3 of Freidlin and Wentzell 2012, p. 170, and observe that if none of the attractor points lies on the boundary of the phenotypic space, then a smooth domain can be drawn around all attractor points, thereby satisfying condition A of Freidlin and Wentzell 2012, p. 150). Therefore, when all attractor points of Equation 10 lie in the interior of the phenotypic space, they correspond to the convergence stable states and are the most likely phenotypic outcomes of evolution. Furthermore, in the infinite size limit (

*N*→ ∞), the stationary distribution becomes fully concentrated around a single of these convergence stable states (Theorem 4.2 of Freidlin and Wentzell 2012, p. 167), which corresponds to the highest peak of the adaptive landscape and the stochastic stable state of the system (

*e.g.*, Foster and Young 1990; Van Cleve and Lehmann 2013).

For an interior point (*z*_{m}, *z*_{f}) of system Equation 10 to be convergence stable, two conditions must be satisfied. First, *dz*_{m}/*dτ* = 0 and *dz*_{f}/d*τ* = 0 must hold, but since and , this condition is equivalent to (11)*i.e.*, that the male and female fitness gradients vanish, which is equivalent to the condition for establishing singular points in deterministic evolution (Leimar 2009). Second, the real part of all the eigenvalues of the Jacobean matrix of Equation 10 must be negative. Combined with Equation 11 and the properties of *b _{uv}* (Leimar 2009), this latter condition is equivalent to the real part of the eigenvalues of (12)being negative at . Conditions (11) and (12) extend the one-dimensional condition of convergence stability for finite populations (Rousset and Billiard 2000; Lehmann 2012), and when

*K*

_{m}=

*K*

_{f}, it is equivalent to the condition for multidimensional convergence stability for populations of infinite size (Leimar 2009), which depends only on the fitness gradients

*G*

_{m}and

*G*

_{f}. When attractor points of Equation 10 lie on the boundary of the phenotypic space or outside of it, the shape of the equilibrium distribution cannot in general be assessed (to the best of our knowledge), and in this case, characterizing convergence stable points is not straightforward.

## Example

In this section, we illustrate a possible application of our model by analyzing the evolution of maternal and paternal care behaviors. The emphasis is on highlighting the effects of selection on reproductive variance in driving the evolution of reproductive traits. We consider a dioecious species with a simple life cycle. An equal number *N* of adult males and females randomly pair up to mate. All females mate once with a single male, but males mate with harems of exactly *H* females. If *H* = 1, the population is monogamous and all males mate. If *H* > 1, then the population is polygynous and some males mate *H* times while others do not reproduce at all. Each female gives birth to exactly *f* offspring with an equal sex ratio. The offspring survive with probabilities that depend on the phenotypes *z*_{m} and *z*_{f} of their parents. Female offspring each survive independently from each other with probability *s*^{f}(*z*_{m}, *z*_{f}). In contrast, the survival of males is strongly correlated within broods and the males of a brood either all survive [with probability *s*^{m}(*z*_{m}, *z*_{f})] or all die [with probability 1 − *s*^{m}(*z*_{m}, *z*_{f})]. This is close to the assumptions of the hermaphroditic model of Proulx (2000). The difference in male and female survival could reflect sex differences in the susceptibility to random variation in the breeding environment provided by the male’s territory, for example. As a consequence of their correlated survival, surviving males are more related than surviving females. The next generation of adults is randomly sampled from the population of surviving offspring as described in the previous section.

The phenotypes that evolve are the level *z*_{m} of paternal care provided by a male and the level *z*_{f} of maternal care provided by a female. The survival probabilities of daughters, *s*^{f}(*z*_{m}, *z*_{f}), and sons, *s*^{m}(*z*_{m}, *z*_{f}), are both increasing functions of *z*_{m} and *z*_{f}. We analyze the fate of four different types of mutations altering the parental-care phenotypes. These types are characterized by the sex of the parent providing the care and the sex of the offspring receiving it. We distinguish between mutations that affect care for sons only and those that affect care for daughters only. For both of those, we consider sex-limited mutations that affect care only in males or only in females. In total, the evolution of four different traits is studied, paternal care for daughters, maternal care for daughters, paternal care for sons, and maternal care for sons.

The covariance between the survival of male offspring depends on the sex of the parent for which care evolves. If female care is evolving, then for each female, her entire male brood either survives or dies, independently from other females, even from those that have mated with the same male. When care is provided by males, then for each male, his entire male brood either survives or dies, independently from other males, but his brood includes all the male offspring he has had with different females. We assume that the effect of mutations is additive and identical for all traits. Thus, compared to the resident homozygote, the level of care is increased by an amount *δ*/2 in heterozygotes and by *δ* in mutant homozygotes.

### Probability of fixation of new mutants

To calculate the fitness of a focal individual (Equation 5) and evaluate the fixation probability of the different mutants, we need to express the means and covariances of offspring production in terms of care phenotypes. We first consider the offspring production of a focal female *j* with phenotype *z*_{f}* _{j}*. Because the evolution of male and female traits are treated separately, mutations for altered maternal care always occur in the presence of constant resident paternal care

*z*

_{m}, and we can therefore omit paternal care from the survival functions. The expected number of daughters and sons the focal female produces are given by (13a) (13b)where the factor 1/2 stems from the equal sex ratio.

The variance terms that contribute to the fitness of a focal female are found by writing *n*_{♀} and *n*_{♂} as the number of daughters and sons, respectively, at birth before survival selection. With the sex ratio being equal, *n*_{♀} ∼ Bin(*f*, 1/2) (and *n*_{♂} = *f* − *n*_{♀}). Given *n*_{♀}, the variance in the number of female offspring born to a female is *n*_{♀}*s*^{f} (*z*_{f}* _{j}*)(1 −

*s*

^{f}(

*z*

_{f}

*)), and the mean is*

_{j}*n*

_{♀}

*s*

^{f}(

*z*

_{f}

*). Therefore, applying the law of total variance, the variance in the number of female offspring born to a focal female is (14)In contrast, since sons do not survive independently from each other, we have, given*

_{j}*n*

_{♂}, that the variance in the number of male offspring is and that the mean is

*n*

_{♂}

*s*

^{m}(

*z*

_{f}

*). Thus, the variance in the number of sons produced by a focal female differs from the variance in the number of daughters and is (15)Finally, since females give birth to and care for their offspring independently of one another, the covariance between the number of offspring of two females is zero .*

_{j}The means and variances of the offspring numbers produced by a focal male *i* are obtained similarly (see *Appendix*, *Calculations for the evolution of parental care*), but in contrast to singly mating females, polygyny leads to a negative covariance between the numbers of offspring sired by two males. The additional variance terms ( for male offspring and for female offspring) must be taken into account when determining the fitness of a focal male and are calculated here. Since maternal care is now constant, we can omit the female care phenotype *z*_{f} from the survival functions. By definition, the covariance between the number of male offspring fathered by the focal male *i* and an “average” male other than *i* is (16)where *E*[*N*_{mat}_{i}N_{mat−}* _{i}*] is the expected product between the number of matings of the focal male

*i*,

*N*

_{mat}

*, and that of another male,*

_{i}*N*

_{mat−}

*. Then, since*

_{i}*N*/

*H*males are chosen at random without replacement to mate exactly

*H*times,

*E*[

*N*

_{mat}

_{i}N_{mat−}

*] = (*

_{i}*N*−

*H*)/(

*N*− 1), and (17)Similarly, the covariance number of females fathered by the focal male

*i*and an average male is (18)As expected, these covariances vanish in a monogamous population (

*H*= 1), but become increasingly negative as fewer males mate. For large

*H*, they contribute significantly to the fitness of a focal male.

#### Fitness gradients:

Substituting Equations 13a and 14 into Equation 5 and deriving according to Equation 2b give the fitness gradient for alleles that code for maternal care of daughters. Similarly, substituting Equations 13b and 15 into Equation 5 and differentiating according to Equation 2b give the fitness gradient for alleles that code for maternal care of sons. The fitness gradients for alleles that code for paternal care are obtained in the same way using Equation 2a.

To identify the different effects of sex-specific reproductive variance, we first consider a population that is strictly monogamous (*H* = 1). If maternal and paternal care have the same effect on offspring survival [*i.e.*, ∂*s ^{v}*(

*z*

_{m},

*z*

_{f})/∂

*z*

_{f}= ∂

*s*(

^{v}*z*

_{m},

*z*

_{f})/∂

*z*

_{m}], the fitness gradients for mutants that increase maternal (

*G*

_{f}) or paternal (

*G*

_{m}) care of daughters are identical and equal to (19)where

*u*∈ {m, f} and

*z*

_{m}and

*z*

_{f}are the resident levels of paternal and maternal care. Likewise, there is no difference between the fitness gradients for mutations that increase maternal (

*G*

_{f}) or paternal (

*G*

_{m}) care of sons, which are both (20)

*u*∈ {m, f}. The (1 − 1/

*N*) term of Equations 19 and 20 describe the balance between the advantage of increasing the expected number of offspring of the focal individual and simultaneously increasing the total expected number of offspring in the population and the level of competition between kin. This is equivalent to the first term of Equation 5. It is equal for gradients describing care for sons and daughters, in line with the fact that the effect of care on the expected numbers of male and female offspring is identical. The second term in the brackets of Equations 19 and 20 captures the increase in fitness due to the reduction in the variance of offspring number (reflecting the remaining terms of Equation 5).

The variance term is greater for mutants that alter the care of sons (Equation 20) because the variance in the number of surviving sons is inherently greater than that of surviving daughters. As a consequence, reducing that variance generates proportionately greater fitness benefits. The benefit of reducing the variance in male offspring number is so large that it completely offsets the reduction in fitness due to kin competition (Equation 20). The benefits of decreasing the variance in the number of surviving daughters vanish as brood size *f* becomes large (Equation 19), due to the fact that daughters survive independently of each other. As a result of these different effects, selection favors the fixation of mutants that increase the care of sons (Equation 20) with greater strength than those that increase the care of daughters (Equation 19), especially when fertility is high.

#### Efficacy of selection:

Differences in the patterns of male and female survival also affect the efficacy of selection on paternal and maternal care, *K*_{m} and *K*_{f}. Both coefficients are calculated using the probabilities of sibship (Equation A2), which are themselves expressed in terms of the moments of offspring production in Table 1. These moments are the same as those appearing in the calculation of focal fitness and in addition include the covariance between the number of male and female offspring produced by a resident individual (*Appendix*, *Calculations for the evolution of parental care*, for calculations). We find that *K*_{m} and *K*_{f} both increase with male and female survival but that their difference (21)[with *f* ∼ *O*(*N*)] depends only male survival rate *s*^{m}(*z*_{m}, *z*_{f}). In the extreme case where all males survive (*s*^{m}(*z*_{m}, *z*_{f}) = 1), selection is as efficient for male and female traits (*K*_{m} = *K*_{f}). But as male survival rate decreases, the efficacy of selection falls more rapidly in males than females. This is caused by the different modes of survival for male and female offspring. Because male offspring tend to be more related than female offspring, genetic variation in males is lower. As a consequence, the efficacy of selection in females is greater than that in males (*K*_{f} > *K*_{m}), and alleles that code for maternal care are under more efficient selection than those for paternal care.

#### Probability of fixation:

Combining the weights *K*_{m} and *K*_{f} (Equation 21) with the fitness gradients *G*_{m} and *G*_{f} (Equations 19 and 20) for the probability of fixation (Equation 1), we find that differences between male and female survival affect the evolution of paternal and maternal care in two ways. First, the correlation in male survival creates a stronger selection pressure on increased care for sons than for daughters. Second, the effect of more stochastic male survival increases drift in males and makes selection on paternal care less efficient than selection on maternal care. As a consequence of these effects, the most probable form of parental care to evolve in our model is maternal care for sons and the least probable is paternal care for daughters (Figure 5 for *H* = 1). Obviously, these conclusions are conditional on the assumptions underlying our analyses, most importantly that mutations affect only male and female care, and do so equally, and that increases in paternal and maternal care result in identical changes in the expected survival of sons and daughters.

#### Mating system:

Polygyny generates a negative correlation between the reproductive output of different males, and this affects the evolution of parental care in two ways. First, reproductive skew in males decreases the strength of selection on male care, due to intensified sib competition among the surviving offspring of a male. This effect can be seen when inspecting the fitness gradients on paternal care for a female offspring (22)and that on a male offspring (23)where both equations are here shown for a population that is strongly polygynous [*H* ∼ *O*(*N*)]. Equations 22 and 23 correspond to the gradients in a monogamous population (Equations 19 and 20), with the exception of the last negative term. This term expresses the cost of intensified sib competition and increases with the level of polyandry *H*. The gradients for the care of female offspring are unaffected by polygyny (see Equations 19 and Equation 20).

Second, polygyny affects the efficacies of selection *K*_{m} and *K*_{f} (*Appendix*, *Calculations for the evolution of parental care*, for calculations). Polygyny, and the associated increase in male reproductive variance, generate additional genetic drift and reduce both *K*_{m} and *K*_{f} relative to monogamy. However, when male brood are brothers through their father, the genetic variance in male offspring decreases with harem size. As a consequence, the depreciation in *K*_{m} with *H* is steeper than in *K*_{f} (24)indicating that the evolution of paternal care is more sensitive to polygyny than the evolution of maternal care. The joint effects of reduced selection and lower efficacy of selection on males compromise the evolution of parental care in polygynous populations. Figure 5 shows analytical predictions of the probability of fixation for varying levels of polygyny. These show that mutants for parental care become less likely to fix as the level of polygyny increases, and this effect is exacerbated for paternal care, demonstrating the double effect of reproductive variance on both the intensity and the efficacy of selection on reproductive traits.

### Long-term evolutionary equilibrium

The probability of fixation captures evolutionary dynamics over a short timescale. But as shown above (Equations 9–12), it can be used to predict long-term phenotypic outcome under a recurrent inflow of mutations. We explore here how reproductive variance affects the long-term evolution of parental care.

In the above parental care example, the evolutionary dynamics governed by the selection gradients (Equations 19, 20, 22, and 23) eventually reach the trivial equilibrium phenotypes of maximum care for sons and daughters [*s*^{m}(*z*_{m}, *z*_{f}) = *s*^{f}(*z*_{m}, *z*_{f}) = 1]. We therefore introduce the realistic assumption that a parent’s resources are limited and that as a consequence, there is a trade-off between the efforts allocated to sons and daughters. The care provided to male offspring by a parent of sex *u* is written *z _{u}* and 1 −

*z*is the care allocated to daughters (with 0 ≤

_{u}*z*≤ 1). As, before, the survival

_{u}*s*(

^{v}*z*

_{m},

*z*

_{f}) of an offspring of sex

*v*is a function of the paternal and maternal care received. Here, a simple additive function is assumed, where the survival of a male offspring is

*s*

^{m}(

*z*

_{m},

*z*

_{f}) = (

*z*

_{m}+

*z*

_{f})/2, while that of a female offspring is

*s*

^{f}(

*z*

_{m},

*z*

_{f}) = (1 −

*z*

_{m}+ 1 −

*z*

_{f})/2. Then, because we consider the evolution of care in one sex while maintaining the care of the other sex constant, Equation 10 shows that the long-term evolution of sex-limited traits can be inferred from selection on that sex alone. In other words, we can predict the phenotypic equilibrium of a male-limited trait (

*b*

_{mf}= 0 and

*b*

_{ff}= 0) from the zeroes of the fitness gradient

*G*

_{m}(

*z*

_{m},

*z*

_{f}) = 0 of males and that of a female-limited trait (

*b*

_{mf}= 0 and

*b*

_{mm}= 0) from the zeroes of

*G*

_{f}(

*z*

_{m},

*z*

_{f}) = 0.

Then, convergent stable states are found in three steps. First, using calculations similar to those used in Equations 13–18 and *Appendix*, *Calculations for the evolution of parental care*, we calculate the moments of reproduction for a focal male and a focal female in the presence of trade-offs to find the fitness components of a focal individual (, , Equation 4). Second, we add the fitness components of a parent of sex *u* derived from male and female offspring to obtain the gradient *G _{u}*(

*z*

_{m},

*z*

_{f}) as in Equation 2. Finally, solving for

*G*(

_{u}*z*

_{m},

*z*

_{f}) = 0, we find the convergence-stable level of investment in sons, is identical for male and female parents when the population is monogamous, and such that male survival is (25)This equation show that the equilibrium investment approaches 1/2 as population size goes to infinity (

*N*→ ∞). This prediction is in line with the fact that kin competition vanishes in infinite populations and with it the selection pressures emanating from reproductive variance. Parents in very large populations are then expected to ensure an even survival of male and female offspring. As the population size decreases, however, reproductive variance starts to affect fitness and it becomes beneficial to invest more in the care of sons to dampen the stochasticity in reproductive output that results from their mortality. The clutch size

*f*has an additional, weaker, effect on equilibrium care. At the extreme of single-offspring clutches,

*f*= 1, the differences between the patterns of male and female survival are irrelevant, and equilibrium care ensures equal survival in males and females when sex ratio is equal . As clutch size increases, the effects of reproductive variance come into play and, for a given population size

*N*, larger clutch sizes result in male bias in care . However, this effect rapidly plateaus with increasing

*f*and is weaker than that of altering population size.

Theory about sex ratio predicts that to minimize the variance in offspring number, hermaphroditic females should make more offspring of the sex that is less variable in survival (Proulx 2000, 2004). In our model population, females should then produce more daughters. Using the same approach as above, we can calculate female fitness when sex ratio *r* at birth (ratio of males to total offspring) is maternally controlled; *i.e.*, the number of females at birth of a focal female with phenotype *z*_{f}* _{j}* now is

*n*

_{♀}∼ Bin(

*f*, 1 −

*r*(

*z*

_{f}

*)). As before, limited care is provided by females, and covariance in survival is greater between related males. Calculating the moments of reproduction as in Equations 13–15 with appropriate sex ratio (*

_{j}*z*

_{f}

*) and survival independent of the trait, we find as expected, (26)that the evolutionary convergent sex ratio is biased toward females .*

_{j}## DISCUSSION

### Capturing sex-specific reproductive variance

It has long been known that reproductive variance impedes adaptation by increasing genetic drift (Wright 1931; Nunney 1993; Nomura 2002; Charlesworth 2009). In parallel, a body of work has shown that reproductive variance is itself under selection, favoring less variable offspring production (Gillespie 1974; Courteau and Lessard 1999; Proulx 2000; Shpak 2005; Shpak and Proulx 2007; Lehmann and Balloux 2007; Rice 2008; Taylor 2009; Proulx and Adler 2010). Together, these studies have provided a solid theoretical basis for understanding the effects of selection on offspring distribution in a natural world that is inherently uncertain within generations. Despite these advances, models for the evolution of reproductive variance and its effects on adaptation have so far ignored biologically realistic cases of sexual reproduction, where the role of the variance can be expected to be most important. Closest to this, Taylor (2009) studied the effect of sex-specific reproductive variance on adaptation, but by modeling mating as the random union of gametes, key features of reproductive biology were neglected, since in most cases it is individuals, not gametes, that unite to mate. Finite numbers of matings and the structure of the mating system have important evolutionary effects. Not only do they generate correlations between the number of offspring of different individuals of the same sex, but they also often underlie disparities between male and female reproductive variance.

In this article, we used an individual-based approach to provide an analytical model for the evolution of male and female reproductive traits within a biologically realistic context of sexual reproduction. First, we calculated the probability of fixation of a mutant that perturbs male and female reproductive phenotypes (Equation 1), taking into account all first and second moments of the probability distribution that describes individual reproduction (Figure 1 and File S1). As a consequence, the fitness gradients, *G*_{m} and *G*_{f} (Equations 2 and 4), which express the direction and intensity of selection on a mutant, reveal the many components of reproductive variance that contribute to fitness and are hence under selection. These include the variance in the reproductive output of a focal individual (; Equations 4 and 5), which decreases fitness (Figure 2A), and the variance in the total reproductive output of the rest of the population (, Equation 4 and 5), which increases fitness (Figure 2B). The impact of these variances on fitness has been accounted for in previous studies (*Appendix*, *Link with previous work*, to see how previous works connect to the model presented here). However, our model also takes into account the covariance between the numbers of juveniles produced by different individuals of the same sex (; Equation 4 and 5), which had been ignored so far and potentially have greater consequences for fitness than the variance alone (Equation 5). This covariance would emerge as a direct consequence to the biological constraints that the number of matings and female fecundity are finite and therefore cannot be ignored.

### Efficacy of selection in males and females

The probability of fixation of a mutant (Equation 1) also depends on the efficacy with which selection can act on mutants. This is represented here by the scaling factors *K*_{m} and *K*_{f}. They measure the degree to which neutral genetic variation results in phenotypic variation is then exposed to selection in males and females (Equation 6). As *K*_{m} and *K*_{f} increase, the probability of fixation of a mutant increasingly reflects the effects it has on male and female fitnesses, respectively.

The scalars *K*_{m} and *K*_{f} express effects similar to those captured by the traditional heritability of a trait (Falconer and Mackay 1996). However, while heritability is a snapshot of a population in time, *K*_{m} and *K*_{f} take into account the segregation of alleles and changes in frequency until loss or fixation of a mutant. This is illustrated by the interpretation of *K*_{m} and *K*_{f} in terms of coalescence times (File S1, Equation SI.37), which can themselves be expressed in terms of probabilities of sibship (File S1, Equation SI.38, and Figures 3 and 4), or how genes coalesce in different individuals of both sexes. The probabilities of sibship depend on the moments of offspring production (Table 1) and thereby establish a link between the mating system of a population and the potential for selection to act on different traits in that population. High probabilities of sibship reflect a situation in which reproduction is monopolized by a subset of individuals. This reproductive skew entails a greater likelihood that a mutant is either transmitted or lost by chance and hence reduces levels of genetic variation. The factors *K*_{m} and *K*_{f} also increase with the dominance coefficient *h* of the mutant. Dominance increases the covariation between genotype and phenotype at low allele frequency, which is the frequency dominating the segregation process of a new mutation, and therefore increases the visibility of mutants to selection.

An important feature of *K*_{m} and *K*_{f} are their sex specificity, respectively scaling on male and female fitness gradients. This reveals that genetic drift can influence male and female evolution with varying strength. Traditionally, population-genetic treatments of evolution in dioecious populations express the effect of genetic drift on the segregation of two alleles simply as the inverse of the overall effective population size or as some mutant frequency-dependent function (Ethier and Nagylaki 1988; Taylor 2009), but in both cases, the effect of genetic drift on male and female selection is the same. This simplification stems from the requirements to obtain a diffusion limit for the segregation process, which ignore some differences between male and female reproduction.

The method we used here to calculate the probability incorporates all second moments of male and female individual reproduction and shows that it is possible for genetic drift to affect selection on males and females differently. When *K*_{f} is larger than *K*_{m}, selection on females contributes more to the probability of fixation than does selection on males, and vice versa. A variety of factors can lead to differences between male and female efficacy of selection. As shown in the *Example* section, discrepancies between *K*_{m} and *K*_{f} can occur as the result of differences between male and female patterns of mortality that generate a greater level of genetic variance in females than in males. This is not only a theoretical possibility. In the house finch, for example, mite ectoparasitism affects related males more strongly than related females, leading to male-biased mortality (Badyaev *et al.* 2006). As a consequence, we expect *K*_{m} to be smaller than *K*_{f} in this species.

### Long-term sex-specific evolution

To predict the joint evolution of male and female phenotypes, we embedded our model into a trait substitution sequence process. We obtained a stochastic model of long-term phenotypic evolution for dioecious populations that allows one to conveniently evaluate convergence-stable states, which correspond to the most likely phenotypic outcomes of evolution at mutation–selection–drift balance (Equations 11 and 12 and *Appendix*, *Diffusion equation for phenotypic evolution in dioecious populations*). When the reproductive variances of males and females are such that there is no difference in the efficacy of selection between the sexes (*K*_{m} = *K*_{f}; Equations 11 and 12), then the conditions for phenotypes to be convergence stable depend only on the fitness gradients, in agreement with previous deterministic models (Leimar 2009). When the efficacy of selection differs between the sexes (*K*_{m} ≠ *K*_{f}), however, they may affect the evolutionary trajectory and change the stability of internal equilibria (Equation 12). Therefore, the most likely phenotypes to be observed in natural populations can be significantly affected by sex-specific reproductive variance.

### Effects of sex-specific variance on parental care

To illustrate the many effects of reproductive variance on the evolution of dioecious species, we calculated the probability of fixation of mutants coding for maternal and paternal care for sons and daughters in a situation where the survival of sons is highly correlated within broods. While very specific, this example allows us to illustrate some of the key effects that our model can capture. First, our results demonstrate how phenotypic evolution can be driven by selection against reproductive variance. Thus, care for sons evolves more readily than care for daughters, because the former alleviates the high degree of reproductive variance that arises as a consequence of correlated male survival (Equations 19 and 20). Second, we showed that the pattern of male survival reduced the efficacy of selection on male traits by decreasing the amount of genetic variation in males (Equation 21). This means that mutants coding for maternal care have a greater probability of fixation than those coding for paternal care, even if the effect of maternal and paternal care on offspring survival is identical. Finally, male polygyny generates a negative correlation between the reproductive outputs of different males, which in turn generates an additional selection pressure on the evolution of paternal care (Equations 23 and 22). These forces mitigate the strength of selection for paternal care because as fewer males monopolize reproduction, kin competition between the offspring of a male increases. The selection pressures generated by (co)variances might be minimal when populations are very large, polygamy extensive and fecundity effectively unlimited. However, in most biologically realistic scenarios, the complicated interactions between the different components of reproductive variance can be expected to affect the evolutionary process through selection and genetic drift.

The phenotypic equilibrium predicted for both parental sexes is, like the probability of fixation, affected by selection against reproductive variance. Thus, fathers and mothers will invest more in the care of sons to mitigate the detrimental effects of their stochastic survival on parental fitness (Equation 25). Interestingly, this prediction contrasts with other results on the evolution of sex-ratio allocation whereby females are expected to produce more daughters when the survival of males within a brood is highly correlated (Equation 26) (Proulx 2000, 2004). Therefore, selection against reproductive variance leads to the counterintuitive equilibrium whereby females produce fewer sons for which they care more. It would be interesting to further explore the effect of sex-specific variance on the evolution of sex allocation. In particular, we expect that the sex ratio at birth would differ according to whether it is controlled by the male or female parent and that the difference between maternally and paternally controlled sex ratio depends on the mating system. For instance, with the life cycle given in the *Example* section of this article, if sex ratio is male controlled and the population is polygynous, then selection on males to minimize their reproductive variance would favor a bias toward females that is even more pronounced than when the sex ratio is female controlled (Equation 26).

Our analysis of the long-term evolution of male and female care behavior showed that both sexes evolve toward the same equilibrium level of care (Equation 25). This contrasts with the predicted short-term dynamics, where greater stochasticity in male survival caused a reduction in the probability of fixation of mutants for paternal care, compared to that seen for maternal care mutants. This discrepancy intuitively implies a lower rate of adaptation in the male than female trait and a longer time to reach the evolutionary equilibrium. In general, the stochastic model of phenotypic evolution (*Appendix*, *Diffusion equation for phenotypic evolution in dioecious populations*) suggests that the rate of adaptation in males and females scales with the efficacy of selection *K*_{m} and *K*_{f}, respectively.

### Outlook

The framework provided in this article is ideal for studying complex social interactions between individuals of sexual populations. Examples of such traits are those involved in evolutionary games between the male and female of a mating pair or strategies in games between individuals of the same sex, for example, in male–male competition for mating and fertilization success. In the latter case, the covariance between the numbers of offspring produced by different males is expected to have important effects. Use of our model to study the social and sex-specific frequency-dependent aspects of reproductive evolution is straightforward because all parameters in Equations 1 and 10 can be derived using only the phenotype of a focal individual and the average male and female phenotypes in the population. Another class of traits for which our model is particularly well suited are sexually antagonistic ones (Parker 1979; Lande 1980; Bonduriansky and Chenoweth 2009; Pennell and Morrow 2013). By taking into account the positive correlation of mutational effects in males and females (*b*_{mf} > 0 in Equation 10), different selection pressures in males and females (*G*_{m} ≠ *G*_{f}), and different levels of reproductive variance in the sexes, the model is well adapted to investigating the evolution of these traits under the simultaneous influences of selection and drift.

With selection on variance being inversely proportional to the population size, selection on the variance will be mostly relevant in small panmictic populations, where genetic drift may therefore mitigate its effects. But if populations are structured into local breeding groups, then selection against reproductive variance is inversely proportional to local patch size (Shpak and Proulx 2007; Lehmann and Balloux 2007), while genetic drift inversely scales with the total population size, which can be very large. All the effects of selection on reproductive variance described in this article may then be particularly relevant in populations that are divided in small patches but are globally large. In fact, if density-dependent regulation takes place before dispersal (soft selection, Roze and Rousset 2003), then selection against reproductive variance is as described by our panmictic model, with fitness given by Equation 5 but with the number of individuals being those in a local patch (Shpak 2005; Lehmann and Balloux 2007).

Explicitly taking spatial structure into account in regimes of soft and hard selection may also reveal interesting examples of sex-specific evolution. In structured and dioecious populations, we expect that sex-specific local competition (*e.g.*, Perrin and Mazalov 2000), but also reproductive variance, will drive the evolution of sex-specific dispersal. In turn, this will generate differences in genetic variation across the sexes (*i.e.*, between *K*_{m} and *K*_{f}), thereby influencing the evolution of sex-specific strategies. It would therefore be particularly interesting to extend the model to explicitly take into account spatial structure to investigate the evolution of sex-specific dispersal strategies and how it interacts with the evolution of other sex-specific traits.

Future development of the model should also accommodate for a greater variety of genetic architecture of traits. Because *K*_{m} and *K*_{f} depend on the covariance between genotype and phenotype, differences in the genetic determination of traits between the sexes would also translate into differences between the efficacy of selection in males and females. This is not unlikely as differences between male and female heritability have been reported for phenotypic traits in animals (*e.g.*, Eisen and Legates 1966; Jensen *et al.* 2003), including humans (Weiss *et al.* 2006) as well as plants (*e.g.*, Ashman 1999). Such differences would naturally arise for sex-linked genes. For instance, in species with an XY sex-determining system, where males are hemizygous for the X chromosome, dominance interactions can occur only between the two X chromosomes of females (Wayne *et al.* 2007). It would therefore be interesting to extend our model for sex-linked genes and test whether the interaction between selection on reproductive variance and the efficacy of selection in males and females lead to different evolutionary dynamics than on autosomes.

To conclude, using a population-genetic approach that takes into account all the relevant moments of reproduction in the two sexes, we have shown that the effect of sex-specific reproductive variance and covariances and selection on it influences the evolution of dioecious species. In particular, we have found that even if the fitness gradients on male and female traits have the same steepness but opposite directions, differences in male and female reproductive variance can lead to selection in one sex dominating selection in the other, and alter the trajectory of long term phenotypic evolution.

## Acknowledgments

We thank Andrew Pomiankowski, Rob Seymour, Lorette Noiret, and Julie Collet for helpful comments on the manuscript. The manuscript also vastly benefitted from the suggestions of two anonymous reviewers. C.M. was supported by a CoMPLEX Ph.D. studentship from the U.K. Engineering and Physical Sciences Research Council, M.R. by funding from the U.K. Natural Environment Research Council (NE/D009189/1, NE/G019452/1), and L.L. by the U.S. National Science Foundation (grant PP00P3-123344).

## Appendix

### Assumption on distribution of juveniles

Given an index set of individuals *i* ∈ ℐ in the population and a corresponding set of powers defined by a mapping *ζ*: ℐ → ℤ^{+}, it is assumed that the following (A1)holds, where |ℐ| is the number of individuals in set ℐ. The remainder terms that appear in *R*, given by the higher-order terms of the Taylor expansion of *F*, are thus of order 1/*N*^{2}.

### Weights for additive mutants

Using Equation SI.38 from File S1, the weights on male and female fitness gradients in the probability of fixation of an additive (*h* = 1/2) mutant are given by (A2)where

### Calculations for the evolution of parental care

Here the remaining components of the probability of fixation of alleles coding for parental care are calculated. Because male survival is different between populations in which maternal and parental care are provided (see main text), it is simpler to consider separately the cases when care is maternal and when care is paternal.

#### Maternal care:

To calculate the weight *K*_{f}, we need the probabilities of sibship (Table 1) in the resident population (*z*_{f}* _{j}* =

*z*

_{f}). To calculate the maternal probabilities of sibship Θ

^{♀}, in addition to Equation 13–15 of the main text, the covariance between number of male and female offspring that a female produces is also required, and it is given by (A3)The paternal probabilities of sibship Θ

^{♂}also influence

*K*

_{f}(Equation A2) and their components are derived below. The expected numbers of male and female offspring of a male are given by Equation 13 evaluated at male phenotype

*z*

_{m}. The variance in the number of females fathered by a male is found by conditioning on the random number of matings

*N*

_{mat}of a male (A4)Because each male is equally likely to mate, and if a male does, it mates exactly

*H*times, we have

*E*[

*N*

_{mat}] = 1 and

*V*[

*N*

_{mat}] =

*H*− 1, so that (A5)The variance in the number of males fathered by a male, given that each male brood entirely survives or dies, reads (A6)To calculate the covariance between number of males and of females produced by a male, we define

*X*as the random product of males and females coming from the same mating. Then, since offspring survival is independent across matings, we have (A7)and since, (A8)Finally, substituting Equations 13–15 and Equations A3–A8 into the probabilities of sibship (Table 1), and in turn, substituting the latter into Equation A2, we find the efficacy of selection on maternal care.

#### Paternal care:

Most of the moments required to calculate focal male fitness and the probabilities of sibship to find *K*_{m} are the same as in the preceding section (evaluated at focal male phenotype *z*_{m}* _{i}* for focal fitness, and at resident male

*z*

_{m}instead of resident female

*z*

_{f}). However, because male survival is different in a population in which parental care is provided by males, the variance in the number of male offspring fathered by the focal male, and the covariance between the number of male and female offspring fathered by a resident male (required for the probabilities of sibship) are different. Using similar arguments as above, we find that they are respectively given by (A9)Calculating the probabilities of sibship (Table 1) using the above and substituting into Equation A2, we find the efficacy of selection on paternal care.

### Diffusion equation for phenotypic evolution in dioecious populations

The substitution rate *k* (Equation 9) is the jump rate of a so-called jump process (Gardiner 2009), which here describes a population “jumping” from a monomorphic phenotypic state to another. The moments of the infinitesimal jump of the evolving phenotypes in each sex (Δ*z*_{m}, Δ*z*_{f}) characterize the distribution of phenotypic changes over an infinitesimally small evolutionary time period, and they are found by integrating the phenotypic effects of a substitution over the p.d.f. of all substitution rates, . To the first order of *δ*, we obtain for the first moment of Δ*z _{u}* for

*u*∈ {m, f}, (A10)for

*v*∈ {m, f},

*v*≠

*u*, and where

*b*=

*C*[

*δ*,

_{u}*δ*] are the second moments of mutant sex-specific effects. Because dominance is independent of

_{v}*δ*

_{m}and

*δ*

_{f}and

*K*is linear in

*h*(File S1, Equation SI.38), is simply evaluated at expected dominance Similarly, it is possible to show that the Equation for

*E*[Δ

*z*] still holds if mutants have sex-specific dominance

_{u}*h*

_{m}and

*h*

_{f}, as long as they are independent of

*δ*

_{m}and

*δ*

_{f}and that they are on average equal . The second moments of infinitesimal phenotypic change of the first order of

*δ*are given by

*E*[(Δ

*z*)

_{u}^{2}] =

*ξb*and

_{uu}*E*[(Δ

*z*

_{m})(Δ

*z*

_{f})] =

*ξb*

_{mf}.

Assuming that the phenotypic changes are continuous in probability, the first two moments of infinitesimal change Δ*z _{u}* can be used to approximate the rate of phenotypic change by a Fokker–Planck equation (Gardiner 2009). This equation characterizes the change of the distribution of male and female phenotypes

*ψ*(

*z*

_{m},

*z*

_{f};

*τ*) in evolutionary time

*τ*. The distribution of phenotypes here is meant over many evolutionary trajectories or experiments, rather than over the population. The population remains monomorphic: fixation of loss of mutants is instantaneous. If the male and female phenotypes (

*z*

_{m},

*z*

_{f}) at evolutionary time

*τ*have p.d.f.

*ψ*(

*z*

_{m},

*z*

_{f};

*τ*), then it satisfies (A11)where the functions

*a*

_{m}(

*z*

_{m},

*z*

_{f}) and

*a*

_{f}(

*z*

_{m},

*z*

_{f}) are the expected infinitesimal changes of male and female phenotypes given in Equation A10 and the stationary distribution is given by

*ψ*(

*z*

_{m},

*z*

_{f}) = lim

_{τ}_{→∞}

*ψ*(

*z*

_{m},

*z*

_{f};

*τ*).

### Link with previous work

Substituting Equation 4 into Equation 2, we find fitness gradients that are consistent with previous work on the evolution of reproductive variance. For instance, Lehmann and Balloux (2007) models the evolution of a helping trait *z*_{f} that disrupts the mean *μ _{f}*(

*z*

_{f}) ∼

*O*(

*N*) and variance in fertility. Mating is random, each female gives birth independently of one another, and sex ratio is equal at birth and in the population. Substituting Equation 4 into Equation 2, we have that the fitness gradient on

*z*

_{f}is proportional to (A12)where is the squared coefficient of variation in fertility in the resident female, in agreement with Equation A37 of Lehmann and Balloux (2007).

The original fitness gradient by Gillespie (1975, Equation 11a) or equivalently, that derived for a dioecious population by Taylor (2009, Equation 14 for an additive mutant) may be found directly from Equation A12. These analyses use the diffusion approximation, which requires that the difference between the mean fertilities of the resident and mutant phenotypes tend to zero as the population size tends to infinity, *i.e.*, that *d* ln *μ _{f}*(

*z*

_{f})/

*dz*

_{f}∼

*O*(1/

*N*). Applying this assumption to Equation A12, we have (A13)where the deleterious effects of sib competition on expected fecundity fall victim to the order condition required by the diffusion approach of Gillespie (1975) and Taylor (2009).

## Footnotes

*Communicating editor: J. Wakeley*

- Received August 3, 2013.
- Accepted October 16, 2013.

- Copyright © 2014 by the Genetics Society of America

Available freely online through the author-supported open access option.