Selection in two-sex stage-structured populations: Genetics, demography, and polymorphism

The outcome of natural selection depends on the demographic processes of birth, death, and development. Here, we derive conditions for protected polymorphism in a population characterized by age- or stage-dependent demography with two sexes. We do so using a novel two-sex matrix population model including basic Mendelian genetics (one locus, two alleles, random mating). Selection may operate on survival, growth, or fertility, any or all of which may differ between the sexes. The model can therefore incorporate genes with arbitrary pleiotropic and sex-specific effects. Conditions for protected polymorphism are expressed in terms of the eigenvalues of the linearization of the model at the homozygote boundary equilibria. We show that in the absence of sexual dimorphism, polymorphism requires heterozygote superiority in the genotypic population growth rate. In the presence of sexual dimorphism, however, heterozygote superiority is not required; an inferior heterozygote may invade, reducing the population growth rate and even leading to extinction (so-called evolutionary suicide). Our model makes no assumptions about separation of time scales between ecological and evolutionary processes, and can thus be used to project sex×stage×genotype dynamics of eco-evolutionary processes. Empirical evidence that sexual dimorphism affects extinction risk is growing, yet sex differences are often ignored in evolutionary demography and in eco-evolutionary models. Our analysis highlights the importance of sexual dimorphism and suggests mechanisms by which an allele can be favored by selection, yet drive a population to extinction, as a result of the structure and interdependence of sex- and stage-specific processes.

A.1 Population projection matrix 655 The population vector is (A-1) The population projection matrixÃ consists of 3×3 blocks, which act on the genotype-657 specific population vectors: where the explicit dependence of B onp has been omitted to avoid a cluttering of brackets.
(A-10) and apply the vec operator to both sides, remembering that asp is a vector, vecp =p, 694 dp(t + 1) = Bdp(t) + vec [I 2ωg (dB)p(t)] . (A-11) Next apply Roth's theorem (Roth, 1934), vecABC = (C T ⊗ A) vecB, to replace the vec 695 operator with the Kronecker product: 696 dp(t + 1) = Bdp(t) + (p T (t) ⊗ I 2ωg ) dvecB. (A-12) Then the first identification theorem and the chain rule together give the following formula for the Jacobian (Magnus and Neudecker, 1985;Caswell, 2007), Our aim is to express the Jacobian matrix M in terms of the genotype specific matrices, The first term in this sum, 1 T 2ωg dÃ p, is equal to zero because and therefore every column in dÃ sums to zero, see equation A-3.

712
Substituting equation (A-21) into equation (A-17) and evaluating at the boundary yields (A-25) Finally substituting the expression above into equation A-14 yields the Jacobian matrix: where we have identified the three terms as A , B , and C .

717
The next task is to work out all the terms in the above expression for the Jacobian. We Next we turn our attention to the second term, B , Using Roth's theorem, (C T ⊗ A) vecB = vecABC, we can simplify as follows Substitute the population vector analyzed at the AA boundary, To derive term C in the Jacobian, we first derive a useful expression for vecÃ in terms substituting this expression into the right hand side of equation (A-44) yields Substitutep T = (p T AA ,p T AA , 0, 0, 0, 0) into the right-hand side of equation (A-46), so that 740 only terms with j = 1 and j = 2 are nonzero, yielding Since none of the A i2 are a function of the frequency vector, Next write down each term in the sum over i and take the derivative of the vecA i1 's to Finally apply Roth's theorem (Roth, 1934), (C T ⊗ A) vecB = vecABC, to the equation 745 above (e.g. C T =p T AA , A = (e 1 ⊗ I ω ), and vecB = vecF AA ) to write this as (A-49) Equation (A-49) requires the derivative of the frequency of allele A in the gamete pool 748 with respect to the population frequency vector: (A-50) Start with equation (5) from the main text: where we can substitute p for n because of homogeneity and where the one norm can 752 be replaced by 1 T 2 W F p because p is nonnegative. For convenience, we will denote the 753 normalizing factor in the denominator with p n , to calculate Similarly for the second term in the sum in equation (A-54), Finally add equations (A-58) and (A-59) to obtain where at the boundary The Jacobian matrix, given by equation (A-64), is upper block triangular, so the eigenvalues 768 of M are the eigenvalues of the diagonal blocks. The largest absolute eigenvalue of the Ja-769 cobian, i.e. the spectral radius ρ(M), determines the stability of the boundary equilibrium. 770 We will denote the three nonzero blocks along the diagonal with M 11 , M 22 , and M 33 (see 771 equation (17)), such that for example the AA boundary, and becausep is stable to perturbations in that boundary, ρ(M 11 ) < 1.

781
The stability ofp is thus determined by which is equation (18) in the main text. The largest absolute value of the eigenvalues of the Jacobian matrix, the dominant eigenvalue, evaluated at the AA boundary, denoted by ζ AA , 784 is therefore equation (21) in the main text. By symmetry, the dominant eigenvalue of the Jacobian 786 matrix evaluated at the aa boundary, denoted byζ aa , is which is equation (22). A boundary equilibrium is unstable to invasion by the other allele 788 if the dominant eigenvalue of the Jacobian evaluated at that equlibrium is larger than one.

B Coexistence conditions under simplifying assumptions
In this section, we consider a simplification that removes sexual dimorphism in survival and 795 transition rates, U i = U i . We consider block M 22 of the Jacobian matrix, equation (A-67) where we use the following notation (equation (19) the main text), Consider an eigenvector of this matrix, with eigenvalue x, which has to satisfy the following equation We have written the eigenvector in terms of a vector associated with the female direction, Moving the terms U Aa u and U Aa u to the right in the top and bottom equations respectively 803 yields, The coexistence conditions are then given by Next we make the additional simplifying assumption that male mating success is propor-817 tional (or equal) to female fertility, F i = βF i . Finally, we also set the primary sex to 818 one, i.e. α = 0.5. We will show that these simplifying assumptions reduce the coexistence The following equalities hold in this case and therefore Substituting equations (B-21)-(B-23) into the D AA matrix yields where p b is the fraction of the total population that is in a breeding stage. Substitute this 854 expression into the definition of D AA , to obtain and equivalently, In the absence of population structure, the demographic matrices U i , U i , F i , and F i all 863 reduce to scalars, which we will label as u i , u i , f i , and f i respectively. Figure A1 shows vector directly to zero (Nussbaum, 1986(Nussbaum, , 1989. The constant population stucture can be (C-12) Although the two models have constant population structures that satisfy the same equation The population projection matrix can be written in terms of nine ω × ω matrices, As before, we define the frequency model as follows , (C-16) sulting derivative at the boundary equilibrium, 928 M = dp(t + 1) dp T (t) p . (C-17) In de Vries and Caswell (2019) it is shown that this method yields the following expression 929 for the Jacobian matrix: where we have identified the three terms as A , B , and C .

931
The next task is to work out all the terms in the above expression for the Jacobian. For

932
A and B we can simply use the results derived in de Vries and Caswell (2019), (C-20) For term C we can use the following result, where we have replaced q b A and q b a from 935 de Vries and Caswell (2019) with q A and q a , (C-21) Equation (C-21) requires the derivative of the frequency of allele A in the gamete pool with 937 respect to the population frequency vector: ( Finally substituting equation (C-28) into equation (C-21) yields  boundary, and becausep is stable to perturbations in that boundary, ρ(M 11 ) < 1.

958
The stability ofp is thus determined by These two sets of coexistence conditions are identical when 1 − α = α, i.e. when α = 0.5.