Mating, births, and transitions: a flexible two-sex matrix model for evolutionary demography.

Models of sexually-reproducing populations that consider only a single sex cannot capture the effects of sex-specific demographic differences and mate availability. We present a new framework for two-sex demographic models that implements and extends the birth-matrix mating-rule approach of Pollak. The model is a continuous-time matrix model that explicitly includes the processes of mating (which is nonlinear but homogeneous), offspring production, and demographic transitions and survival. The resulting nonlinear model converges to exponential growth with an equilibrium population composition. The model can incorporate age- or stage-structured life histories and flexible mating functions. As an example, we apply the model to analyze the effects of mating strategies (polygamy or monogamy, and mated unions composed of males and females, of variable duration) on the response to sex-biased harvesting. The combination of demographic complexity with the interaction of the sexes can have major population dynamic effects and can change the outcome of evolution on sex-related characters.

Introduction

Models of sexually-reproducing populations that consider a only single sex (typically females) implicitly assume that both sexes have identical vital rates and that the availability of the neglected sex (typically males) does not affect fertility (Pollard 1974; Caswell 2001; Iannelli et al. 2005). In reality, both these assumptions are frequently violated. Males and females often differ significantly in terms of fertilities and mortalities (Kuczynski 1932; Pollak 1990; Jenouvrier et al. 2010), developmental schedules (Caswell 2001), behavioral interactions (Rankin and Kokko 2007), dispersal patterns (Miller et al. 2011), and selective harvest pressures (Ginsberg and Milner-Gulland 1994). Additionally, ecological, environmental, and evolutionary changes may alter the most limiting sex over time (Hardy 2002; Miller and Inouye 2011). One-sex models miss, and cannot explore, important components of population dynamics in all these cases. Here, we introduce a flexible framework that can explore many of these factors.

In response to growing concerns over discrepancies in male and female reproductive rates (Karmel 1947), early dynamical models with sex structure were introduced in the late 1940s. Pollard (1948) used coupled Lotka integral equations, considering female births to males and male births to females in order to reconcile the growth rates of both sexes. Kendall (1949) introduced a system of ordinary differential equations for males and females (and later, married couples), and was the first to incorporate nonlinear interactions between the sexes via a mating term. Subsequent models considered other nonlinear mating functions (Pollard 1974; Yellin and Samuelson 1974) and couple dissolution through death and divorce (Hadeler et al. 1988).

Extensions to age-structured populations were made by Fredrickson (1971) and Hoppensteadt (1975), who allowed birth and death rates, as well as couple formation and divorce rates, to depend on age and sex. Hadeler later included the age (i.e., duration) of married pairs (Hadeler et al. 1988) and maturation delays (Hadeler 1993). Age-structured mating functions have similarly been proposed (Martcheva and Milner 2001). Many of these models use continuous-time equations that incorporate age structure through coupled McKendrick–von Foerster partial differential equations (Fredrickson 1971; Keyfitz 1972; Hoppensteadt 1975; Hadeler 1989, 1993), though discrete-time two sex matrix models have also been developed (Caswell and Weeks 1986).

A powerful conceptual approach to two-sex models was proposed by Pollak (1986, 1987, 1990) and called the birth matrix-mating rule (BMMR) model. It proposes that mating, births, and life cycle transition processes repeat periodically, one after the other. It contains three main components:BMMR is a useful approach for describing two-sex populations because it can specify age (and, more generally, stage) structure over all parts of the life cycle. This structure, in turn, can have significant effects on two-sex population dynamics (e.g., Sundelöf and Åberg 2006, where the addition of size-specific birth functions affects growth rate and reproductive output) and recommended management strategies (e.g., Ginsberg and Milner-Gulland 1994, where incorporating age-specific fecundity changes the outcomes of sex-biased harvest).

A mating rule function that gives the number of matings between males of age i and females of age j.

A birth matrix whose entries are the expected number of male and female offspring produced by a mating of a male of age i and a female of age j.

Sex-specific mortality schedules, which project surviving individuals to the next age class, or, in our generalization, include other stage-specific life cycle transitions.

