Two roads to two sexes: unifying gamete competition and gamete limitation in a single model of anisogamy evolution

Abstract

Recent studies have revealed the importance of self-consistency in evolutionary models, particularly in the context of male–female interactions. This has been largely ignored in models of the ancestral divergence of the sexes, i.e., the evolution of anisogamy. Here, we model the evolution of anisogamy in a Fisher-consistent context, explicitly taking into account the number of interacting individuals in a typical reproductive group. We reveal an interaction between the number of adult individuals in the local mating group and the selection pressures responsible for the divergence of the sexes. The same underlying model can produce anisogamy in two different ways. Gamete competition can lead to anisogamy when it is relatively easy for gametes to find each other, but when this is more difficult and gamete competition is absent, gamete limitation can provide another route for anisogamy to evolve. In line with earlier models, organismal complexity favors anisogamy. We argue that the early contributions of Kalmus and Scudo, largely dismissed as group selectionist, are valid under certain conditions. Linking their work with the contributions of Parker helps to explain why precisely males keep producing more sperm than can ever lead to offspring: sperm could evolve to provision zygotes but this brings little profit for the effort required, because sperm would have to be equipped with provisioning ability before it is known which sperm will make it to the fertilization stage. This insight creates a logical link between paternal care under uncertain paternity (where again investment is selected against when some investment never brings about genetic benefits) and gamete size evolution.

Appendix 1

