Local-regional richness relationships and alternative stable states in metacommunities with local facilitation

Abstract

Local-regional species richness relationships have been used to infer relative contributions of local and regional forces to determining the richness of local communities. Although most previous research assumed competition as major local species interactions, growing empirical evidence suggests that facilitation is also an important driver of local community dynamics. Here, I explore how facilitation affects the shape of local-regional richness relationships, by incorporating local facilitation into a patch-occupancy model of metacommunity dynamics. I find that facilitation can generate local-regional richness relationships with the alternative stable states of mean local richness at intermediate to high levels of regional richness. These alternative stable states tend to occur in a metacommunity in moderately harsh environments. This result cautions against assuming that only competition can be primarily important local interactions when interpreting the shapes of local-regional richness relationships. Moreover, the possibility of alternative stable states suggests that gradual decline of regional species diversity might cause a sudden collapse of metacommunities with local facilitation.

Appendices

Appendix A

Approximate equations for the mean and variance of local richness

I derive two equations describing the dynamics of the mean and variance of local richness, $L = \sum\nolimits_{m = 0}^R {{{m}}{{{p}}_m}}$ and $V = \sum\nolimits_{m = 0}^R {{{\left( {m - L} \right)}^2}{p_m}}$, respectively, from the full equations (8a). To obtain the equation for L, multiply both sides of all R + 1 equations of (8a, b, c) by m. Then adding the left- and right-hand sides of these equations yields

$$ \frac{{dL}}{{{dt}}} = h\left( {\sum\limits_{m = 0}^R {{{m}}{{{g}}_m}{p_m}} } \right)\sum\limits_{m = 0}^R {\left( {1 - \frac{m}{R}} \right){s_m}{p_m}} - \sum\limits_{m = 0}^R {{{m}}{{{e}}_m}{p_m}} . $$

(A1a)

The dynamic equation for V is similarly obtained by multiplying both sides of all R + 1 equations in (8a–c) by (m–L)², and adding all terms on each side across all equations.

$$ \frac{{dV}}{{{dt}}} = h\left( {\sum\limits_{m = 0}^R {{{m}}{{{g}}_m}{p_m}} } \right)\sum\limits_{m = 0}^R {\left\{ {2(m - L) + 1} \right\}\left( {1 - \frac{m}{R}} \right){s_m}{p_m}} - \sum\limits_{m = 0}^R {\left\{ {2(m - L) - 1} \right\}{{m}}{{{e}}_m}{p_m}} . $$

(A1b)

These equations are not closed in terms of L and V. To obtain closed forms, the functions g _m , S _m , and e _m are Taylor expanded with respect to m around L, and plugged into equations (A1a and b). Neglecting the third- or higher-order terms, closed equations are

$$ \frac{{dL}}{{{dt}}} = h\left\{ {{{L}}{{{g}}_L} + V\left( {g_L^\prime + \frac{1}{2}{{Lg}}_L^{\prime \prime }} \right)} \right\}\left( {{s_L} - \frac{{{{L}}{{{s}}_L}}}{R} + V\frac{{ - 2s_L^\prime - {{Ls}}_L^{\prime \prime } + {{Rs}}_L^{\prime \prime }}}{{2R}}} \right) - \left\{ {{{L}}{{{e}}_L} + V\left( {e_L^\prime + \frac{1}{2}{{Le}}_L^{\prime \prime }} \right)} \right\}, $$

(A2a)

$$ \begin{array}{*{20}{c}} {\frac{{dV}}{{{dt}}} = h\left\{ {{{L}}{{{g}}_L} + V\left( {g_L^\prime + \frac{1}{2}Lg_L^{\prime \prime }} \right)} \right\}\left( {{s_L} - \frac{{{{L}}{{{s}}_L}}}{R} + V\frac{{ - 4{s_L} - 2s_L^\prime - 4{{Ls}}_L^\prime + 4{{Rs}}_L^\prime - {{Ls}}_L^{\prime \prime } + {{Rs}}_L^{\prime \prime }}}{{2R}}} \right)} \\{ - \left\{ { - {{L}}{{{e}}_L} + V\left( {2{e_L} - e_L^\prime + 2{{Le}}_L^\prime - \frac{1}{2}{{Le}}_L^{\prime \prime }} \right)} \right\}} \\\end{array}, $$

(A2b)

where $g_L^\prime$, $g_L^{\prime \prime }$, $s_L^\prime$, $s_L^{\prime \prime }$, $e_L^\prime$, and $e_L^{\prime \prime }$ are the first and second derivatives of g _m , s _m , and e _m , respectively, with respect to m evaluated at L. Once the functions g _m , s _m , and e _m are specified, equations (A2) can be used to find approximate values of equilibrium local richness and variance. I confirmed that the approximate equilibrium values show good correspondence to the equilibrium values based on equations (7b).