Conceptually appealing as it is, the BMMR model has yet to be fully incorporated into the framework of stage-classified matrix population models. Our goal here is to do so, providing a general model that can be adapted to a wide range of life cycles and mating strategies. We do so by a novel extension of periodic matrix models to continuous time, based on transition rate matrices. In this paper, we focus on ecological implications of the model, but the approach can also be applied to study the evolution of sex-related traits using methods from adaptive dynamics (Shyu and Caswell 2016a, b).

Model development

For reasons that will become apparent, we have incorporated sex structure, stage structure, and life cycle processes into a continuous-time, rather than a discrete-time matrix population model. The mating, birth, and transition processes from the BMMR framework are described by separate rate matrices. The mating process introduces nonlinearity into the model because reproduction depends on the relative proportions of males and females, of appropriate stages, in the population. This dependence can be flexibly modeled by generalized weighted mean mating functions, which satisfy the biological criteria of sexual reproduction (Iannelli et al. 2005). The resulting BMMR matrix models are nonlinear but homogeneous (i.e., frequency-dependent). Models of this form generally converge to an exponential growth rate and a stable stage frequency distribution (Martcheva 1999).

Incorporating sex and stage structure

The model classifies individuals into stages based on age, developmental state, sex, reproductive status, or other variables of interest. Stage densities are projected forward in time by a projection matrix that contains the demographic rates characterizing survival, reproduction, and transitions between stages. The properties of this projection matrix provide information about the population as a whole, thereby linking individual-level life cycle information (i.e., the stage-specific vital rates in the matrix entries) to population-level properties important for ecology and evolution (e.g., growth rates or stage distributions).

A population with s stages is described by a population vector , the entries of which are the densities of each stage at time t. In a two-sex population, would contain male stages, female stages, and mated stages (unions) that could include married couples or breeding harems.

For example, the population vector for a two-sex population with mating adults and nonmating juveniles could have the form:The dynamics of the population vector are given by a system of ordinary differential equationswhere the entries of are either transition rates or rates of offspring production, and we have indicated that they may depend on the population vector.

Incorporating the BMMR processes

The BMMR framework incorporates mating, birth, and transition processes. We describe each of these processes by a separate matrix:Other life cycle processes can be included with additional matrices.

The mating (union formation) process, where adult males and females organize into reproductive unions, is described by the matrix .

The birth process, where unions produce new offspring, is described by the matrix .

The transition process, which includes other life cycle events like mortality, maturation, or divorce, is described by the matrix .

A discrete-time model would incorporate these processes as a periodic matrix product (Caswell and Shyu 2012). For example, the product would describe mating, followed by the production of offspring from the matings, followed by survival and transitions of the resulting individuals. In continuous-time, we conceptualize the processes as occurring simultaneously. It can be shown (Appendix 1) that the projection matrix in the continuous-time matrix model (Eq. 2) is the average of the transition rate matrices, e.g.,:

Modeling the mating process

The mating process, as described by the union formation matrix , depends on the relative numbers of males and females in the population, not all of which may be mature enough or available for breeding (Pollard 1974; Iannelli et al. 2005). As a result, depends on the population’s sex and stage composition, making a function of the population vector .

The total mating function gives the rate of union formation (total number of unions formed per unit time in a population of composition ). This embodies the mating rule in the BMMR framework. Here, “unions” refers to any mated, reproducing units in the population, including both one-to-one male-female pairs and harems with multiple individuals of the one sex.

We convert the total mating function into the per capita mating rates and (the average mating rates per available males m or females f respectively), The total population mating rate is .

Many commonly used mating functions are generalized weighted means (Hölder means) of the formwhere and are constants; , (Hadeler 1989; Bessa-Gomes et al. 2010; Iannelli et al. 2005). Figure 1 shows several generalized weighted mean mating functions and biologically desirable criteria that they satisfy (McFarland 1972; Pollard 1974; Yellin and Samuelson 1974).

Fig. 1

Mating functions from the generalized weighted mean family (Eq. 6) with a check to indicate which of the biologically desirable mating function criteria they satisfy

If multiple male and female stages interbreed to form different types of unions, stage-specific mating preferences can also be integrated into this mating function (Shyu and Caswell 2016b; see  Appendix 4).