Our aim is to find the fitness gradients $\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.$ and $\frac{{\partial {{\hat{W}}_x}}}{{\partial {{\hat{m}}_x}}}\left| {_{{{{\hat{m}}_x} = {m_x}}}} \right.$, the second derivatives $\frac{{{\partial^2}{{\hat{W}}_y}}}{{\partial \hat{m}_y^2}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.$ and $\frac{{{\partial^2}{{\hat{W}}_x}}}{{\partial \hat{m}_x^2}}\left| {_{{{{\hat{m}}_x} = {m_x}}}} \right.$ and the matrix ${\hbox{C}} = \left( {\begin{array}{*{20}{c}} {\frac{\partial }{{\partial {m_x}}}\left( {\frac{{\partial {{\hat{W}}_x}}}{{\partial {{\hat{m}}_x}}}\left| {_{{{{\hat{m}}_x} = {m_x}}}} \right.} \right)\quad \frac{\partial }{{\partial {m_y}}}\left( {\frac{{\partial {{\hat{W}}_x}}}{{\partial {{\hat{m}}_x}}}\left| {_{{{{\hat{m}}_x} = {m_x}}}} \right.} \right)} \hfill \\{\frac{\partial }{{\partial {m_x}}}\left( {\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.} \right)\quad \frac{\partial }{{\partial {m_y}}}\left( {\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.} \right)} \hfill \\\end{array} } \right)$ for evaluating the direction of evolution, evolutionary stability, and convergence stability respectively. We derive here the equations for the case where the mutant individual is of mating type y. Equivalent equations for type x are obtained simply by exchanging the indices x and y in all equations. A “hat” (such as in ${\hat{m}_y}$) is used to denote the mutant individual. For gamete numbers we will use the notation N _x , N _y , and ${\hat{N}_y}$ instead of the longer forms ${N_x}({m_x},{m_y},{\hat{m}_y})$, ${N_y}({m_x},{m_y},{\hat{m}_y}),$ and ${\hat{N}_y}({m_x},{m_y},{\hat{m}_y})$. However, it is important to keep in mind that these are functions of all three variables when carrying out the differentiations below.

Note that although in the main text we use fixed functions for gamete production rate, gamete mortality rate, and zygote survival probability, the equations below are given in a general form, and can flexibly accommodate alternative forms for these functions. The model also allows independent values for A _x and A _y , but in the main text we limit ourselves to cases with an even adult sex ratio (A _x  = A _y ).

All the equations that follow are derived from the time dynamics Eqs. (A1a–A1c). For their biological interpretation, see the main text concerning Eqs. (3a–3b) and (4a–4c). Despite the fact that some of the resulting equations are relatively complex and difficult to interpret by themselves, the full biology of the system is already contained in Eqs. (A1a–A1c) together with the functions for gamete production, gamete mortality, and zygote survival. The following equations are simply tools to determine the behavior of the system.

$$ \left\{ {\begin{array}{*{20}{c}} {\frac{{d{N_x}}}{{dt}} = {A_x}H({m_x}) - \mu ({m_x}){N_x} - \gamma {N_x}({N_y} + {{\hat{N}}_y}) = 0} \hfill \\{\frac{{d{N_y}}}{{dt}} = ({A_y} - 1)H({m_y}) - \mu ({m_y}){N_y} - \gamma {N_x}{N_y} = 0} \hfill \\{\frac{{d{{\hat{N}}_y}}}{{dt}} = H({{\hat{m}}_y}) - \mu ({{\hat{m}}_y}){{\hat{N}}_y} - \gamma {N_x}{{\hat{N}}_y} = 0} \hfill \\\end{array} } \right. $$

(A1a–A1c)

The fitness of the mutant individual is given by its rate of fertilizations ($\gamma {N_x}{\hat{N}_y}$, or alternatively $H({\hat{m}_y}) - \mu ({\hat{m}_y}){\hat{N}_y}$ from Eq. A1c) multiplied by the survival probability $f({m_x},{\hat{m}_y})$ of its zygotes:

$$ {\hat{W}_y} = \gamma {N_x}{\hat{N}_y}f({m_x},{\hat{m}_y}) = (H({\hat{m}_y}) - \mu ({\hat{m}_y}){\hat{N}_y})f({m_x},{\hat{m}_y}) $$

(A2)

Therefore the selection differential is

$$ \frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right. = \left[ {(H({{\hat{m}}_y}) - {{\hat{N}}_y}\mu ({{\hat{m}}_y}))\frac{{\partial f({m_x},{{\hat{m}}_y})}}{{\partial {{\hat{m}}_y}}} + f({m_x},{{\hat{m}}_y})\left( {\frac{{dH({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} - {{\hat{N}}_y}\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} - \mu ({{\hat{m}}_y})\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right)} \right]\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right. $$

(A3)

and the second derivative

$$ \begin{array}{*{20}{c}} {\frac{{{\partial^2}{{\hat{W}}_y}}}{{\partial \hat{m}_y^2}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.} = {\left[ {(H({{\hat{m}}_y}) - {{\hat{N}}_y}\mu ({{\hat{m}}_y}))\frac{{\partial {f^2}({m_x},{m_y})}}{{\partial \hat{m}_y^2}} + 2\frac{{\partial f({m_x},{{\hat{m}}_y})}}{{\partial {{\hat{m}}_y}}}\left( {\frac{{dH({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} - {{\hat{N}}_y}\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} - \mu ({{\hat{m}}_y})\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right) + \ldots } \right.} \hfill \\{\left. { \ldots f({m_x},{{\hat{m}}_y})\left( {\frac{{{d^2}H({{\hat{m}}_y})}}{{d\hat{m}_y^2}} - {{\hat{N}}_y}\frac{{{d^2}\mu ({{\hat{m}}_y})}}{{d\hat{m}_y^2}} - 2\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}}\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}} - \mu ({{\hat{m}}_y})\frac{{{\partial^2}{{\hat{N}}_y}}}{{\partial \hat{m}_y^2}}} \right)} \right]\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.} \hfill \\\end{array} $$

(A4)

Everything in equations (A3) and (A4) is straightforward to derive, with the exception of $\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}},\;\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}},\;\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}\;{\hbox{and }}\frac{{\partial {{\hat{N}}_y}}}{{\partial \hat{m}_y^2}}$. These can be found by using implicit differentiation with respect to ${\hat{m}_y}$ on equations (A1a–A1c). Differentiating once leads to the equations