Moreover, these approximate equations can be used to study analytically the changes of mean and variance of local richness along the gradient of regional richness. These results provide an analytical proof for the existence of alternative stable states in facilitative metacommunities. There are at least five cases (I-V) that are analytically tractable. Case I assumes neutral local interactions, in which I set g _m  = G, S _m  = S, and e _m  = E (when α = β = 0 and S = T). Case II assumes competitive local interactions, in which the immigration rate is inversely related to local richness (g _m  = G/m, s _m  = S, and e _m  = E, when α = −1, S = T, and β = 0). Case III also assumes competitive interactions, in which the extinction rate is linearly related to local richness (g _m  = G, s _m  = S, and e _m  = Em, when α = 0, S = T, and β = 1). Case IV assumes facilitative local interactions, in which the immigration rate is a linear function of local richness (g _m  = Gm, s _m  = S, and e _m  = E, when α = 1, S = T, and β = 0). Case V also assumes facilitative interactions, in which the extinction rate is inversely related to local richness (g _m  = G, s _m  = S, and e _m  = E/m, when α = 0, S = T, and β = −1). Approximate equations and the equilibrium values of L and V for cases (I-V) are summarized in Table 1.

Appendix B

A simulation model incorporating species difference in facilitative performance

The patch-occupancy model in the main text assumes no difference in species. Here I develop another simulation model to confirm that alternative stable states can occur in metacommunities with local facilitation even when species differ in their facilitative performance. I model a metacommunity consisting of a number of local communities linked by dispersal. For simplicity, I assume that there is no difference in environmental conditions among local communities. I formulate a discrete-time system, in which the following sequence of events occurs within a unit time step: (1) within-community local interactions, (2) local extinction due to low densities, (3) local extinction due to local disturbance, and (4) between-community dispersal.

Denote z _ij (t) as the density of the j-th species in the i-th local community at time t. Local interactions change the density z _ij (t) to w _ij (t) according to the following equation.

$$ {w_{ij}}(t) = {z_{ij}}(t)\exp \left[ {r - \frac{m}{{1 + \sum\limits_{j = 1}^R {{a_{jk}}{z_{ik}}(t)} }} - \sum\limits_{j = 1}^R {{b_{jk}}{z_{ik}}(t)} } \right], $$

(B1)

where r is the basic birth rate and m is the basic death rate. a _jk represents the per capita facilitative effect of the k-th species on the j-th species (a _jk  ≥ 0 and a _jj  = 0). Different values of a _jk confer that different species have different facilitative effects on one another. The term $m/\left( {1 + \sum\limits_{j = 1}^R {{a_{jk}}{z_{ik}}(t)} } \right)$ means that the death rate decreases as the richness and/or densities of facilitative species becomes greater, and/or as the per capita facilitative effects become stronger. b _jk (>0) represents the per capita competitive effect of the k-th species on the j-th species (when j ≠ k) and self-limitation (when j = k).

If local species densities fall below a certain threshold (Ω), I assume that such species become extinct from local communities. I denote as x _ij (t) the species densities after this extinction due to low densities. Further, I assume that local disturbance, such as physical destruction, physiological stress, predation, or disease outbreak, can occur at a given rate (∆) randomly on every species in every community. When this local disturbance strikes a local population, it drives the population into extinction. I denote the species densities after this disturbance-driven extinction as y _ij (t).

Finally, a fixed fraction (Ψ) of the member of each local population disperse. The total number of dispersers of the j-th species is $a\sum\limits_{i = 1}^M {\Psi {y_{ij}}(t)}$. The probability that a local community receives dispersers are assumed as $\sum\limits_{i = 1}^M {\Psi {y_{ij}}(t)} /\left( {H + \sum\limits_{i = 1}^M {\Psi {y_{ij}}(t)} } \right)$, where H is a fixed constant. The number of dispersers is equally distributed among the local communities that receive dispersers. Dispersal completes the sequence of events within a unit time-step, and yields the species densities of each local community at the next time-step.

Simulation of this model is performed for 200 time-steps, and then the mean of local richness is calculated. Survey of parameter space has found sets of parameter values with which simulations starting from different initial species densities lead to either positive or zero mean local richness after 200 time-steps. The probability of attaining positive mean local richness increased as initial local richness increased (Fig. 5). This suggests the existence of alternative stable states of mean local richness, implying that alternative stable states may be a general characteristic of facilitative metacommunities regardless of the assumption of species difference.