It is difficult to distinguish between mating functions in populations where the sex ratio does not vary significantly (Keyfitz 1972). However, recent empirical studies (Miller and Inouye 2011) support the harmonic mean as a mating function. Because the harmonic mean satisfies reasonable biological criteria for mating functions (Caswell and Weeks 1986; Iannelli et al. 2005), and captures the qualitative properties of other generalized means, we will use a harmonic mean mating function for our analyses.

Frequency-dependent dynamics

The mating process often depends on the relative frequencies, rather than absolute abundances, of males and females. As a result, although the mating function is nonlinear, it is homogeneous of degree 1 in . That is:for any positive constant c.

As a result, the per capita mating functions (Eqs. 4 and 5) are homogeneous of degree 0 in , so that:If all entries in the projection matrix in Eq. 3 are also homogeneous of degree 0, the system is said to be frequency-dependent. This means that can be written as a function of the population frequency vector:where is the 1-norm of .

Frequency-dependent models of this type usually converge asymptotically to an equilibrium population structure (the stable stage distribution) in which all stage frequencies are constant (e.g., see Yellin and Samuelson 1974; Hadeler 1989; Martcheva 1999). The population then grows or decays exponentially at a long-term growth rate given by the dominant eigenvalue of .

To find the equilibrium stage distribution and population growth rate , it is sufficient to consider the dynamics of . It can be shown (Appendix 2) that:where is a identity matrix and is a vector of ones. One can integrate Eq. 10 with a numerical differential equation solver until population frequencies converge to . Then is the dominant eigenvalue of . The dominant right eigenvector of is proportional to the stable stage distribution .

A 5-stage BMMR matrix model

We now present an example of a BMMR matrix model with five stages: juvenile males and juvenile females , adult males and adult females , and reproducing unions u that consist of one adult male and one adult female. Single adult males and females interact to form unions, which then produce new juvenile offspring (Fig. 2). A summary of the variables, parameters, and matrices in this model is provided in Table 1.

Fig. 2

Life cycle diagram for a 5-stage population with juvenile males and juvenile females , adult males and adult females , and reproducing unions u. The functions and parameters shown here, as described in Table 1, appear in the union formation matrix (Eq. 14) (red), birth matrix (Eq. 15) (green), or transition matrix (Eq. 16) (blue). From Shyu and Caswell (2016a)

Table 1

A summary of the variables, parameters, matrices, and population properties in the 5-stage BMMR matrix model

Matrices and vectors
A	Projection matrix
B	Birth matrix (15)
T	Transition matrix (16)
U	Union matrix (14)
n	Population density vector (11)
p	Population frequency vector (9)
p^ or w	Equilibrium stage structure
Population properties
λ	Long-term population growth rate, dominant eigenvalue of A[p^]
s_1	Primary sex ratio (proportion of offspring that are born male)
s_2	Secondary sex ratio (proportion of adults that are male)
f_1,f_2	Juvenile, adult female density
m_1,m_2	Juvenile, adult male density
u	Union (mated pair) density
Life cycle parameters
α_m,α_f	Male, female maturation rates
d	Divorce rate (rate at which a male-female pair bond breaks)
h	Average harem size
k	Union reproductive rate
μ_f1,μ_f2	Juvenile, adult female mortality rates
μ_m1,μ_m2	Juvenile, adult male mortality rates
M	Total mating rate (total unions formed per time)
U_m,U_f	Per capita mating rates (4) and (5)
Harvest parameters
E	Total harvest rate in (18)
s_h	Harvested sex ratio (proportion of harvest that targets males) in Eq. 18

As in Eq. 1, we write the population vector asUsing the harmonic mean mating function, the total and per capita mating functions areAgain, we consider the life cycle in terms of mating, birth, and transition processes, which are described by matrices , , and respectively.The advantages of the continuous-time formulation are apparent at this point. The rates of maturation, mortality, and divorce in Eq. 16 combine in a simple additive fashion. Transition probabilities would not combine additively; they would involve sums of products of conditional probabilities, over all the possible events. The number of these products increases dramatically when more stages, and hence more mating combinations, are included. See Shyu and Caswell (2016b) and Appendix 4 for incorporation of multiple mating stages in the continuous-time model.

The union formation matrix contains the per capita mating functions from Eq. 12. Note that and must halved in the last row of to avoid double counting the unions formed from both male and female partners.