$$ \left\{ {\begin{array}{*{20}{c}} {\gamma ({N_y} + {{\hat{N}}_y})\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}} + \mu ({m_x})\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}} + \gamma {N_x}\left( {\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} + \frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right) = 0} \hfill \\{\gamma {N_y}\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}} + (\gamma {N_x} + \mu ({m_y}))\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} = 0} \hfill \\{{{\hat{N}}_y}\left( {\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} + \gamma \frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}} \right) + (\gamma {N_x} + \mu ({{\hat{m}}_y}))\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}} = \frac{{dH({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}}} \hfill \\\end{array} } \right. $$

(A5a–A5c)

and after differentiating twice we have

$$ \left\{ {\begin{array}{*{20}{c}} {2\gamma \frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}\left( {\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} + \frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right) + \gamma ({N_y} + {{\hat{N}}_y})\frac{{{\partial^2}{N_x}}}{{\partial \hat{m}_y^2}} + \mu ({m_x})\frac{{{\partial^2}{N_x}}}{{\partial \hat{m}_y^2}} + \gamma {N_x}\left( {\frac{{{\partial^2}{N_y}}}{{\partial \hat{m}_y^2}} + \frac{{{\partial^2}{{\hat{N}}_y}}}{{\partial \hat{m}_y^2}}} \right) = 0} \hfill \\{2\gamma \frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} + \gamma {N_y}\frac{{{\partial^2}{N_x}}}{{\partial \hat{m}_y^2}} + \left( {\gamma {N_x} + \mu ({m_y})} \right)\frac{{{\partial^2}{N_y}}}{{\partial \hat{m}_y^2}} = 0} \hfill \\{2\left( {\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} + \gamma \frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}} \right)\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}} + {{\hat{N}}_y}\left( {\frac{{{d^2}\mu ({{\hat{m}}_y})}}{{d\hat{m}_y^2}} + \gamma \frac{{{\partial^2}{N_x}}}{{\partial \hat{m}_y^2}}} \right) + (\gamma {N_x} + \mu ({{\hat{m}}_y}))\frac{{{\partial^2}{{\hat{N}}_y}}}{{\partial \hat{m}_y^2}} = \frac{{{d^2}H({{\hat{m}}_y})}}{{d\hat{m}_y^2}}} \hfill \\\end{array} } \right. $$

(A6a–A6c)

Now we can solve $\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}},\;\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}}{\hbox{, and }}\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}$ from Eqs. (A5a–A5c):

$$ \left\{ {\begin{array}{*{20}{c}} {\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}} = \frac{{\gamma {N_x}\left( {\gamma Nx + \mu \left( {{m_y}} \right)} \right)\left( { - \frac{{dH\left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}} + {{\hat{N}}_y}\frac{{d\mu \left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}}} \right)}}{{{\gamma^2}N_x^2\mu \left( {{m_x}} \right) + \left( {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu \left( {{m_x}} \right)} \right)\mu \left( {{m_y}} \right)\mu \left( {{{\hat{m}}_y}} \right) + \gamma {N_x}\left( {\gamma {N_y}\mu \left( {{m_y}} \right) + \gamma {{\hat{N}}_y}\mu \left( {{{\hat{m}}_y}} \right) + \mu \left( {{m_x}} \right)\left[ {\mu \left( {{m_y}} \right) + \mu \left( {{{\hat{m}}_y}} \right)} \right]} \right)}}} \hfill \\{\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} = - \frac{{\gamma^2 {N_x} {N_y} \left( { - \frac{{dH\left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}} + {{\hat{N}}_y}\frac{{d\mu \left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}}} \right)}}{{{\gamma^2}N_x^2\mu \left( {{m_x}} \right) + \left( {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu \left( {{m_x}} \right)} \right)\mu \left( {{m_y}} \right)\mu \left( {{{\hat{m}}_y}} \right) + \gamma {N_x}\left( {\gamma {N_y}\mu \left( {{m_y}} \right) + \gamma {{\hat{N}}_y}\mu \left( {{{\hat{m}}_y}} \right) + \mu \left( {{m_x}} \right)\left[ {\mu \left( {{m_y}} \right) + \mu \left( {{{\hat{m}}_y}} \right)} \right]} \right)}}} \hfill \\{\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}} = - \frac{{\left( {\gamma {N_x}\left[ {\gamma {{\hat{N}}_y} + \mu \left( {{m_x}} \right)} \right] + \left[ {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu \left( {{m_x}} \right)} \right]\mu \left( {{m_y}} \right)} \right)\left( { - \frac{{dH\left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}} + {{\hat{N}}_y}\frac{{d\mu \left( {{{\hat{m}}_y}} \right)}}{{d{{\hat{m}}_y}}}} \right)}}{{{\gamma^2}N_x^2\mu \left( {{m_x}} \right) + \left( {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu \left( {{m_x}} \right)} \right)\mu \left( {{m_y}} \right)\mu \left( {{{\hat{m}}_y}} \right) + \gamma {N_x}\left( {\gamma {N_y}\mu \left( {{m_y}} \right) + \gamma {{\hat{N}}_y}\mu \left( {{{\hat{m}}_y}} \right) + \mu \left( {{m_x}} \right)\left[ {\mu \left( {{m_y}} \right) + \mu \left( {{{\hat{m}}_y}} \right)} \right]} \right)}}} \hfill \\\end{array} } \right. $$

