Research articlesModelling the impacts of male alternative reproductive tactics on population dynamics
Abstract
Observations of male alternative reproductive tactics (ARTs) in a variety of species have stimulated the development of mathematical models that can account for the evolution and stable coexistence of multiple male phenotypes. However, little attention has been given to the population dynamic consequences of ARTs. We present a population model that takes account of the existence of two male ARTs (guarders and sneakers), assuming that tactic frequencies are environmentally determined and tactic reproductive success depends on the densities of both types. The presence of sneakers typically increases overall population density. However, if sneakers comprise a sufficiently large proportion of the population—or, equivalently, if guarders are sufficiently rare—then it is possible for the total population to crash to extinction (in this extreme regime, there is also an Allee effect, i.e. a threshold density below which the population will go extinct). We apply the model to the example of the invasive round goby (Neogobius melanostomus). We argue that ARTs can dramatically influence population dynamics and suggest that considering such phenotypic plasticity in population models is potentially important, especially for species of conservation or commercial importance.
1. Introduction
Alternative reproductive tactics (ARTs) represent a taxonomically widespread biological phenomenon characterized by the coexistence of two or more discrete phenotypes that achieve reproduction in very different ways. The existence of two phenotypically distinct male reproductive types (with behavioural, morphological and physiological differences) was first observed in field crickets, Gryllus integer [1] and bluegill sunfish, Lepomis macrochirus [2]. Males often use large size, showy colours or other costly displays to compete for, or court, females. In some species, there is a distinct male type that lacks these morphological, physiological or behavioural traits, and—using ‘sneaky’ or coercive interloper tactics—exploits males that court and/or provide parental care [3].
The coexistence of two or more competing reproductive tactics in a single population challenged the assumption that there is one ‘best’ reproductive tactic, and changed the view of how reproduction evolves [4,5]. As examples of discontinuous behavioural and morphological variation in reproduction, ARTs provide an excellent opportunity to develop an evolutionary understanding of phenotypic plasticity, and to shed light on evolutionary and ecological processes in general [3,5–9].
Past efforts to model ARTs have focused predominantly on understanding why and how selection might favour the evolution of ARTs, and how these alternative phenotypes are maintained within a population (reviewed in [3,5,8]). Game theory has been employed to explore conditions for long-term evolutionary stability and to understand the evolutionary trade-offs between alternative tactics [4,5,10–13]. Phenotypic variation in tactics is thought to arise either (i) via polymorphic genotypes (with equal fitness) coexisting due to balancing selection or (ii) via a monomorphic genotype in which condition or status dictates which of the possible tactics (with unequal fitness) will be adopted [3,9,14,15]. Conditional strategies are common in nature; often an environmentally determined trait such as growth rate or body size at a critical age or time will cue an individual to employ one tactic or another [3,5,15–17].
In situations where ARTs arise from genetic polymorphism, the evolutionary dynamics can be important in the short term. For example, Myers [10] modelled the evolutionary dynamics of ARTs in Atlantic salmon, in which the two types are sexually precocious small parr (sneakers) and larger adult males. Because the evolutionarily stable strategy (ESS) depends on mortality rate, changes in fishing pressure can shift the ESS so the frequency of precocious parr changes, with potentially substantial consequences for the yield of a salmon fishery [10,18].
By contrast, if ARTs arise from a conditional strategy that is genetically monomorphic in the population, evolutionary changes will occur on much longer timescales, so evolution can be ignored when considering short-term population dynamics. To our knowledge, the short-term population dynamics of species with condition-dependent ARTs have not previously been modelled.
In this study, we develop a simple ordinary differential equation (ODE) model to describe the population dynamics of a species in which there are two competing reproductive tactics, the frequencies of which are determined by (slowly changing) environmental conditions. Our aim was to improve general understanding of the role of ARTs in population dynamics, and also to enhance our specific understanding of how ARTs may affect the population dynamics of one particular species, the round goby (Neogobius melanostomus), an invasive fish in the Laurentian Great Lakes and parts of Europe. We consider how the presence of more than one reproductive male tactic influences the expected overall population density, and how this overall density can be expected to change in response to environmentally induced changes in the proportions of the population displaying each of the ARTs.
Understanding the short-term impacts of ARTs on the population biology of a commercially important species, or a species of conservation concern, can improve predictions of abundance patterns and extinction risk (as has been explored for a species with alternative life-history strategies that do not involve distinct reproductive tactics [19,20]). If, for example, male fish employing different reproductive tactics are distinguishable by radically different body sizes, then selectively harvesting one male type (either accidentally or intentionally) could have unexpected, undesirable or even irreversible effects on the population (as Myers [10,18] inferred for Atlantic salmon based on an evolutionary analysis).
Our model can be used to obtain qualitative insights about the population dynamics of many species that express similar ARTs, the frequencies of which are environmentally determined. However, our work was initially motivated by the invasive round goby and we will illustrate our results using parameter values estimated for this (or related) species.
2. The round goby
Round goby, originally from the Ponto-Caspian region, were accidentally released from ship ballast water into Lake St Clair, which is attached to Lake Erie (one of the Laurentian Great Lakes), around 1990 [21]. Since then, round goby have rapidly expanded in both range and abundance to all five of the Great Lakes [22,23] and (independently) to Western Europe [24]. The successful invasion of the round goby constitutes a triple threat: (i) they can out-compete native fish species for food, shelter and breeding habitat [25–30], (ii) they eat eggs and young of other species [31,32], and (iii) by virtue of their capacity to eat bivalves, they appear to be contributing to toxicant transfer in areas of contamination [29,33–36]. Round goby are thought to have contributed to the decline of many native species and to the deterioration of ecosystem health in general [24,28,37]. A great deal of time, effort and money has been spent on preventing round goby from expanding further [38–40].
Round goby exhibit male ARTs [41,42]. Smaller sneaker males exploit the effort of larger, nest-guarding males by sneaking into nests and fertilizing the eggs within, thus avoiding the energy expenditures of both guarding the nest and caring for young [41,43]. These ARTs have been studied in round goby in both fresh and brackish waters, e.g. Lake Ontario [41,42,44,45], the Rhine, Elbe and Danube rivers, and the Bay of Lübeck [46–48]. ARTs have also been investigated in many other gobies, including the common goby, Pomatoschistus microps [49], the black goby, Gobius niger [50] and the sand goby, Pomatoschistus minutus [51]. However, to our knowledge, the impacts of ARTs on population dynamics—of goby and most other species with ARTs—have not previously been investigated.
3. Model of population dynamics with alternative reproductive tactics
We refer to the males that provide parental care as guarders and the other males as sneakers. We denote the density of guarders by G, the density of sneakers by S, and the total (male) density by F = G + S. We use the symbol F because we are motivated by a fish species; we do not model females directly, but assuming a 1 : 1 sex ratio, F is also the density of females.
We are interested in situations in which reproductive tactic choice is determined by environmental variables, such as water temperature, food availability or the degree of hormone-altering contamination. We introduce the model parameters below.
3.1. Intrinsic reproductive rate of guarders
If a single female lays her clutches at a guarder’s nest, and the expected clutch size and interspawn interval are E and T, respectively, then the (per capita) instrinsic reproductive rate of guarders—in the absence of sneakers—will be
3.2. Sneaker probability
We suppose that when a male enters the breeding population he will become a sneaker with probability σ (or a guarder with probability 1 − σ). We assume that σ is environmentally determined (not genetically determined and not determined by population density), so it can be considered constant in our population dynamic model. If environmental conditions were to change (either naturally or as a result of human impacts), the probability σ would change.
3.3. Mortality rate
We assume for simplicity that guarders and sneakers have the same natural death rate (μ), which ensures that the proportion of the population that is made up of sneakers (the sneaker prevalence) will quickly converge to σ (since sneakers enter the population in this proportion). Thus, we can assume that
In practice, guarders might live longer than sneakers and have a different age of maturity [52] (but see [53]); however, the two male types may nevertheless have similar reproductive lifespans, which is what actually affects the dynamics of our model.
3.4. Sperm competition from sneakers
Sneakers hinder reproduction of guarders through sperm competition [6,54,55]. The proportion of fertilizations obtained by sneakers depends on the number of sneakers at a guarder’s nest, and at the population level on the average number of sneakers per guarder (S/G, equation (3.4)). We refer to the proportion of fertilizations that are obtained by guarders as the guarders' share and denote it by . Exactly how depends on the number of sneakers per guarder does not affect our mathematical analysis, but when relating results to sneaker prevalence we typically assume that more sneakers per guarder reduces the guarders' share (i.e. decreases with S/G or σ).
It is important to bear in mind that when we refer to sperm competition in the context of our model, the competition is between sperm of males that differ only in genes that do not influence reproductive strategies. The model assumes that offspring display the same (conditional) reproductive strategy regardless of whether they were fathered by guarders or sneakers.
3.5. Competition between guarding males for breeding space
If the density of guarders [G = (1 − σ)F] increases, then their individual reproductive success can be expected to decrease due to competition for suitable breeding habitat. For sufficiently high densities, competition for space will be so severe that the density of guarders will decrease. We formalize the density dependence of guarder reproduction with a (standard) factor
3.6. Guarder reproductive rate
Taking account of the intrinsic reproductive rate (νg), sperm competition from sneakers () and competition among guarders for nest space (3.6), the per capita rate at which guarders produce new male offspring (a proportion σ of which will become sneakers) is
3.7. Sneaker reproductive rate
The expected share of eggs obtained by an individual sneaker is the sneakers’ total share of eggs () divided by the expected number of sneakers at a nest (S/G, equation (3.4)). Thus, the rate at which a single sneaker is expected to produce offspring from a given nest is
3.8. Model equation
Combining equations (3.8) and (3.12), and recalling that G = (1 − σ)F and S = σF (equation (3.3)), the total reproductive rate is
Considering the linearization of equation (3.16) for small populations (x ≪ 1),
3.9. Parameter estimates
Table 1 lists the parameters of the model (3.13), together with our best estimates of their values for the round goby in Hamilton Harbour. Similarly, table 2 lists the three parameters of the dimensionless version of the model (3.16). Our parameter estimation is described in appendix A.
| parameter | symbol | estimate | definition |
|---|---|---|---|
| intrinsic reproductive rate of guarders | νg | 1.47 ± 0.98 yr−1 | equation (3.1); rate at which guarding males produce offspring that survive to maturity (without resource limitations or competition) |
| mortality rate | μ | 0.65 (0.46, 0.85) yr−1 | deaths per unit time per capita (same rate for both male types) |
| sneaker proportion | σ | 0.33 (0.15, 0.50) | fixed proportion of males that become sneakers, independent of tactic of parent |
| guarders’ share | 0.92 (0.76, 0.98) | proportion of eggs fertilized by guarders | |
| threshold guarder density | — | equation (3.6) | |
| breeding habitat capacity (guarder carrying capacity) | Kg | — | equation (3.9); if G > Kg then guarder density declines in the absence of sneakers |
| challenge to sneaker success | equation (3.11); guarding male density required such that the probability that sneakers find nests to invade is 1/2 |
| parameter | symbol | estimate | definition |
|---|---|---|---|
| guarder intrinsic lifetime reproductive success | ν | 2.25 ± 1.64 | equations (3.15b), (A 3) |
| guarders’ share | 0.92 (0.76, 0.98) | proportion of eggs fertilized by guarders | |
| challenge to sneaker success (scaled version of ) | 1.47 ± 0.96 | equations (3.15a), (A 5), (A 6) |
4. Model analysis
In order to understand how the expected population dynamics depends on parameter values, we examine the equilibria of the model (3.16) and their stability. Our analysis is summarized in table 3.
| properties of equilibria | region of parameter space |
dynamics | bifurcation | ||
|---|---|---|---|---|---|
| 1 | x− < 0 < x+ | 1 < ν, | certain persistence | ||
| 2 | ν+ ≤ ν, | ||||
| 3 | 0 < x− < x+ | 1 ≤ ν ≤ ν+, | Allee effect | ||
| 4 | , | ||||
| 5 | x− = 0 < x+ | , | persistence (marginal) | transcritical | |
| 6 | 0 = x− = x+ | , | extinction (marginal) | transcritical, saddle node | |
| 7 | x− < 0 = x+ | , | extinction (marginal) | transcritical | |
| 8 | 0 < x− = x+ | 1 < ν ≤ ν+, | Allee effect (marginal) | saddle node | |
| 9 | 0 < x− = x+ | 1 < ν ≤ ν+, | Allee effect (marginal) | saddle node | |
| 10 | x− = x+ < 0 | ν− ≤ ν < 1, | certain extinction | saddle node | |
| 11 | x± complex | ν− ≤ ν ≤ ν+, | certain extinction | ||
| 12 | x− < x+ < 0 | , | certain extinction | ||
4.1. Existence of equilibria
We first note that the model is well posed, i.e. that it does not predict negative populations: since the population state variable x is a factor of the right-hand side of equation (3.16), the rate of change of x when x = 0 cannot be negative, so positive initial states always yield positive solutions.
To find equilibria, we set dx/dτ = 0 and solve for x, which yields three potential equilibrium points , where
In order to specify the various cases concisely in table 3, we define
4.2. Stability of equilibria
Because our model (3.16) is one-dimensional and dx/dτ is a continuous function, a complete dynamical stability analysis is straightforward. Local stability of an equilibrium at a point is determined by the sign of dx/dτ on either side of (in situations like ours in which dx/dτ is differentiable, we can make use of the sign of d2x/dτ2 at ). Local stability of all equilibria determines all basins of attraction (which comprised the segments between equilibria). This is how we have characterized the dynamics in the various parameter regions listed in table 3.
Figure 1 shows dx/dτ as a function of x. The nine panels—all of which were drawn using the same values of ν and —differ in the guarders’ share , and illustrate how changes in can alter the equilibria and dynamics. Where dx/dτ is positive (negative) the population will increase (decrease). Regardless of parameter values, sufficiently large populations will decrease (in equation (3.16), dx/dτ < 0 for all x > x+ if x+ is positive and for all x > 0 if x+ is non-positive or complex).
Figure 1. The possible dynamical behaviours of our population model (3.16), which incorporates the effects of two alternative reproductive tactics (ARTs) in males, i.e. guarding or sneaking. In the panel labels, [i] refers to row i in table 3. Increases and decreases in population density x (horizontal axis) are determined by the sign of dx/dτ (vertical axis). The nine panels correspond to different guarders’ share ; the other parameters are fixed at ν = 3.2 and = 1.86. The hatched region in each panel is not biologically relevant (x < 0), but is shown so changes in the global dynamics and bifurcations of the model are more evident. Grey arrows indicate the direction in which x is changing. The axis scales are the same across each row of panels, but differ in each row.
The marginal cases shown in the middle column of figure 1 are important for dynamical understanding since they correspond to bifurcations of the model; however, because they correspond to parameter combinations for which two equilibria exactly coincide they will not occur in practice. Parameter combinations that yield bifurcation points are indicated in the final column of table 3. In general, there are three biologically distinct dynamical possibilities:
4.2.1. Case 1: certain persistence (figure 1a,b)
If a single positive equilibrium exists (x− ≤ 0 < x+), then all initially positive populations tend to x+ ( is unstable or at least unstable from above). (The collision of x− with 0 causes a transcritical bifurcation [58–60].)
4.2.2. Case 2: bi-stability (Allee effect; figure 1c,d,i)
If two positive equilibria exist (0 < x− < x+) then there is a critical population density (x−), below which extinction is certain ( is stable) and above which persistence and convergence to the largest equilibrium (x+) is certain. Thus, there is an Allee effect [61]. There is also a marginal Allee case, in which the two non-zero equilibria are equal and positive (0 < x− = x+); in this situation (figure 1e,h), the unique positive equilibrium is semi-stable: all populations larger than x+ tend to this point, and all others go extinct. (The collision of x− with x+ causes a tangent or saddle-node bifurcation [58–60,62].)
5. Model interpretation
If the basic reproduction number (3.18) is greater than one (), then persistence is guaranteed even if it is extremely challenging for sneakers to reproduce (i.e. even if (3.15a) is very large). If then extinction is possible, and the more challenging it is for sneakers to reproduce (the larger ) the larger the parameter region in which the population will certainly go extinct. Similarly, the less challenging it is for sneakers to reproduce (the smaller ) the easier it is for sneakers to compensate for the fact that guarders cannot sustain the population on their own, and hence the larger the parameter region in which the population can persist via an Allee effect.
Figure 2 displays the different possible dynamical regimes in the versus ν parameter plane. In the main (top left) panel, the challenge to sneaker success is set to its estimated value (, table 2). The other (smaller) panels show how the parameter regions that yield different dynamics vary with . The region of certain persistence (green) does not depend on . However, if sneaking is easy (small ) then there is a large parameter region where an Allee effect occurs (blue), whereas if sneaking is hard (large ) then the Allee region is very small and the certain extinction region is correspondingly large. Figure 2 can be understood in biological terms as follows:
| — | In each panel of figure 2, the right boundary corresponds to the situation where guarders obtain all the eggs (). On this boundary, the basic reproduction number (3.18) is simply , so the population persists if ν > 1 (guarders more than replace themselves) and the population goes extinct if ν < 1 (guarders do not manage to replace themselves). | ||||
| — | Away from the right boundary, sneakers obtain some of the eggs (), so the condition ν > 1 is not sufficient to guarantee that the population as a whole replaces itself. Instead, the certain persistence condition () can be written , so the left boundary of the green region is the curve . Within the green region, even though guarders obtain only a share of the eggs, their intrinsic lifetime reproductive success ν is sufficient that the population will be sustained no matter how small the sneaker reproductive rate. | ||||
| — | Outside of the green region, extinction is always possible because . In the grey region (certain extinction), sneakers are unable to compensate for the inability of guarders to replace themselves. In the blue region (Allee effect), the total population will go extinct if it starts from too low a density, but it can persist if the population density exceeds x− (equation (4.1)). | ||||
Figure 2. Dynamical behaviour of model (3.16) as a function of the guarders’ share (, §3.4) and the guarder intrinsic lifetime reproductive success (ν, equation (3.15b)). Different areas of the (, ν) parameter plane correspond to one of three possible behaviours described in §4.1: certain persistence, an Allee effect (persistence only above a threshold density), or certain extinction. Each panel has the challenge to sneaker success (, equation (3.15a)) fixed at a different value (larger x1/2 values imply fewer sneaking opportunities). The black dot and error bars in panel (a) show our estimate of (, ν) for the round goby in Hamilton Harbour (cf. table 2). To connect with the analysis summarized in table 3, within the region in which extinction is certain, we use a dotted curve to indicate the subregion where x± are complex. The dynamics of the model are driven ultimately by the environmentally determined proportion of the population that employs the sneaker tactic (σ, §3.2). The sneaker proportion σ determines the proportion () of eggs that are fertilized by guarders. if is sufficiently small () then the population can collapse and go extinct (see §5).
Figure 3 shows how the scaled equilibrium population densities (x±) vary as a function of the guarders’ share . Each panel corresponds to a different horizontal line (fixed ν) in figure 2, and colours of curves indicate the corresponding stability region in figure 2.
Figure 3. Equilibrium densities (x−, x+, and 0) as functions of the guarders’ share , with ν and fixed (see equation (4.1)). Solid curves denote stable equilibria; where positive, they are coloured according to the corresponding region in figure 2. Dotted black curves denote unstable equilibria. In all six panels, there is a transcritical bifurcation, where the stable equilibrium at x = 0 changes to unstable (solid grey line at 0 changes to dotted; also marked by a dot on the upper (x+) curve). In panels (a), (b) and (c), there are also two saddle node bifurcations (stable and unstable equilibria collide and disappear: blue or green solid lines vanish when they intersect the dotted black lines below them). The transition from panel (c) to (d) illustrates a saddle node bifurcation that occurs as ν changes with fixed. Each panel corresponds to the horizontal line at the indicated value of ν in figure 2a, and together they illustrate the possible ways that total population density can change as sneaker proportion (σ) changes (which causes the guarders’ share to change; e.g. (5.2)). Our best estimates for the round goby (table 2) yield panel (a).
It is easier to interpret the results if we express equilibria in units of the total male density F (guarders plus sneakers) (obtained from equation (3.14a)) relative to the sneaker-free equilibrium Kg (equation (3.9)),
Figure 4. Equilibrium densities as functions of the guarders’ share , as in figure 3, but plotted using the more easily interpretable scale of total male density (F) relative to the guarder carrying capacity (Kg). To compute the equilibrium curves, we assumed the specific relationship between and the sneaker prevalence (σ) given in equation (5.2), and set in equations (5.2) and (5.3). Line styles and colours have the same meanings as in figure 3.
From figure 4, we can infer population density effects that can be expected to occur if environmental changes cause the prevalence of sneakers (σ) to increase (and the guarders’ share to decrease according to equation (5.2)). If there are no sneakers (σ = 0, ), the equilibrium density will be F+/Kg = 1, whereas if environmental conditions change in a way that causes at least some individuals to opt to sneak (σ > 0, ) then the equilibrium total density F+ will increase, as is evident from the green curves in all panels of figure 4. Exactly what is expected as sneaker prevalence is increased further depends on the guarder intrinsic lifetime reproductive success ν (it also depends on the challenge to sneaker success , which has the same value in all panels of figure 4). Eventually, if sneaker prevalence continues to increase to very high levels, one of two things can happen: either the population will certainly crash and go extinct (figure 4a–c) or the total population density will continue to rise, but with the danger of extinction due to an Allee effect (figure 4d–f).
In the interesting—if unlikely—extreme of very large sneaker prevalence (hence near 0), where there is an Allee effect, figure 4 shows that the total density F can greatly exceed the guarder carrying capacity Kg. However, the density of guarders is only (1 − σ)F, which will be much smaller.
5.1. Application to Hamilton Harbour
The parameter estimates listed in table 2 place the Hamilton Harbour round goby system in the certain persistence region (Case 1 of §4.2 with x− ≃ −1.40 and x+ ≃ 0.54); see the heavy dot in figure 2a. Thus, if we take our best estimates (table 2) at face value then the model (3.16) predicts that the round goby population in Hamilton Harbour will persist and should approach an equilibrium density at which the total male density F+ (guarders plus sneakers) (equation (3.14a)) relative to the sneaker-free equilibrium Kg (equation (3.9)) is
Of course, tables 1 and 2 also indicate large uncertainties in estimated parameter values. Considering these uncertainties, we can at best suggest that the Hamilton Harbour goby population probably corresponds to some point in the region defined by the error bars shown in figure 2a. This plausible parameter region does dip into the (grey) certain extinction region in figure 2a (and, in fact, within the grey region of all panels of figure 2, since the boundary of the green region is independent of ). However, a much larger proportion of the area encompassed by the error bars is green; consequently, within the limitations of the model, it is reasonable to conclude that—in the absence of environmentally induced changes in or ν—the round goby is likely to persist in Hamilton Harbour.
5.2. Affecting population size and dynamics
Our approach to modelling the population dynamics of a species with ARTs was motivated by the prevailing view that expression of a tactic is state- or condition-dependent, controlled by the environment (e.g. by food or nest availability, or by temperature; [5,63]). Our model (3.16) provides a way to forecast (or at least qualitatively understand) the possible population dynamic effects of environmental changes that are either naturally, accidentially, or intentionally induced.
5.2.1. Observed population dynamics
The population dynamics of the round goby in Hamilton Harbour have been studied for more than two decades and have been described previously [42,64,65]. The population density of round goby appeared to be declining from 2002 to 2008, but more recently seems to have stabilized and might reflect a (noisy) equilibrium. The initial apparent decline over the study period might reflect existing predators adjusting to the presence of round goby in the lake and increasing their consumption of this new prey species. If environmental conditions were roughly constant after the decline—and hence parameters of our model could be considered to be unchanging—our model would predict convergence to an equilibrium, which is roughly consistent with the noisy apparent equilibrium that is observed.
5.2.2. Population dynamic effects of pollution
Previous work has indicated that the prevalence of sneakers is greater in more contaminated sites [36], perhaps as a result of endocrine disruption. Suppose that contaminant exposure controls the proportion of sneakers in the population (σ) without affecting any other parameters of our model, and that σ always increases as a function of contaminant concentration. If a goby population is currently at the equilibrium associated with our best estimates for the parameters (table 2) then the possible population density changes that can be induced by changes in contaminant levels correspond to moving left from (towards a smaller guarders’ share ) in figures 2a, 3a or 4a. In particular, figure 4a indicates that continually increasing pollution should eventually cause the population to decline and go extinct. If, rather than focusing on our best estimate, we were to consider the uncertainty in ν (table 2), then figure 4f could be more relevant; in that case, pollution could push the population into the (blue) Allee region, implying that we would expect the population to persist and continually increase in density. In this scenario, high contaminant levels could induce a crash to extinction only if random fluctuations (which we have not modelled) caused the population density to fall below the persistence threshold (density x−).
5.2.3. Intentional population control
Since the round goby is invasive in Lake Ontario, if human actions were to cause the species to go extinct, that outcome would probably be considered beneficial. Of course, achieving that goal by polluting the lake is not desirable! Figures 2 and 4 show that there could be great value in identifying strategies that would not harm the ecosystem and yet could (i) greatly increase the prevalence (σ) of sneakers, (ii) reduce the expected lifetime reproductive success (ν) of guarders in the absence of sneakers, and/or (iii) make it harder for sneakers to succeed (increase ). All three of these conditions would probably be satisfied if a commercial fishery were established, since there would be a strong bias for removing larger individuals. A second approach would be to encourage sport fishing, but—unlike usual catch-and-release practices—to require permanent removal of large fish.
6. Discussion
To our knowledge, this study represents the first attempt to model the population dynamic impact of alternative reproductive tactics (ARTs), the frequencies of which are environmentally determined. Our model (3.16) is sufficiently simple that we were able to carry out a complete dynamical analysis, which showed that the presence of sneakers can substantially alter the expected population density and susceptibility to extinction.
The outcomes predicted by our model (figure 2) depend on three parameters: the intrinsic lifetime reproductive success of guarders in the absence of sneakers (ν, equation (3.15b)), the proportion of fertilizations obtained by guarders (, §3.4), and how challenging it is to sneak successfully (, equation (3.15a)). The population will persist provided the basic reproduction number (, equation (3.18)) is greater than one. The proportion of the male population that adopts the sneaker tactic (σ, §3.2) or, equivalently, the proportion (1 − σ) that provides parental care, ultimately determines the guarders’ share () and consequently controls persistence versus extinction. Our mathematical analysis (§4) does not depend on how is related σ, but by making an assumption about this relationship we can infer how the prevalence of sneakers—or of caring males—affects overall population density.
Assuming that a guarder will fertilize fewer eggs ( smaller) if there are more sneakers at the nest (σ larger), we find that greater sneaker prevalence (larger σ) is typically associated with increased total population density, unless sneaker prevalence is so high—i.e., guarder prevalence is so low—that susceptibility to extinction is greatly increased (figure 4). Sufficiently high sneaker prevalence—high enough that , so and guarders cannot sustain themselves (§5)—induces extinction risk: extinction either becomes certain (grey in figure 2) or becomes possible at very high sneaker prevalence due to an Allee effect (i.e. a density threshold that must be exceeded for the population to persist; blue in figure 2).
In our model, effects of high sneaker prevalence are equivalent to effects of low guarder prevalence, but similar effects in real systems could be generated by guarders abandoning nests that are overwhelmed by sneakers, or simply by guarders investing less in parental care because they perceive lower paternity [66], e.g. if they recognize fewer offspring as their own. (Note that in the context of our model, recognition would have to be based on genes unrelated to reproductive strategies, since we assume the population is monomorphic with respect to reproductive strategies.) The relationship between paternity and parental care is complex [67]. In the specific context of the round goby, offspring recognition and the impact of paternity on parental care have not been studied, but they have been explored in other gobies: sand goby have been shown to recognize whether offspring are their own and alter care accordingly [68,69], whereas the common goby has been found not to have this ability [70].
Consistent with our model (figure 4), a positive correlation between sneaker prevalence and total population density could arise from fertility enhancement; more eggs might be successfully fertilized in nests with multiple males [71–73]. At the same time, an Allee effect could be generated by high sneaker prevalence (large σ) simply because the proportion of males providing parental care to young (1 − σ) is small.
Our model predicts conditions under which populations persist or go extinct (table 3). Our estimates of the model parameters for the round goby in Hamilton Harbour (table 2) suggest strongly that if environmental conditions do not change then this population will persist (§5.1).
6.1. Limitations and future directions
Our modelling approach has several limitations that could be addressed in future studies:
6.1.1. Density-dependent tactic choice
We have assumed that the sneaker probability σ can be treated as constant, meaning that changes in σ occur slowly compared with the rate at which equilibrium density is reached. This assumption is reasonable for the kinds of situations that motivated us, namely tactic choice determined by environmental conditions that change on timescales of many generations of the focal organism.
One might argue that tactic choice should depend on population density, since at low density it will be more difficult for sneakers to find a guarder with a nest to parasitize. However, any density-dependent challenge in finding guarders will be faced equally by females. Consequently, the expected number of sneakers per guarder—or, equivalently, the expected number of sneakers per nest in which a female has laid eggs—might be similar at high and low density, implying that the density changes that occur as equilibrium is approached might not affect tactic choice.
However, some goby species have been observed to make density-dependent tactic choices [74,75]. For these species at least, the probability (σ) of adopting the sneaker tactic does depend on density (and hence so does the guarders' share, ). There are some species that display an ontogenetic shift, whereby (small) young males start as sneakers and then grow into (larger) guarders [76–78]. In the case of the round goby, it is unknown if tactic is set for life or if a sneaker may become a guarder in the next season [45]. While our model provides a useful starting point—and begins to close the gap in the theory of how ARTs influence population dynamics—addressing the effects of density-dependent tactic choice will be important to consider in future modelling work.
6.1.2. Combined effects of genes and environment
We have assumed that the probability that an individual becomes a sneaker is determined entirely by the environment. However, whether a male ends up as a sneaker or guarder morph may depend on genetics, the environment or (most likely) on some combination of the two [5,17,79]. It would be valuable to consider a model in which tactics are inherited from parents (e.g. in ruffs, Philomachus pugnax, and marine isopods, Paracerceis sculpta, ARTs have been shown to have an inherited component [80,81]), or emerge as an outcome of some kind of gene-environment interaction [17] where the switch is genetically controlled [82].
6.1.3. Parameter estimation
In addition to a more complex and more realistic model, having firmer parameter estimates for round goby would yield greater predictive power and increase confidence about the hypothesized fate of the population in Hamilton Harbour. In particular, the estimate of degree of paternity loss to sneakers (i.e. , where is the guarders’ share) could be improved as this was based on estimates for a relative, the sand goby, and not the round goby itself [83]. It is worth noting that our estimate for this parameter () is similar to the degree of paternity gained by plainfin midshipman sneakers (; [84]), but much lower than the degree of paternity gained by bluegill sunfish sneakers (; [85–87]) or by mature salmon parr (; [88–90]). Further investigation is needed to determine the average level and range of fertilizations obtained by round goby guarders and sneakers in various habitats.
6.1.4. Seasonality and spatial dynamics
Our current, simple model (3.16) does not allow for variance in reproductive success across time or space. Including seasonal forcing in the model could capture the real variation observed in breeding success across the season (April to October in the case of the round goby [42,64,65]). More generally, exploring how both spatial and temporal variation in reproductive success influences population parameters would be of great interest, especially for a fish species like the round goby, which can tolerate many different environments and has rapidly expanded its range into many different ecological niches and habitats [24,91,92].
6.1.5. Invasion dynamics
Previous modelling of round goby population dynamics [93] has focused on invasion dynamics in the absence of ARTs [94–96]. Here, we have focused on the effects of ARTs on the dynamics of established populations. Potential effects of ARTs on the invasion process remain to be investigated.
6.2. Conclusion
In general, there is a growing need to incorporate our knowledge of mating systems, adaptive phenotypic plasticity, and variation in mating behaviour into population models. Modelling these behavioural processes will help us answer fundamental questions in ecology and might lead to better control and conservation management strategies.
Ethics
This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility
All data used to parametrize the model are taken from previous publications, as described in appendix A.
Declaration of AI use
We have not used AI-assisted technologies in creating this article.
Authors' contributions
J.A.M.Y.: conceptualization, data curation, formal analysis, investigation, methodology, software, writing—original draft; S.B.: conceptualization, data curation, funding acquisition, investigation, project administration, resources, supervision, writing—review and editing; D.J.D.E.: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing—original draft, writing—review and editing.
The study was conceived jointly by the three authors. J.A.M.Y. wrote the first draft, an early version of which appeared in her MSc thesis (which was jointly supervised by D.J.D.E. and S.B.). J.A.M.Y. tragically died before completing a sequence of planned revisions. D.J.D.E. developed the model and analysis further, and revised the manuscript. D.J.D.E. and S.B. discussed and edited all versions. D.J.D.E. and S.B. gave final approval for publication and agreed to be held accountable for the work performed therein.
Conflict of interest declaration
We declare we have no competing interests.
Funding
We were funded by NSERC discovery grants (to S.B. and D.J.D.E.), and project grants from Fisheries and Oceans Canada (to S.B.).
Acknowledgements
We are grateful to the large number of students and colleagues with whom we have collaborated on round goby research over the years, but especially to Marten Koops, Julie Marentette, Erin McCallum, Adrienne McLean, Natalie Sopinka and Sina Zarini. We are also grateful to Erol Akçay and Erin Okey for valuable comments on the manuscript.
Footnotes
1 τ is time expressed in units of mean lifetime. To see that the mean lifetime is 1/μ, note that if the intrinsic reproductive rate were zero (νg = 0), so the only process occurring was natural death, then equation (3.13) would reduce to dF/dt = −μF, the solution of which is F(t) = F0 e−μt. Thus, if F0 individuals enter the population at time 0, at any future time t there will be F0 e−μt still alive, implying that the distribution of lifetimes is exponentially distributed with mean 1/μ.
Appendix A. Parameter estimates
To determine the dynamics that we might expect to observe in a round goby population, we estimate a mean and plausible range for each of the model parameters (table 1). Whenever possible we use data from our own studies of the round goby population in Hamilton Harbour, Lake Ontario [36,42,64]. For remaining parameters we use estimates from other published work, using studies of the round goby when possible, and studies of other related goby species otherwise.
A.1. Guarder intrinsic reproductive rate in the absence of sneakers (νg)
We are interested in the number of offspring that survive to maturity per caring male per unit time. We need to estimate E, T, p, b and r in equation (3.1).
A.1.1. Clutch size, E
MacInnis & Corkum [97] counted ripe eggs in the ovaries of 136 females in the Upper Detroit River, and found that mean fecundity was 198 eggs, with a range of 84 to 606 eggs.
A.1.2. Interspawn interval, T
Eggs develop in 14–15 days at , or in 18–20 days at 17−19°C [98]. The mean near-shore temperature, averaged over data from the 2004–2019 breeding seasons in Hamilton Harbour, was approximately , and so we estimate—via linear interpolation—that eggs in Hamilton Harbour take an average of 17.5 days to hatch, which we interpret as the mean interspawn interval. We have not incorporated recovery or courting time between clutches, so we may be overestimating the reproductive potential of guarding males.
A.1.3. Survival probability, p
Charlebois et al. [98] estimated that 95% of the eggs in a clutch are fertilized, and 95% of those are successfully hatched. We estimate that the probability that juveniles survive to maturity is 0.0053, following Table 1 of Vélez-Espino et al. [93], who collected abundance data in Hamilton Harbour between 2002 and 2008 and used it to estimate the mean annual juvenile survival at four different sites in Hamilton Harbour. Hence, we find
A.1.4. Breeding season length as a proportion of the year, b
In Hamilton Harbour, males exhibit reproductive characteristics from late April to late October [36,42,64], suggesting approximately 190 days of suitable breeding conditions. Therefore, we estimate b = 190/365 ≃ 0.52.
A.1.5. Proportion of males that are reproductive, r
Reproductive status of males in Hamilton Harbour from 2006 to 2017 was analysed by McCallum et al. [42]; based on Table S1 from that paper, we estimate that the annual proportion of all males that showed reproductive status was r = 0.29 ± 0.08.
Inserting the above estimates for E, T, p, b and r in equation (3.1), we find νg = 1.47 male offspring per guarder yr−1. To obtain an error Δνg, we apply the Delta method [56,57] to equation (3.1),
A.2. Mortality rate (μ)
Vélez-Espino et al. [93] estimated the mean annual survival (as a proportion) for the round goby population in Hamilton Harbour from 2002 to 2008 to be a = 0.52 (range 0.42–0.62). Converting to an instantaneous mortality rate via standard survival analysis (cf. endnote 1) we have , i.e. μ = log(1/a) and Δμ = −(1/a)Δa, from which we infer μ = 0.65 yr−1 (range 0.46–0.85 yr−1).
A.3. Sneaker proportion (σ)
From 2006 to 2017 round goby were sampled in Hamilton Harbour, and males were categorized as guarder, sneaker or non-reproductive [41,42,64]. Of reproductive males, the mean annual proportion of sneaker males over the six sites was 0.33 (range 0.15–0.50), which we use as an estimate of σ.
A.4. Dimensionless guarder intrinsic lifetime reproductive success (ν)
Inserting our estimated μ in equation (3.15b), we find ν = 2.25. Using the Delta method,
A.5. Guarders’ share
Jones et al. [83] used microsatellite DNA analysis to determine rates of sneaking in a natural sand goby population (Pomatoschistus minutus). Of nests that had sneakers present, the average proportion of assayed eggs that each guarder male had successfully fertilized was (range 0.76–0.98). (Jones et al. [83] observed considerable variance in the number of sneakers in individual nests; in their DNA analysis they found four instances of one clutch having three fathers, i.e. the parental male, plus two sneaker males.)
A.6. Challenge to sneaker success (G½, x½)
If the probability that a sneaker successfully finds a nest to parasitize (3.11) is P, and we take the threshold guarder density (equations (3.6) and (3.15a)) to be , and the system is in equilibrium (4.1), then