The birth matrix contains k, the average reproductive rate of a union, and the primary sex ratio , the proportion of offspring that are male.

The transition matrix contains the male mortality rates ( for juveniles, for adults) and female mortality rates ( for juveniles, for adults), the rates of maturation from juveniles to adults ( for males, for females), and the divorce rate d (rate at which the male-female pair bond breaks). Note that unions may also dissolve due to partner death, and that union dissolution through both death and divorce returns surviving males and females to the single adult stages.

As per Eq. 3, the average of these three matrices is the continuous-time projection matrix . The corresponding equation for the proportional structure (Eq. 10) is thus:where is given by Eq. 14, is given by Eq. 15, and is given by Eq. 16.

As shown in Fig. 3, the population vector ultimately grows exponentially, while the frequency vector ultimately converges to the constant distribution of stages . To determine the equilibrium stage distribution , we integrated Eq. 17 with the Matlab ODE45 differential equation solver until population frequencies converged (e.g., until vector entries do not change significantly over consecutive integration intervals).

Fig. 3

Dynamics of the 5-stage BMMR model with monogamous unions and no harvest. The population consists of juvenile males and juvenile females , adult males and adult females , and reproducing unions u. a) Growth of the population density vector (11), where dynamics are given by (3). b) Convergence of the population frequency vector (9), where dynamics are given by (17). Parameters are fixed at , .

The Matlab code used for all the following examples is provided in the Electronic Supplementary Material.

Mating systems and the effects of sex-biased harvest

As an example of the use of the two-sex model, we consider the effectes of sex-biased harvesting. In sport or trophy hunting, harvest is often male-biased and significantly exceeds natural mortality (Festa-Bianchet 2003; Milner et al. 2007). Age or size bias is also common, as larger or older males with well-developed adult characteristics (e.g., large antlers or horns) are usually the most desirable targets. This unnatural selection may alter population structure, reproductive strategies, body morphology, and developmental timing (e.g., Ashley et al. 2003; Festa-Bianchet 2003; Allendorf and Hard 2009). The population response depends on multiple demographic factors, including the mating system.

The mating system determines how males and females organize for breeding and is thus a key component of two-sex population structure (Emlen and Oring 1977). Some species form monogamous, one-to-one pair bonds between males and females. Other species have polygynous mating systems in which one male mates with multiple females (e.g., lions, deer, seals), or, more rarely, polyandrous systems where one female mates with multiple males (e.g., jacanas, pipefish). Unions formed by mating may be transient and limited to a single breeding episode (e.g., lek systems) or may persist over multiple breeding seasons (e.g., lion harems) or even until partner death (e.g., albatrosses and other species with high mate fidelity) (Cézilly and Danchin 2008).

These factors motivate the use of a demographic two-sex model in analyzing harvest effects. To this end, we will use our BMMR matrix framework to explore the effect of mating systems on the response to sex-biased harvest. A range of mating systems will be approximated by varying two model parameters, d and h (Table 2).

Table 2

Mating systems corresponding to different values of the divorce rate d and harem size h

	Low d (persistent unions)	High d (transient unions)
h=1 (monogamy)	Persistent pair bonds, high mate fidelity (e.g., albatross)	Serial pair bonds (e.g., humans, Emperor penguins)
h>1 (polygyny)	Persistent harems (e.g., lion prides)	Leks, scramble competition (e.g., grouse, cod, horseshoe crabs)

The divorce rate d. This is a measure of union transience. Unions with higher values of d are more likely to dissolve after a given mating, while unions with lower values of d are more likely to persist over multiple breeding interactions.

The harem size h. This is a measure of polygamy. Unions with are monogamous and consist only of one-to-one male-female pair bonds, while unions with are polygamous groups of size . As polyandrous mating systems are relatively rare (Cézilly and Danchin 2008), we will consider only the polygynous form of polygamy, where one male mates with h females.

Harvest strategies are characterized by the overall harvest rate E and the harvested sex ratio (proportion of the total harvest rate that targets males). Assuming that only adults are harvested, the adult mortality rates are modified by harvest as follows:To determine how various mating systems, as characterized by different values of h and d, respond to sex-biased harvest, we will examine harvest effects on the long-term population growth rate and the secondary sex ratio (proportion of all adults that are male). We will assume that males and females have the same baseline vital rates, and that the primary sex ratio (proportion of males at birth) is 0.5. Thus, the main sex-specific differences we consider are sex-biased harvest pressures and male versus females roles within the polygynous mating systems.

