Neutral hybridization can overcome a strong Allee effect by improving pollination quality

Abstract

Small populations of plant species can be susceptible to demographic Allee effects mainly due to pollen limitation. Although sympatry with a common, co-flowering species may somewhat alleviate the problem of pollinator visitation (pollination quantity), the interspecific pollen transfer, IPT, (pollination quality) may remain a barrier to reproduction in small populations such as new introductions. However, if the two species are crosscompatible, our hypothesis is that neutral hybridization can help the small founding population overcome the Allee effect by improving the quality of pollination. We tested this hypothesis by using a novel modelling approach based on the theory of kinetic reactions wherein pollinators act as enzymes to catalyse the reaction between the two substrates: pollen and unselfed ovule. Using a single locus, two-allele genetic model, we developed a generic model that allows for hybridization between the invading and the native genotypes. Analysing the stability properties of the trivial equilibria in hybridization model as compared with the single genotype invasion model, we found that hybridization can either remove or reduce the Allee effect by making an otherwise stable trivial equilibrium unstable. Our study suggests that hybridization can be neutral but still be the key driver of a successful invasion by mediating pollen limitation. Conservation programmes should therefore account for this cryptic role that hybridization could play in plant invasions.

Appendices

Appendix A: Derivation of the models

In this section, we use a quasi-steady state approximation followed by rescaling of the dependent and independent variables to obtain a dimensionless ordinary differential equation that describes the rate of change in seed number as a function of pollen and ovule density; the dynamics are defined by Eqs. 1, 2 and 3.

A.1 Kinetic reaction theory and quasi-steady-state approximation

The time-evolution of N can be obtained through a quasi-steady-state approximation as described in Keener and Sneyd (2009) and Murray (2002). Let s ₁=[S ₁], s ₂=[S ₂], e=[E], c ₁=[C ₁], c ₂=[C ₂], n=[N] be the abundance of pollen, ovules, pollinator without pollen, pollinators with pollen, pollinator with pollen on ovules and seeds, respectively. We derive, from the kinetic reaction (1), (2) and (3), the following system of ODEs:

$$\left\{\begin{array}{lllllllll} \frac{ds_{1}}{dt}=-k_{1}s_{1}e,\\ \frac{ds_{2}}{dt}=-k_{0}s_{2}-k_{3}c_{1}s_{2}+(1-\nu)k_{4}c_{2},\\ \frac{de}{dt}=-k_{1}s_{1}e+k_{2}c_{1},\\ \frac{dc_{1}}{dt}=k_{1}s_{1}e-k_{2}c_{1}-k_{3}c_{1}s_{2}+k_{4}c_{2},\\ \frac{dc_{2}}{dt}=k_{3}c_{1}s_{2}-k_{4}c_{2},\\ \frac{dn}{dt}=k_{0}s_{2}+\nu k_{4}c_{2}. \end{array}\right. $$

(16)

Note that

• n will be known once c ₂ and s ₂ are known and it does not affect the dynamics of the other variables in the system so the equation can be removed from the system,

• $\displaystyle {\frac {de}{dt}+\frac {dc_{1}}{dt}+\frac {dc_{2}}{dt}=0}$ thus e + c ₁ + c ₂ = e(0) + c ₁(0) + c ₂(0) = e ₀ and e can be deduced from c ₁ and c ₂, so we can remove e from the system by replacing it with e ₀−c ₁−c ₂. Rewriting the equations for s ₁, s ₂, c ₁ and c ₂, we obtain the following new system

  $$\left\{\begin{array}{llll} \frac{ds_{1}}{dt}=-k_{1}s_{1}e_{0}+k_{1}s_{1}c_{1}+k_{1}s_{1}c_{2},\\ \frac{ds_{2}}{dt}=-k_{0}s_{2}-k_{3}c_{1}s_{2}+(1-\nu)k_{4}c_{2},\\ \frac{dc_{1}}{dt}\,=\,k_{1}s_{1}e_{0}-k_{1}s_{1}c_{1}-k_{1}s_{1}c_{2}-k_{2}c_{1}-k_{3}c_{1}s_{2}+k_{4}c_{2},\\ \frac{dc_{2}}{dt}=k_{3}s_{2}c_{1}-k_{4}c_{2}. \end{array}\right. $$

  (17)

The quasi-steady state approximation assumes that the rate of formation and breakdown of complexes are essentially equal at all times, that is d c ₁/d t=0 and d c ₂/d t=0. This approximation can be derived through non-dimensionalization of the system, choosing for example the following dimensionless variables:

$$\tau=k_{1}e_{0}t,\: u_{1}=\frac{s_{1}}{s_{0}},\: u_{2}=\frac{s_{2}}{s_{0}},\: v_{1}=\frac{c_{1}}{e_{0}},\: v_{2}=\frac{c_{2}}{e_{0}},$$

where e ₀ and s ₀ are the initial number of enzyme (pollinator) and substrate (pollen and ovule), respectively.

Now, writing the equations based on the scaled variables, one gets:

$$\left\{\begin{array}{lllllll} \frac{du_{1}}{d\tau}=-u_{1}+u_{1}v_{1}+u_{1}v_{2},\\ \frac{du_{2}}{d\tau}=-\frac{k_{0}}{k_{1}e_{0}}u_{2}-\frac{k_{3}}{k_{1}}v_{1}u_{2}+(1-\nu)\frac{k_{4}}{k_{1}s_{0}} v_{2},\\ \epsilon\frac{dv_{1}}{d\tau}=u_{1}-u_{1}v_{1}-u_{1}v_{2}-\frac{k_{2}}{k_{1}s_{0}}v_{1}-\frac{k_{3}}{k_{1}}v_{1}u_{2} +\frac{k_{4}}{k_{1}s_{0}}v_{2},\\ \epsilon\frac{dv_{2}}{d\tau}=\frac{k_{3}}{k_{1}}v_{1}u_{2}-\frac{k_{4}}{k_{1}s_{0}}v_{2}, \end{array}\right. $$

(18)

with $\epsilon =\frac {e_{0}}{s_{0}}$.

Notice that in our framework 𝜖<<1 as that there is much more pollen and ovule (substrate) than pollinator (enzyme) in the beginning.

Thus, letting 𝜖→0 while considering the last two equations in Eq. 18, we have

$$\left\{\begin{array}{lllllll} 0=u_{1}-u_{1}v_{1}-u_{1}v_{2}-\frac{k_{2}}{k_{1}s_{0}}v_{1}-\frac{k_{3}}{k_{1}}v_{1}u_{2}+\frac{k_{4}}{k_{1}s_{0}}v_{2},\\ 0=\frac{k_{3}}{k_{1}}v_{1}u_{2}-\frac{k_{4}}{k_{1}s_{0}}v_{2}. \end{array}\right.$$

which after solving for v ₁ and v ₂, we obtain the following:

$$\left\{\begin{array}{lllllll} v_{2}=\frac{k_{3}s_{0}}{k_{4}}v_{1}u_{2},\\ v_{1}=\frac{u_{1}}{u_{1}+\frac{k_{3}s_{0}}{k_{4}}u_{1}u_{2}+\frac{k_{2}}{k_{1}s_{0}}}. \end{array}\right. $$

(19)

Scaling back to the initial variables c ₁ and c ₂, Eqs. 19 can be written

$$\left\{\begin{array}{lllllll} c_{1}=\frac{e_{0}s_{1}}{\frac{k_{2}}{k_{1}}+s_{1}+\frac{k_{3}}{k_{4}}s_{1}s_{2}},\\ c_{2}=\frac{k_{3}}{k_{4}}c_{1}s_{2}. \end{array}\right. $$

(20)

Note that these two equations could have been obtained directly by setting d c ₁/d t=0 and d c ₂/d t=0. Using the last equation in (16):

$$\frac{dn}{dt}=k_{0}s_{2}+\nu k_{4}c_{2}.$$

and after some substitutions and simplifications, we obtain the following ODE that describes the dynamics of seed production n as related to the density of pollen, s ₁, and ovules, s ₂:

$$ \frac{dn}{dt}=k_{0}s_{2}+\frac{V_{\max}s_{1}s_{2}}{Q_{1}Q_{2}+Q_{2}s_{1}+s_{1}s_{2}}, $$

(21)

where

$$V_{\max}=\nu e_{0}k_{4},\:Q_{1}=\frac{k_{2}}{k_{1}},\:Q_{2}=\frac{k_{4}}{k_{3}}.$$

We also obtained equations for the temporal dynamics of s ₁ and s ₂ but our focus here is on the role of the ovules and pollen on the production of seeds and not the dynamics of these reactants per se.

The above Eq. 21 models the birth rate of seeds only: to account for death of individuals, we incorporated a density-dependent mortality component as follows:

$$\frac{dn}{dt}=k_{0}s_{2}+\frac{V_{\max}s_{1}s_{2}}{Q_{1}Q_{2}+Q_{2}s_{1}+s_{1}s_{2}}-d_{1}n -d_{2}n^{2}, $$

(22)

with d ₁ and d ₂ being the density-independent and the density-dependent death rates, respectively (Harper and McNaughton 1962).

This last Eq. 22 describes the rate of change in n, the density of seeds, as a function of n itself but also that of the pollen and ovule densities. The first term in the right hand side of the Eq. 22 represents reproduction through self-fertilization, whereas the second term indicates seed production via outcrossing. The last term on the right-hand side of Eq. 22 accounts for death that can happen at any stage of the plant’s life cycle over the whole growing season. Note that Eq. 22 models changes in plant density over time as a single stage, so that from a mathematical perspective, n can be envisioned both as a seed in the seedbank or as a mature plant above the ground. That is, the rate of change in n is simply the result of the difference between birth and death rates. This simplification is reasonable as our model is based on the life cycle of an annual plant with no persistent seedbank and thus with no overlap between generations. We can therefore further simplify the model by assuming that the other two dynamic variables, i.e. s ₁ and s ₂ in Eq. 22, are proportional to n (this results in the final model (4) as will be discussed later).

A.2 Non-dimensionalization

To reduce the number of free parameters in the system, we scaled the variables and reduced the model to a dimensionless problem (Buckingham 1914; Segel 1972) by strategically choosing

$$t^{*}=td_{1},\quad n^{*}=\frac{d_{2}}{d_{1}}n,\quad s_{1}^{*}={c_{s}^{1}}\frac{k_{1}}{k_{2}}s_{1},\quad s^{*}_{2}={c_{s}^{2}}\frac{k_{3}}{k_{4}}s_{2},$$

where $c_{s}^{1,2}$ are positive constants that will be determined later. Equation 22, after substituting the scaled variables, becomes the following:

$$\frac{dn^{*}}{dt^{*}}\,=\,\frac{d_{2}k_{0}k_{4}}{d_{1}^{2}c_{s}^{2}k_{3}}s_{2}^{*}+\nu\frac{k_{4} e_{0}d_{2}}{d_{1}^{2}}\frac{s_{1}^{*}s_{2}^{*}}{c_{s}^{1}c_{s}^{2}\,+\,c_{s}^{2}s_{1}^{*}\,+\,s_{1}^{*}s_{2}^{*}}- n^{*}(1+n^{*}). $$

(23)

For simplicity, we drop the ^∗, but hereafter the variable n, s ₁, s ₂ and t refer to the non-dimensionalized version of the same variables. We obtain the non-dimensional version of Eq. 22:

$$ \frac{dn}{dt}=A\frac{s_{2}}{c_{s}^{2}}+B\frac{s_{1}s_{2}}{c_{s}^{1}c_{s}^{2}+c^{2}_{s}s_{1}+s_{1}s_{2}}-n-n^{2}, $$

(24)

where

$$ A=\frac{d_{2}k_{0}k_{4}}{d_{1}^{2}k_{3}}, $$

(25)

and

$$ B=\frac{\nu k_{4}e_{0}d_{2}}{d_{1}^{2}}. $$

(26)

The ratio d ₂/d ₁ represents the relative contributions of density- and non-density-related factors to mortality, which are not easy to separate and can vary depending on the habitat type and species (Harper 1977). In disturbance-prone habitats such as coastal systems, for example, mortality is mainly due to erosion and is less affected by the density of plants (e.g. Keddy 1981). However, in tropical habitats and forests, where competition could become intensive, density can be a driving force of mortality (e.g. Lambers et al. 2002). Here, we assumed that mortality is mainly caused by density-independent regulators (as in Watkinson and Harper1978; Keddy1981; Watkinson et al. 1989), though the opposed scenario can also be used without the loss of the generality of results. In the annual species Vulpia fasciculata plant, for example, mortality due to density-independent factors ranged from 7 to 40 % in the study of Watkinson and Harper (1978). Similarly, fire, as a density-independent factor, killed ∼40 % of flowering Sorghum intrans plants (Watkinson et al. 1989).

If we assume 30 % of plants die due to factors unrelated to the density while competition, i.e. density-dependent processes between plants, accounts for 9 % death over the growing season of 100 days, we will obtain a $d_{2}/{d_{1}^{2}} = 100$ day per unit density. We also assume that pollen deposition and the breakdown of the complex C ₂ happen at the same rate, so that k ₃ = k ₄. This gives k ₃/k ₄=1 per unit density. The selection of these values is somewhat arbitrary and for mathematical convenience, e.g. parameter A, which is dimensionless, is now reduced to A = k ₀×100 day where 100 is the length of growing season assumed for an annual plant: a longer or shorter growing period will not qualitatively affect the results. Parameter A now represents the selfing rate (k ₀) over the season. Nevertheless, the above mortality values are still realistic and fall within the range of reported values. Using the above argument, one also obtains that $B=\nu \cdot k_{4} \cdot e_{0}\cdot 100 \: \mathrm { day} \cdot \mathrm {m}^{2}.$ The parameter e ₀ represents the density of pollinators at the time of the introduction and was assumed to be 1 m⁻², which lies within the range of pollinator abundance in several surveys (e.g. Grixti and Packer 2006; Janovsky et al. 2013; Moreira et al. 2015; Westphal et al. 2008). The density of pollinators, e ₀, can change (positively) with the plant density (Kunin 1993), but as we are mainly interested in the potential effect of hybridization on the quality of pollination, we assumed that this parameter is constant and independent of plant density. Now, the dimensionless parameter B becomes B = ν⋅k ₄⋅100 day,which represents the fertile outcrossing rate as function of cross pollination rate, k ₄, and fertilization success rate, ν, in a growing season.

Equation 24 is thus a dimensionless equation that describes the rate of change in n as a function of A and B, representing the selfing and outcrossing rates over the season, respectively. This equation was then used to derive two different models: a model with a single invading species and the other with two hybridizing species.

Equation 22 and its dimensionless version Eq. 24 describe the multiple temporal processes within the life cycle of a plant simultaneously as a single phenomenon where n can be any expression of a plant’s developmental stages from seed to seedling to mature plant. We assume that, on average, the ovule density and pollen load density are proportional to the plant density, which, in turn, is proportional to the matured seed density. We choose the constants of proportionality to be the undetermined constants from the non-dimensionlization, ${c_{s}^{1}}$ and ${c_{s}^{2}}$, so that $s_{1}={c_{s}^{1}} n$ and $s_{2}={c_{s}^{2}} n$.

Appendix B: Model of two non-hybridizing genotypes: competition model

B.1 The model

In this section, we expand our single species model (4) to a system of two competing genotypes (species), namely the native (resident) genotype, denoted by N, and a new incoming (invading) genotype, denoted by I. The two genotypes interact with each other through density-dependent regulation, e.g. death rate is a function of the total population size N + I, but they do not intercross, i.e no hybridization occurs between them. The two genotypes are assumed to vary in their selfing, parameter A, and outcrossing, parameter B, rates. Extending the Eq. 4 to include two genotypes results in the system of two ODEs:

$$ \left\{\begin{array}{lllll} n_{I}^{\prime}(t)=n_{I}(A_{I}+B_{I}\frac{n_{I}}{1+n_{I}+n_{I}^{2}}-(1+n_{I}+n_{N}))\\ n_{N}^{\prime}(t)\,=\,n_{N}(A_{N}+B_{N}\frac{n_{N}}{1+n_{N}+n_{N}^{2}}-(1+n_{N}+n_{I})),\!\!\! \end{array}\right. $$

(27)

where n _I and n _N denote the densities of invading and native genotypes, respectively. From Eq. 27, it can be seen that the existence of another genotype will always reduce the growth of either genotype. Now, suppose that the native genotype is at equilibrium (i.e. $n_{N}\equiv \bar {n}_{N}$) at the time the propagules of the new genotype arrive at the location. As the density of the native is kept constant, the dynamics of invading can be studied through the following:

$$n_{I}^{\prime}(t)=n_{I}\left( A_{I}+B_{I}\frac{n_{I}}{1+n_{I}+n_{I}^{2}}-(1+n_{I}+\bar{n}_{N})\right)=\bar{F}(n_{I}). $$

(28)

As we assumed that the native genotype always remains at equilibrium and that it cannot interbreed with the invading genotype, the differential model of the invader in the presence of a resident species is similar to that of the single invading genotype, with the only exception being density-dependent mortality formulation.

B.2 Analysis of the model

Similar to the single species dynamics (Appendix C), bifurcation analysis of the above model gave rise to four outcomes, contingent on the value of A _I with respect to the bifurcation point $(1+\bar {n}_{N})$ and the value of B _I with respect to 1 (Fig. 7). When $A_{I}>(1+\bar {n}_{N}),$ population grows either with a weak Allee effect (if B _I >1) or with no Allee effect (if B _I <1). In the case of $A_{I} <(1+\bar {n}_{N}),$ the invading population will have no growth and go extinct (if B _I is small) or experience strong Allee effects (if B _I becomes large).

Fig. 7

Bifurcation diagram (27) depicting four dynamic outcomes for the invading genotype as related to parameters A (representing selfing rate) and B (representing outcrossing rate). See the “B.2 Analysis of the model” section for more details

Appendix C: Analysis of the dynamics of invasion in the single species model

C.1 Single species model

In this section, we derive the analysis performed on the function F of Eq. 4, in order to characterize its shape and the type of growth for n.

When A > 1, using Descartes’ rule of sign (Appendix C.2 or in Pearson (1990), chapter 1), one can conclude that there exists a unique positive root n ^∗ for F. Given that when n > 0 is small, F(n)>0 whereas $F(n)\to -\infty$ as $n\to +\infty$ we obtain

$$F(n)>0 \text{ if } 0<n<n^{*} \text{ and } F(n)<0 \text{ if } n>n^{*}.$$

We can conclude that for A>1, the growth function F is always positive and there exists a stable equilibrium n ^∗>0. Moreover, the per capita growth rate

$$\frac{F(n)}{n}=A-1-n+\frac{Bn}{1+n+n^{2}}$$

is maximal at 0 if and only if B≤1. The model will thus give rise to a weak Allee effect if and only if B>1 (Fig. 1).

When A=1, the model reduces to

$$ F(n)=n^{2}\left( \frac{B}{1+n+n^{2}}-1\right), $$

(29)

which is positive for some values of n>0 if and only if B>1. In this case, there exists $n^{*}_{B}>0$ such that

$$F(n)>0, \quad \text{ if}~0<n<n^{*}_{B} \text { and}\;F(n)<0 \text{ if}\; n>n^{*}_{B}.$$

Lastly, if we assume that A<1, the use of Descartes’ rule of sign affirms the growth function is always negative if B≤2−A. To investigate the model behaviour for situations where the population may grow, i.e. B>2−A, we use Cardan’s method (Appendix C.3; also see Pearson (1990), chapter 1). We found that the function F can either be bistable or negative depending on the sign of the so-called discriminant Δ (see Appendix C.3). Accordingly, when Δ<0, there exists two positive roots $0<n^{*}_{1}<n^{*}_{2}$ such that

$$F(0)=F(n^{*}_{1})=F(n^{*}_{2})=0,$$

$$F(n)<0,\quad \text{if}\;0<n<n^{*}_{1} \text{ or } n>n^{*}_{2}$$

and

$$F(n)>0,\quad\: \text{if}\;n^{*}_{1}<n<n^{*}_{2}.$$

Whereas 0 and $n^{*}_{2}$ equilibria are stable, the $n^{*}_{1}$ equilibrium is unstable. The function F is thus bistable when Δ<0 and there is a strong Allee effect. When Δ≥0, F(n)≤0 for all n≥0 and the only non-negative stable equilibrium is 0, implying that there is no growth.

C.2 Descartes’ rule of sign

Descartes’ rule of sign is used to obtain an upper bound on the number of positive roots of a polynomial (Pearson 1990, chapter 1). That is, if the polynomial is ordered by descending variable exponent, then the number of sign changes between consecutive non-zero coefficients is greater or equal to the number of positive roots of the polynomials. Moreover, if the number of positive roots is not equal to the number of sign changes of non-zero consecutive coefficient, then it is smaller than it by an even number. Let Z denotes the number of sign changes of consecutive non-zero coefficients and P the number of positive roots, then

$$P\geq 0 \text{ and } P=Z-2k, \text{ for some } k\in\mathbb{N}.$$

This means for example that if Z is odd then there is at least one positive root and if Z=1, then we know that there is exactly one positive root.

This rule can be extended to negative roots also by considering −x in the polynomial instead of x and again using the above rule.

In our framework we write F as follows:

$$ F(x)\,=\,\frac{x(-x^{3}\,+\,(A-2)x^{2}\,+\,(A\,-\,2\,+\,B)\,+\,(A-1))}{1+x+x^{2}}. $$

(30)

As the denominator in (30) is positive when x is non-negative, using Descartes’ rule of signs for the third order polynomial in the numerator, we found that if A>1 then there is always only one sign change between the coefficients and thus there is exactly one positive root x ^∗ > 0.

C.3 Cardan’s method

A brief description of Cardan’s method can be found in Pearson (1990), chapter 1. Here, we provide more details on the method. We are interested in solving the general cubic equation, using Cardan’s method

$$x^{3}+ax^{2}+bx+c=0. $$

(31)

We first reduce this equation in a degenerate cubic equation of the form

$$x^{3}+px+q=0. $$

(32)

This can be done by the following substitution:

$$x=t-\frac{a}{3}.$$

Then Eq. 31 becomes

$$ t^{3}+pt+q=0, $$

(33)

with p=(3b−a ²)/3 and q=(2a ²−9a b+27c)/27.

• If p = q=0, then the only solution is t=0, which is equivalent to x=−a/3,

• If q ≠ 0 and p=0, then $t=\{\sqrt [3]{-q}, \sqrt [3]{-q}(-\frac {1}{2}-\sqrt {3}\frac {i}{2}),\sqrt [3]{-q}(-\frac {1}{2}+\sqrt {3}\frac {i}{2})\}$, and x = t−a/3,

• If p ≠ 0 and q=0, then $t=\{0, -\sqrt {-p}, \sqrt {-p}\}$ and x = t−a/3,

• If p ≠ 0 and q ≠ 0, then we introduce u and v such that t = u−v. Given that

  $$ (u-v)^{3}+3uv(u-v)-(u^{3}-v^{3})=0 $$

  (34)

  we obtain p=3u v and q=−(u ³−v ³). Note that u ≠ 0 and v ≠ 0. This implies that $v=\frac {p}{3u}$ and

  $$u^{3}-\left( \frac{p}{3u}\right)^{3}+q=0\Leftrightarrow u^{6}+qu^{3}-\frac{p^{3}}{27}=0,$$

  thus, u ³ solves a quadratic equation and

  $$u^{3}=\frac{-q\pm\sqrt{\Delta}}{2},$$

  with Δ = q ²+4p ³/27. As v ³−u ³ = q, it follows that

  $$v^{3}= \frac{q\pm\sqrt{\Delta}}{2}.$$

  Notice that as x = t−a/3 = u−v−a/3, the same results would be obtain if we change $+\sqrt {\Delta }$ to $-\sqrt {\Delta }$, we can thus consider, without loss of generality, only the case with $+\sqrt {\Delta }$ and

  $$u^{3}=-\frac{q+\sqrt{\Delta}}{2}, \: v^{3}=\frac{q+\sqrt{\Delta}}{2}.$$

  Depending on the sign of Δ, now we can reach conclusion about the nature of the roots x:

  • If Δ=0, then $u=\{\sqrt [3]{-\frac {q}{2}},\sqrt [3]{-\frac {q}{2}}(-\frac {1}{2}-i\frac {\sqrt {3}}{2}),\sqrt [3]{-\frac {q}{2}}(-\frac {1}{2}+i\frac {\sqrt {3}}{2})\}$ and using the fact that $3uv=p\in \mathbb {R,}$ we know that u and v must be conjugate. Substituting for x, we get

    $$x=\left\{\sqrt[3]{\frac{q}{2}}-\frac{a}{3}, 2\sqrt[3]{-\frac{q}{2}}-\frac{a}{3}\right\}.$$

    and the roots are real with $\left (2\sqrt [3]{-\frac {q}{2}}-\frac {a}{3}\right )$ being a simple root and $(\sqrt [3]{\frac {q}{2}}-\frac {a}{3})$ being a double root.

  • If Δ>0, then using the same arguments, we obtain

    $$u^{3}=\frac{-q+\sqrt{\Delta}}{2}\in\mathbb{R}, \: v^{3}=\frac{q+\sqrt{\Delta}}{2}\in\mathbb{R},$$

    and

    $$u_{1}= \sqrt[3]{\frac{-q+\sqrt{\Delta}}{2}} \in\mathbb{R},\: v_{1}=\sqrt[3]{\frac{q+\sqrt{\Delta}}{2}}\in\mathbb{R},$$
    $$u_{2}=\sqrt[3]{\frac{-q+\sqrt{\Delta}}{2}}\left( -\frac{1}{2}+i\frac{\sqrt{3}}{2}\right),\: v_{2}=\sqrt[3]{\frac{q+\sqrt{\Delta}}{2}}\left( -\frac{1}{2}-i\frac{\sqrt{3}}{2}\right),$$
    $$u_{3}=\sqrt[3]{\frac{-q+\sqrt{\Delta}}{2}}\left( -\frac{1}{2}-i\frac{\sqrt{3}}{2}\right),\: v_{3}=\sqrt[3]{\frac{q+\sqrt{\Delta}}{2}}\left( -\frac{1}{2}+i\frac{\sqrt{3}}{2}\right).$$

    We can then deduce $x_{1}\in \mathbb {R}$, $x_{2},\:x_{3}\in \mathbb {C}$.

  • If Δ<0, then

    $$u^{3}=\left( \frac{-q+i\sqrt{-{\Delta}}}{2}\right)\in\mathbb{C},\: v^{3}=\left( \frac{q+i\sqrt{-{\Delta}}}{2}\right)\in\mathbb{C}.$$

    Using the trigonometric representation of a complex number, we can write

    $$u^{3}=r(cos(\varphi)+i\sin(\varphi)),$$

    with $r^{2}=\left (-\frac {q}{2}\right )^{2}-{\Delta }=\left (-\frac {p}{3}\right )^{3}$ and $\cos (\varphi )=-\frac {q}{2\sqrt {(-\frac {p}{3})^{3}}}$. We know that

    $$u=r^{1/3}\left( cos\left( \frac{\varphi}{3}\right)+isin\left( \frac{\varphi}{3}\right)\right)$$

    and similarly

    $$v=r^{1/3}\left( -cos+isin\left( \frac{\varphi}{3}\right)\right).$$

    Then, $x=\{2r^{1/3}cos(\frac {\varphi }{3})-\frac {a}{3},2r^{1/3}cos(\frac {\varphi +2\pi }{3})-\frac {a}{3},2r^{1/3}cos(\frac {\varphi +4\pi }{3})-\frac {a}{3}\}$ and the roots of Eq. 31 are distinct and real.

Now, we apply this derivation to our equations:

$$ F(n)\,=\,\frac{-n(n^{3}\,-\,(A\,-\,2)n^{2}\,-\,(A\,+\,B\,-\,2)n\,-\,(A\,-\,1))}{1+n+n^{2}}. $$

(35)

We want to find the roots (or at least their sign) of the cubic polynomial on the numerator:

$$ P(n)=n^{3}-(A-2)n^{2}-(A+B-2)n-(A-1), $$

(36)

when A<1 and B>2−A. First, using the Descartes’ rule of signs, we know that this polynomial has either two or no positive roots. Moreover, we know that $P(n)\to +\infty$ (respectively $-\infty$) when $n\to +\infty$ (respectively $-\infty$); P(0)=−(A−1)>0, $P^{\prime }(0)=-(A+B-2)<0$. Now we apply Cardan’s method. Using the same notation given in the derivation above, we have

$$a\,=\,-(A-2)\!>\!0,\: b\,=\,-(A+B-2)\!<\!0,\: c\,=\,-(A-1)\!>\!0$$

and thus

$$q=\frac{2a^{2}-9b+27c}{27}>0, \quad p=\frac{3b-a^{2}}{3}<0.$$

Denoting Δ as the determinant defined in the derivation above, we have:

• If Δ>0, then there is only one real root and thus it cannot be positive. Moreover, as P(0)>0 and P(n) do not change sign for n>0, we conclude that P(n)>0 for all n>0,

• If Δ=0, there are only two real roots and the simple one is $2\sqrt [3]{-\frac {q}{2}}-\frac {a}{3}<0$. Denoting $n^{*}_{2}$ as the second real root, we have $P\left (n^{*}_{2}\right )=P^{\prime }\left (n^{*}_{2}\right )=0$ and P only changes sign at the simple root $2\sqrt [3]{-\frac {q}{2}}-\frac {a}{3}<0$. As P(0)>0, we have P(n)≥0 when n>0,

• If Δ<0, we know that the polynomial has three real roots $n_{0}^{*},\:n_{1}^{*}$ and $n_{2}^{*}$ such that

  $$ P(n)=(n-n_{0}^{*})(n-n^{*}_{1})(n-n^{*}_{2}), $$

  (37)

  and knowing that P(0)>0 and $P^{\prime }(0)<0$; this implies that $n^{*}_{0}<0<n^{*}_{1}< n^{*}_{2}$, P(n)>0 when $n\in (0,n^{*}_{1})$, P(n)<0 when $n\in \left (n^{*}_{1},n^{*}_{2}\right )$ and P(n)>0 when $n>n^{*}_{2}$.

As $F(n)=\frac {-nP(n)}{1+n+n^{2}}$, we obtain the required conclusion about the shape of F.

Appendix D: Analysis of the dynamics of hybrid and invading genotypes for the hybridization scenario

D.1 Stability analysis of trivial equilibria in the hybridization model (11)–(12)

In this section, we detail the analysis of the two-dimensional ODE system (11)–(12) and more precisely the stability analysis of the trivial equilibrium (0,0). We first need to compute the Jacobian of the system at the trivial equilibrium:

$$ J(0,0)=J^0=\left( \begin{array}{ccc} J^0_{11} & J^0_{12}\\ J^0_{21} &J^0_{22}\end{array}\right) $$

(38)

where

$$J^{0}_{11}\,=\,A_{I} -(1+\bar{n}_{N}),\!\quad \!J^{0}_{12}\,=\,\frac{A_{H}}{4},\!\quad \!J^{0}_{21}\,=\,\frac{B_{IN}\bar{n}_{N}}{1+\bar{n}_{N}}+B_{IN}\bar{n}_{N}$$

and

$$J^{0}_{22}=\frac{A_{H}}{2}+\frac{B_{HN}\bar{n}_{N}}{2}+\frac{B_{HN}\bar{n}_{N}}{2(1+\bar{n}_{N})}-(1+\bar{n}_{N}).$$

The stability properties of (0,0) can be inferred from the sign of the trace and the determinant of the Jacobian at (0,0):

$$ tr(J^0)=J^0_{11}+J^0_{22},\quad det(J^0)=J^0_{11} J^0_{22}-J^0_{12} J^0_{21}. $$

(39)

Indeed, the trivial equilibrium is stable if t r(J ⁰)<0 and d e t(J ⁰)>0. It is an unstable node or unstable spiral if d e t(J ⁰)>0 and t r(J ⁰)>0 and an unstable saddle if d e t(J ⁰)<0. We know from the positivity of the parameters that $J^{0}_{12},J^{0}_{21}>0$, thus

$$ \left( tr(J^0)\right)^{2}-4det(J^0)=(J^0_{11}-J^0_{22})^{2}+4J_{12}^0J^0_{21}>0 $$

(40)

and the eigenvalues of (0,0)

$$\lambda_{1,2}^0=\frac{tr(J^0)\pm\sqrt{[tr(J^0)]^{2}-4det(J^0)}}{2}\in\mathbb{R} $$

(41)

Both eigenvalues are real, which means that the trivial equilibrium is either a node (stable or unstable) or a saddle.

Note that the trace and the determinant of the Jacobian at (0,0) only depend on B _{I N} and B _{H N} , given that A _H = A _I = A were chosen to be fixed. Our goal is to determine the parameter space under which either the trace and the determinant are positive (and (0,0) will then be an unstable node) or the determinant is negative (and (0,0) will be an unstable saddle): if any of these conditions (t r(J ⁰),d e t(J ⁰)>0 or d e t(J ⁰)<0) are met, the trivial equilibrium (0,0) will become unstable. For the first case (i.e. both trace and det >0), we have the two following inequations:

$$\begin{array}{@{}rcl@{}} tr(J^0)&>&0\Leftrightarrow B_{HN}\\ &>&-\frac{2(\bar{n}_{N}+1)}{\bar{n}_{N}(2+\bar{n}_{N})}\left( \frac{3A}{2}-2(1+\bar{n}_{N})\right)\\ &=&T^0>0\end{array} $$

(42)

and

$$\begin{array}{@{}rcl@{}} det(J^0)&>&0\Leftrightarrow B_{IN}<\frac{4(1+\bar{n}_{N})(A-(1+\bar{n}_{N}))}{A\bar{n}_{N}(2+\bar{n}_{N})}\\ &&\times\left[\frac{A}{2}-(1+\bar{n}_{N})+\frac{\bar{n}_{N}(2+\bar{n}_{N})}{2(1+\bar{n}_{N})}B_{HN}\right].\end{array} $$

(43)

But, when d e t(J ⁰)>0 then t r(J ⁰)<0. A d e t(J ⁰)>0 implies that

$$B_{HN}<\left( (1+\bar{n}_{N})-\frac{A}{2}\right)\frac{2(1+\bar{n}_{N})}{\bar{n}_{N}(2+\bar{n}_{N})}<T^{0},$$

as we assumed that A−1<0 in Eq. 14, to meet a stable trivial equilibrium in Eq. 4. Thus, when hybridization is allowed, the trivial equilibrium (0,0) is either a stable node (d e t(J ⁰)>0 and t r(J ⁰)<0) or an unstable saddle (d e t(J ⁰)<0). The trivial equilibrium (0,0) is unstable and thus facilitation is warranted if, and only if, the determinant becomes negative, i.e when B _{I N} and B _{H N} are such that

$$\begin{array}{@{}rcl@{}} B_{IN}\!\!&>&\!\!\frac{4(1+\bar{n}_{N})(A-(1+\bar{n}_{N}))}{~}\\ &&\!\!\times{A\bar{n}_{N}\!(2\!\,+\,\bar{n}_{N})}\!\left[\!\frac{A}{2}\,-\,(1\!\,+\,\bar{n}_{N})\,+\,\frac{\bar{n}_{N}(2\!\!+ \!\bar{n}_{N})}{2(1\,+\,\bar{n}_{N})}B_{HN}\!\!\right]\!. \end{array} $$

(44)

Changing the values of the parameters B _{I N} and/or B _{H N} to move from a negative determinant to a positive determinant thus generates a saddle node bifurcation at (0,0).

D.2 Analysis of the null clines of the hybridization model (11)–(12)

In this section, we examine the mechanisms underlying the saddle node bifurcation of (0,0) by computing the null clines of Eqs. 11–12, i.e. the points (n _I ,n _H ) in the phase plane (n _I ,n _H ) such that

$$ n_{I}^{\prime}\equiv 0 \text{ and } n_{H}^{\prime}\equiv0. $$

(45)

Note that the intersection of these null clines gives the equilibrium point(s). We computed these null clines for n _I ,n _H >0 (Fig. 8) with parameters A and B in Eq. 13 chosen so that 0 is stable in the model with a single species (here A=0.8, B=3). As can be inferred from Eq. 11, as well as shown in Fig. 8, the null cline for n _I , i.e the curve such that $n_{I}^{\prime }\equiv 0$, neither depends on B _{H N} (Fig. 8a) nor on B _{I N} (Fig. 8b). Therefore, if we change the values of these two parameters to change the stability of the trivial steady state (0,0) of Eqs. 11–12, only the null cline associated with the hybrid population n _H changes. When B _{I N} and B _{H N} are very small (e.g. =0.5), there is no non-trivial equilibrium, i.e. the null clines of the two genotypes only intersect at the stable equilibrium (0,0). With this parameter setting, the invading genotype will go extinct even if it can hybridize with the native genotype. A small increase in B _{I N} or B _{H N} (e.g. =1) gives rise to three equilibrium points: (0,0), $n^{1}=({n_{I}^{1}},{n_{H}^{1}})$ and $n^{2}=({n_{I}^{2}},{n^{2}_{H}})$, where the null clines of the two genotypes intersect (Fig. 8). In this case, the trivial equilibrium is still stable (Fig. 3) but the stability of the two non-trivial equilibria requires further analysis (see below). Here, a low density of the invading genotype will again decline to 0. When large values of B _{I N} and B _{H N} were chosen (e.g. =2), the model resulted in one non-trivial equilibrium n ² only (the intersection of solid black line and the dotted grey line in Fig. 8). In this case, (0,0) is unstable (Fig. 3) but n ² is stable and the densities of invading and hybrid genotypes will converge to the non-trivial equilibrium n ² (Appendix D.3). To understand what happens if the density of invading and hybrid genotype does not converge to the trivial equilibrium (0,0), we need to determine the stability of the non-trivial equilibria. To do so, we analysed the vector field of Eqs. 11–12 in more detail, as shown in Appendix D.3. The transition from a solution with two non-trivial equilibria to a solution with a single non-trivial equilibrium (Fig. 8) corresponds to a saddle node bifurcation at (0,0) (Fig. 3). As the parameters B _{I N} and B _{H N} increase, the intermediate equilibrium n ¹ approaches the trivial equilibrium (0,0) (Fig. 8), leading to a saddle node bifurcation. When B _{I N} and B _{H N} were large enough so that d e t(J ⁰)<0, only two equilibria result, the trivial (0,0) which is now a saddle and the non-trivial n ², which is stable.

Fig. 8

Null clines for invading, genotype, $n_{I}^{\prime }\equiv 0$ (solid black line) and hybrid genotype $n_{H}^{\prime }\equiv 0$ (grey lines) for different B _{H N} values (a) or B _{I N} values (b) of 0.5 (solid grey line), 1.5 (dashed grey line), 2 (for B _{H N} ) or 3 (for B _{I N} ) (dotted grey line). The growth rate of the invading genotype is positive in the areas above its null cline ($n^{\prime }_{I}\equiv 0$: solid black curve) but negative below it. For the hybrid genotype, positive growth rate occurs in the areas under the null clines ($n^{\prime }_{H}\equiv 0$: grey lines) while it is negative above them. The solid circles show the intersections of the two null clines (black and grey curves) which correspond to different equilibria. Other parameters were B _{I N} =1 (a) or B _{H N} =1 (b), A=0.8 and B=3

D.3 Stability analysis of non-trivial equilibria in the hybridization model (11)–(12)

To investigate the stability of the non-trivial equilibrium point(s) of Eqs. 11–12, we analyse the vector field around the equilibria, as outlined in Murray (2002), chapter 7.3. We only provide the details for one scenario, the intermediate equilibrium n ¹ (dotted line in Fig. 8), as the procedure is the same, to a large extent, for other non-zero equilibria. At this point n ¹:

$$ \frac{dn_{H}}{dn_{I}}|_{F_{I}\equiv0}>\frac{dn_{H}}{dn_{I}}|_{F_{H}\equiv0}>0, $$

(46)

where F _I and F _H are the “growth” functions for the invading and hybrid genotypes in Eqs. 11–12, that is Eqs. 11 and 12 can be written in term of F _I and F _H :

$$\left\{\begin{array}{llllll} n_{I}^{\prime}(t)=F_{I}(n_{I},n_{H})\\ n_{H}^{\prime}(t)=F_{H}(n_{I},n_{H}) \end{array}\right. $$

(47)

Given that F _I is positive above and negative below its null cline, $n^{\prime }_{I}\equiv 0$, (the black curve in Fig. 8) whereas F _H is negative above and positive below its null cline, $n^{\prime }_{H}\equiv 0$, (the grey curve in Fig. 8) we can conclude that at the very vicinity of n ₁,

$$ \frac{\partial F_{I}}{\partial n_{I}}|_{F_{I}\equiv0}<0, \: \frac{\partial F_{I}}{\partial n_{H}}|_{F_{I}\equiv0}>0 $$

(48)

and

$$ \frac{\partial F_{H}}{\partial n_{I}}|_{F_{H}\equiv0}>0, \: \frac{\partial F_{H}}{\partial n_{H}}|_{F_{H}\equiv0}<0. $$

(49)

This implies that $det(J({n^{1}_{I}},{n^{1}_{H}}))<0$ so that the non-trivial equilibrium n ¹ is a saddle. Using the same analysis for the other non-trivial equilibrium n ² yields

$$ \frac{dn_{H}}{dn_{I}}|_{F_{I}\equiv0}<\frac{dn_{H}}{dn_{I}}|_{F_{H}\equiv0}, $$

(50)

and we can conclude that n ² is stable. We also conclude from Fig. 8 that there exists a positively invariant set containing the equilibria and observe that there exists no limit cycle in this invariant domain. Using Poincaré-Bendixon theory (as in Murray (2002), chapter 7.3), we conclude that the solutions n _I and n _H of Eqs. 11–12 converge to an equilibrium which can only be (0 , 0 ) when it is stable or n ² (when it exists).