(A7a–A7c)

and $\frac{{{\partial^2}{{\hat{N}}_y}}}{{\partial \hat{m}_y^2}}$ from equations (A6a–A6c):

$$ \begin{array}{*{20}{c}} {\frac{{{\partial^2}{{\hat{N}}_y}}}{{\partial \hat{m}_y^2}} = } \hfill \\{\frac{{ - (\gamma {N_x}\left[ {\gamma {{\hat{N}}_y} + \mu ({m_x})} \right] + \left[ {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu ({m_x})} \right]\mu ({m_y}))\left( { - \frac{{{d^2}H({{\hat{m}}_y})}}{{d\hat{m}_y^2}} + {{\hat{N}}_y}\frac{{{d^2}\mu ({{\hat{m}}_y})}}{{d\hat{m}_y^2}} + 2\left( {\frac{{d\mu ({{\hat{m}}_y})}}{{d{{\hat{m}}_y}}} + \gamma \frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}} \right)\frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right) + 2{\gamma^2}{{\hat{N}}_y}\frac{{\partial {N_x}}}{{\partial {{\hat{m}}_y}}}\left( {\gamma {N_x}\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} + \mu ({m_y})\left( {\frac{{\partial {N_y}}}{{\partial {{\hat{m}}_y}}} + \frac{{\partial {{\hat{N}}_y}}}{{\partial {{\hat{m}}_y}}}} \right)} \right)}}{{{\gamma^2}N_x^2\mu ({m_x}) + \left( {\gamma {N_y} + \gamma {{\hat{N}}_y} + \mu ({m_x})} \right)\mu ({m_y})\mu ({{\hat{m}}_y}) + \gamma {N_x}(\gamma {N_y}\mu ({m_y}) + \gamma {{\hat {N}}}\mu ({{\hat{m}}_y}) + \mu ({m_x})\left[ {\mu ({m_y}) + \mu ({{\hat{m}}_y})} \right])}}} \hfill \\\end{array} . $$

(A8)

Next, we solve equations (A1a–A1c) for N _x , N _y , and ${\hat{N}_y}$ when ${\hat{m}_y} = {m_y}$:

$$ \left\{ {\begin{array}{*{20}{c}} {{{\left. {{N_x}} \right|}_{{{{\hat{m}}_y} = {m_y}}}} = \frac{{{A_x}\gamma H({m_x}) - {A_y}\gamma H({m_y}) - \mu ({m_x})\mu ({m_y}) + \sqrt {{4{A_x}\gamma H({m_x})\mu ({m_y}) + {{\left[ { - {A_x}\gamma H({m_x}) + {A_y}\gamma H({m_y}) + \mu ({m_x})\mu ({m_y})} \right]}^2}}} }}{{2\gamma \mu ({m_x})}}} \hfill \\{{{\left. {{N_y}} \right|}_{{{{\hat{m}}_y} = {m_y}}}} = \left( {{A_y} - 1} \right)\left( {{{\hat{N}}_y}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.} \right)} \hfill \\{{{\left. {{{\hat{N}}_y}} \right|}_{{{{\hat{m}}_y} = {m_y}}}} = \frac{{ - {A_x}\gamma H({m_x}) + {A_y}\gamma H({m_y}) - \mu ({m_x})\mu ({m_y}) + \sqrt {{4{A_x}\gamma H({m_x})\mu ({m_x})\mu ({m_y}) + {{\left[ { - {A_x}\gamma H({m_x}) + {A_y}\gamma H({m_y}) + \mu ({m_x})\mu ({m_y})} \right]}^2}}} }}{{2Ay\gamma \mu ({m_y})}}} \hfill \\\end{array} } \right.. $$

(A9a–A9c)

Finally we get $\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.$ by sequentially plugging Eqs. (A9a–A9c) into (A7c), and then plugging this result together with (A9c) into (A3). Similarly $\frac{{{\partial^2}{{\hat{W}}_y}}}{{\partial \hat{m}_y^2}}\left| {_{{{{\hat{m}}_y} = {m_y}}}} \right.$ is calculated by the sequence (A9) → (A7) → (A8) → (A4).

An analytic expression for the convergence stability matrix C can be found by following steps similar to the ones described above, but both the methods and the resulting equations are more complex. We omit them from this appendix, but the equations are available from the authors upon request.

We now have the tools to find the equilibrium points and evolutionary trajectories for gamete sizes, and to determine their stability. The existence of a closed form analytic solution for the equilibrium points depends on the form of zygote fitness, gamete mortality, and gamete production functions used. With the functions we use in the main text, analytic solutions are not possible. However, we can use the following equation to find the equilibriums and plot the evolutionary trajectories to an arbitrary level of accuracy:

$$ \left( {\begin{array}{*{20}{c}} {{m_{{{x_{{i + 1}}}}}}} \hfill \\{{m_{{{y_{{i + 1}}}}}}} \hfill \\\end{array} } \right) = \left( {\begin{array}{*{20}{c}} {{m_{{{x_i}}}}} \hfill \\{{m_{{{y_i}}}}} \hfill \\\end{array} } \right) + \left( {\begin{array}{*{20}{c}} {\frac{{\partial {{\hat{W}}_x}}}{{\partial {{\hat{m}}_x}}}\left| {_{{{{\hat{m}}_x} = {m_{{{x_i}}}}}}} \right.} \hfill \\{\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = {m_{{{y_i}}}}}}} \right.} \hfill \\\end{array} } \right)\Delta $$

(A10)

where Δ is a sufficiently small number.

Once the equilibria $m_x^{*}$ and $m_y^{*}$ are found, the following conditions must be satisfied:

• 1. $\frac{{\partial {{\hat{W}}_x}}}{{\partial {{\hat{m}}_x}}}\left| {_{{{{\hat{m}}_x} = m_x^{*}}}\;} \right. = 0,\;\frac{{\partial {{\hat{W}}_y}}}{{\partial {{\hat{m}}_y}}}\left| {_{{{{\hat{m}}_y} = m_y^{*}}}\;} \right. = 0;$ this means there is no directional selection on gamete size

• 2. $\frac{{{\partial^2}{{\hat{W}}_x}}}{{\partial \hat{m}_x^2}}\left| {_{{{{\hat{m}}_x} = m_x^{*}}}\;} \right. \leqslant 0,\;\frac{{{\partial^2}{{\hat{W}}_y}}}{{\partial \hat{m}_y^2}}\left| {_{{{{\hat{m}}_y} = m_y^{*}}}\;} \right. \leqslant 0;$ this shows that the equilibrium is an ESS.

• 3. The real parts of the eigenvalues of the matrix C must be negative to show convergence stability.

Table 1 Conditions for the equilibria in Figs. 3 and 5