Monogamy ()

Consider a monogamous two-sex population with juveniles and adults. The mating process forms unions that are one-to-one pair bonds of adult males and females. This scenario uses the rate matrices , , and as given by Eqs. 14, 15, and 16 respectively, and the mating functions in Eq. 12. We vary the divorce rate d to explore the effects of transient vs. persistent pair bonds.

As shown in Fig. 4a, the proportion of adults in the reproductive union stage (mated adults) decreases as d increases. The unharvested population growth rate similarly decreases (Fig. 4b), because populations with more transient couples (fewer mated adults) cannot produce offspring as rapidly as populations with more persistent couples (more mated adults). Because males and females have the same baseline vital rates, the secondary sex ratio remains unbiased (, not shown) for all values of d.

Fig. 4

Population growth rates, structure, and responses to adult harvest in the monogamous () model. a) Proportion of mated adults (adults in stage u, rather than or ) for an unharvested population, and the b) corresponding population growth rates , as functions of the divorce rate d. c) The change in and d) the change in secondary sex ratio , when a proportion of harvest targets males. Without harvest, for all values of d. Other parameters are fixed at

Figure 4c shows how increasingly sex-biased harvest strategies affect population growth. Both biased and unbiased harvest strategies most strongly reduce growth in populations with lower divorce rates, as adult mortality will also disrupt pair bonds. The greatest decreases in occur when harvest is strongly sex-biased, i.e., is close to 0 (only females are harvested) or 1 (only males are harvested). This suggests that monogamous populations with high fidelity pair bonds will be the most impacted by sex-biased harvest, and that concentrating harvest on a single sex will more greatly reduce population growth.

Figure 4d shows how the same harvest strategies decrease the secondary sex ratio from equality. Unsurprisingly, male-biased strategies reduce , female-biased strategies increase , and unbiased strategies () have relatively minimal effect. Populations with greater divorce rates experience larger reductions in , possibly due to their lower growth rates (Fig. 4a) being unable to replenish harvested individuals as rapidly.

Polygyny ()

Now consider a polygynous population that forms unions consisting of one male with a harem of females. Because the death or departure of a single female changes the harem’s size and reproductive rate, we must now account for multiple union (harem) stages.

The stage represents polygynous unions consisting of one male and a harem of i females. When h is the maximum harem size, the population vector contains h union stages, which range from a union with a harem of size 1 (, equivalent to a monogamous couple) to a union with a harem of size h (, the largest possible union):Assume that when males and females mate, they form the largest possible union . Their union formation rate is still given by the harmonic mean mating function Eq. 12, but with the number of single females now replaced by the number of prospective harems:Note that the total union formation rate is maximized when the sex ratio of single adults isThus, as the harem size h increases, a higher proportion of single females is needed to maximize the mating rate.

If an individual female has a reproductive rate of k, a harem with i females has a total reproductive rate of ik; larger harems are thus more productive. Each union , regardless of harem size, can change in three possible ways (Fig. 5):As a result, a union may shrink (but not grow) in size due to the departure or death of its members (Fig. 6). After a union shrinks to the smallest possible size , or if the male harem leader dies, the union dissolves and its members return to the stages for unmated adults.

Fig. 5

Reproduction and three possible transitions for , a union with harem size i

Fig. 6

Stages in a population with maximum harem size h, which include juveni le males and juvenile females , adult males and adult females , and reproducing unions (one male with a harem of i females). Adults form harems of size h when mating, and these harems can shrink in size, but not grow, over time

The male harem leader dies (with mortality rate ). This returns i adult females to the stage .

A female harem member dies (with mortality rate ). This shrinks the union from to . For the union , which has only one female, corresponds to the single adult male stage (i.e., the death of a wife returns her husband to the pool of unmated singles).

A female harem member departs from the union, with divorce rate d. We assume that only females leave, which seems biologically plausible. Her departure returns one female to and shrinks the union from to . For the union , divorce dissolves the union and returns one male to and one female to .

Appendix 3 shows how to write the rate matrices , , and for a polygynous system with maximum harem size h. The population vector and matrices for the case where are as follows:

Figure 7 shows how the unharvested population rate of growth and secondary sex ratio vary across different values of d and h. As in the monogamous case, lower divorce rates result in more mated reproducing adults and, thus, higher population growth. Larger harems lead to more rapid population growth, possibly because of their higher total reproductive rates. Unions with a high maximum harem size are more resilient to divorce and female mortality, because they can lose more females before shrinking to a and then dissolving.

Fig. 7

Unharvested population dynamics in the polygynous () union model, as functions of the divorce rate d and harem size h. a) Population growth rate . b) Secondary sex ratio (proportion of males in the adult population). Other parameters are the same as in Fig. 4

Even without sex-biased harvest, the secondary sex ratio is slightly female-biased (), but varies only a few tenths of a percentage point across a wide range of h and d. Populations with high h and low d (large, persistent harems) are the most biased.

Figure 8 demonstrate how female-biased (), unbiased (), and male-biased () harvest strategies affect population growth rates and secondary sex ratios.As expected, the secondary sex ratio increases during female-biased harvest, decreases during male-biased harvest, and undergoes only minimal changes when harvest is unbiased (Fig. 8, right). Populations with high d and low h experience the largest sex ratio shifts under biased harvest. This may be because the smaller growth rates of high d, low h populations are less effective in offsetting harvest-induced sex ratio biases.

Fig. 8

Population responses to harvest that is a female-biased (), b unbiased (), and c male-biased () in the polygynous () model. (Left) The change in population growth rate . (Right) The change in secondary sex ratio . Other parameters are the same as in Fig. 4

Female-biased harvest (Fig. 8a) most strongly reduces growth in populations with large d and small h, the same populations with the lowest unharvested growth rates (Fig. 7a). Smaller harems are less resilient to female-biased harvest for the same reason they are less resilient to divorce and female mortality - because they cannot lose as many females before dissolving. Female-biased harvest may reduce the average harem size, and also makes it difficult for large harems to form or reform after breaking up. When h is high, a higher proportion of females is needed to maximize the mating rate, in accordance with Eq. 21, but these females are depleted by harvest. Increasing the divorce rate increases the rate at which harems dissolve. This more drastically reduces growth for larger harems (contours are steeper at large h), because they need more females to reform.

Unbiased harvest (Fig. 8b) yields similar qualitative trends. Again, populations with higher divorce rates experience greater reductions in growth, and larger harems are more affected by divorce. The effect of increasing divorce is not as pronounced as with female-biased harvest (contours are flatter overall), as less female harvest makes it easier for harems to reform. At low h, however, populations with lower d are actually more impacted by harvest. Low h unions have only a few females and are more likely to dissolve from increased mortality. Unions with high divorce rates are already dissolving quickly, regardless of harvest mortality. Unions with low divorce rates, in contrast, break up much more frequently once harvest mortality occurs. As a result, low d, low h populations experience the largest decreases in growth.

Male-biased harvest (Fig. 8c) reverses the effects of increased divorce rate. Focusing harvest on males is more likely to dissolve harems by killing their male leaders. Populations with low d experience the largest reductions in growth, because male-biased harvest makes these unions dissolve more frequently than they normally would (similar to the low d, low h case for unbiased harvest). As in the previous scenarios, the growth of large h populations is less affected by harvest. Even though male-biased mortality causes unions to break up more frequently, it also returns (potentially many) females to the pool. This may be beneficial when h is high, as a higher proportion of females is needed to maximize the mating rate in accordance with Eq. 21.

Figure 9 shows how the growth rates of the mating systems in Table 2 vary with harvest bias and intensity. Again, we see that high h, low d populations (large, persistent harems) have the largest growth rates of all the mating systems, even under harvest. Low h, high d populations (small, transient harems) have the smallest growth rates.

Fig. 9

Growth rates as a function of the total harvest rate E for populations with various mating systems. The four types of mating systems shown correspond to different harem sizes h and divorce rates d (Table 2); in this example, low , high , low , and high . Harvest may be female-biased (), unbiased (), or male-biased (). Other parameters are the same as in Fig. 4

Increasing the total harvest rate E in Eq. 18 amplifies the differences between female-biased, unbiased, and male-biased harvest strategies. Female-biased harvest (red) decreases population growth more severely than male-biased harvest (blue) does, even across populations with different h and d. This may be because there is an excess of single males waiting to become harem leaders, whereas single females are usually in shorter supply (especially when h is large). Additionally, the death of a male harem leader immediately dissolves his union; this can return many females to the singles stage, allowing new, full-sized harems to reform with a new male leader. The death of female harem members, in contrast, does not necessarily cause the union to dissolve, and may instead generate small, stunted harems with reduced productivity.

Depending on the mating system, unbiased harvest (black) can decrease growth rates more or less than female-biased harvest does. In populations with low h (small harems), female-biased harvest has the most drastic impacts of all the harvest strategies — again, perhaps, because small harems cannot afford to lose as many females before dissolving. In populations with high h (large harems), however, unbiased harvest may be just as, if not more, detrimental to population growth.

Discussion

The framework we present here is a tool for studying the effects of sexual reproduction, mating systems, and life histories on population dynamics. Our implementation of the BMMR approach places no limitation on the complexity of the life cycle, and is flexible enough to accommodate age or stage structure, diverse mating systems, and mating preferences. Because it is formulated as a matrix model, powerful sensitivity analysis tools are available (see Shyu and Caswell 2016a, b).

There are interesting open problems to which the approach can be applied. As formulated here, the model describes a constant environment. Time-varying (seasonal or stochastic) models would be interesting and challenging extensions. Spatial models connecting multiple subpopulations might have important effects on mate limitation. And, because the model is formulated in terms of pair formation, it will be interesting to see how it might apply to species that do not form pairs, but in which sex ratio may influence population dynamics (e.g., pollen limitation in plants).

In our application to sex-biased harvest we found that mating factors, including harem size and union duration, affect not only unharvested population growth, but also the responses of growth rate and sex ratio to sex-biased harvest. In unharvested populations, high rates of divorce, which lead to more transient unions, tend to reduce population growth, especially when harem size is small. Sex-biased harvest affects not only population sex ratios, but also long-term growth rates, with effects depending on sex bias, harem size, and divorce rate. These complex, and sometimes counterintuitive, nuances would be impossible to capture without a demographic two-sex model like this, motivating the use of such models in ecological studies.

Our two-variable depiction of the mating system, defined by harem size and union duration, can be extended to other factors. Mating systems may differ in parental investment by males and females. While polygynous males tend to provide minimal parental care, monogamous males invest on par with their female partners (Emlen and Oring 1977; Cézilly and Danchin 2008). Sex-biased harvest may have different consequences for offspring survival in these mating systems. Other species have additional nuances in how they respond to sex-biased harvest; African lions, for instance, commit infanticide when male harem leaders are killed (Whitman et al. 2004), which exacerbates the effects of male harvest on population growth.

How populations respond to selective harvest also has important consequences for evolution. Growing evidence suggests that evolutionary considerations are relevant to sustainable long-term management (Ashley et al. 2003), and that human-induced selection is especially important for harvested species. As harvest mortalities are often more severe and selective than natural mortalities, they may drive evolution in directions that would not occur under natural conditions.

Because it integrates life cycle structure, sex ratio, and a mating function, the approach introduced here is ideally suited to studying the evolution of sex ratio as a component of life history evolution. Sex ratio evolution has a long and distinguished history in evolutionary demography (e.g., Darwin 1871; Fisher 1930; Trivers 1972; Charnov 1982, among many others). Many of these discussions invoke demographic factors, such as relative mortality of males and females, but do so without a demographic model that incorporates those factors.

Our approach provides such a model, and can be analyzed using the framework of adaptive dynamics (e.g., Geritz et al. 1998), which accounts for the frequency-dependence inherent in the mating process. We have applied this to the evolution of sex ratio as influenced by sex-biased offspring costs (Shyu and Caswell 2016a) and multiple maternal conditions (Shyu and Caswell 2016b). In these studies, the dynamics of the population frequency vector are used to evaluate the possibility for a mutant phenotype to invade a resident phenotype; phenotypes that are not invasible are singular strategies; their stability properties lead to the identification of evolutionarily stable (ESS) or convergence stable strategies. See Shyu and Caswell (2016a, b) for details.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary material 1 (ZIP 354 KB)