Differences among species in seed dispersal and conspecific neighbor effects can interact to influence coexistence

Abstract

Seed dispersal is a critical mechanism for escaping specialist natural enemies. Despite this, mean dispersal distances can vary by an order of magnitude among plant species in the same community. Here, we develop a theoretical model to explore how interspecific differences in seed dispersal alter the impact of specialist natural enemies, both on their own and though a trade-off between seed dispersal and enemy susceptibility. Our model suggests that species are more able to recover from rarity if they have high dispersal because (1) seedlings are more likely to escape their parent’s natural enemies, (2) adults are more spread out, reducing the chance that a seed will disperse near conspecifics, and (3) seedlings compete less with kin for open gaps. Differences in dispersal do not produce stabilizing mechanisms—species with low dispersal are purely at a disadvantage and do not gain a novel niche opportunity. However, dispersal-susceptibility trade-offs will be equalizing, as species disadvantaged by low dispersal will benefit from being less susceptible to specialist natural enemies. This mechanism, unlike most mechanisms of dispersal-mediated coexistence, does not require that there is an abundance of empty space: high-dispersers gain an advantage by escaping from their enemies, not by colonizing empty habitat. Our study therefore suggests that differences in dispersal are unlikely to promote diversity on their own, but may strengthen other coexistence mechanisms.

Appendices

Appendix A1: Deriving the formula for ${\tilde {\lambda }}_{j}^{\prime }(t)$

A1.1 Size assumptions for analysis

Our analysis is based on Taylor Series approximations. These work well for a parameter range where the populations are approximately quadratic, though may fail if parameters are too extreme. Thus, to bound the amount of variation in our model, we define a “small” parameter 𝜖. We write O(𝜖) to indicate that a parameter is of the same order as 𝜖 (for a technical definition, see Chesson (1994)). We assume that for all species, $Y_{j}/Y_{k} = 1 + O(\epsilon ) I_{j/_k} = 1+O(\epsilon )$, p_j is O(𝜖), and ${d_{j}^{2}}(z)/{d_{k}^{2}}(z) = 1+O(\epsilon )$, and $(d_{\text {{p}}-j}\left (z\right ))^{2}/(d_{\text {{p}}-k}\left (z\right ))^{2} = 1+O(\epsilon )$ for all z. Terms that are smaller than these, such as (p_j)², are mostly left out of our approximation, leaving the growth rate of species as ${\tilde {\lambda }}_{j}^{\prime }(t)=O(\epsilon )$. We write ≈ to indicate that numbers are identical up to their O(𝜖) terms.

In a few cases in the text, we do not approximate our model as much as we could. For example, in Table 2, we mention that the impact of seedling mortality caused by a seedling’s parentS,Δκ_Pj, contains the term $\left (p_{j} \widetilde {d_{j}d_{\text {{P}},j}} - \overline {p_{j}\widetilde {d_{j}d_{\text {{P}},j}}}\right )$. We could approximate $p_{j}\widetilde {d_{j}d_{\text {{P}},j}}$ using a Taylor series around $p_{j}=\overline {p_{j}}$ and $\widetilde {d_{j}d_{\text {{P}},j}}=\overline {\widetilde {d_{j}d_{\text {{P}},j}}}$ (where $\overline {p_{j}}$ and $\overline {\widetilde {d_{j}d_{\text {{P}},j}}}$ are the means of p_j and $\widetilde {d_{j}d_{\text {{P}},j}}$ across species) and get that

$$ p_{j} \widetilde{d_{j}d_{\text{{P}},j}} \approx \overline{p_{j}}\overline{\widetilde{d_{j}d_{\text{{P}},j}}} + (p_{j} - \overline{p_{j}})\overline{\widetilde{d_{j}d_{\text{{P}},j}}} + (\widetilde{d_{j}d_{\text{{P}},j}} - \overline{\widetilde{d_{j}d_{\text{{P}},j}}}) \overline{p_{j}} $$

(20)

$\overline {p_{j}}$ and $(\widetilde {d_{j}d_{\text {{P}},j}} - \overline {\widetilde {d_{j}d_{\text {{P}},j}}})$ are both O(𝜖); therefore, the third term is O(𝜖²). Thus, we could ignore that term, and our approximation simplifies to $p_{j} \widetilde {d_{j}d_{\text {{P}},j}} \approx p_{j} \overline {\widetilde {d_{j}d_{\text {{P}},j}}}$. However, we found that this leads to poorer approximations. Thus, there were some parts of our analysis where some O(𝜖^1.5) and O(𝜖²) terms were included implicitly. We feel that the quality of our approximations (Appendix A4) justifies their inclusion and note that previous work has suggested higher order terms will at times be necessary (Chesson 1994; Barabás et al. 2018).

A1.2 Deriving ${\tilde {\lambda }}_{j}^{\prime }(t)$

Equation 1 gives the fitness of an individual of species j. The effect of competition, C(x, t), is the number of seedlings in site x at time t. For a seed that was produced at site y to become a seedling at x, it must disperse there (which it does with probability d_j(|x − y|)), and then it must survive its natural enemies. ThuS,C(x, t) can be calculated by summing this parameter over all sites and all species. Each adults of species j will disperse Y_jd_j(|x − y|) seeds to x, and therefore, a total of $Y_{j} {\sum }_{y=1}^{\infty } N_{j}(y,t) d_{j}(|x-y|)$ seeds will be dispersed to x (where the summation is over all sites y). Each of those seeds will then become seedlings with probability S_j(x, t). Thus,

$$ C(x,t)=\sum\limits_{j=1}^{R} S_{j}(x,t) Y_{j} \sum\limits_{y=1}^{L^{2}} N_{j}(y,t) d_{j}(|x-y|) $$

(21)

where ${\sum }_{j=1}^{R}$ is the summation over all species.

For the text that follows, we write the growth rate as $\lambda _{j}^{\prime }(x,t)=(\lambda _{j}(x,t)-1)/\delta$, so that fitness is defined as 0 when there is no population growth and it is in per-generation time scale. Thus,

$$ {\lambda}_{j}^{\prime}(x,t)=Y_{j} \sum\limits_{y=1}^{L^{2}} \frac{d_{j}(|y-x|) S_{j}(y,t)}{C(y,t)} - 1. $$

(22)

Note that ${\lambda }_{j}^{\prime }(x,t)$ is defined in terms of x (the position of the adult), as it will vary across space (due to the distribution of conspecific and heterospecific adults).

The summation can be seen as a weighted average,

$$ \begin{array}{@{}rcl@{}} &&Y_{j} \sum\limits_{y=1}^{L^{2}}\frac{d_{j}(|y-x|) S_{j}(y,t)}{C(y,t)}\\ &=& L^{2} \mathbf{E}\left[{Y_{j}\frac{d_{j}(|y-x|) S_{j}(y,t)}{C(y,t)} }\right]\\ &=& L^{2} \mathbf{E}\left[{d_{j}(|y-x|)}\right] \mathbf{E}\left[{Y_{j}\frac{S_{j}(y,t)}{C(y,t)} }\right]\\ &+& L^{2} Y_{j} \text{\textbf{cov}}\left( d_{j}(|y- x|), Y_{j}\frac{S_{j}(y,t)}{C(y,t)} \right), \end{array} $$

(23)

where the covariance is over all sites y. If S_j(x, t) = I_j(1 + O(𝜖)) and Y_j/Y_k = 1 + O(𝜖) for all j and k (as we assumed), then $\frac {Y_{j}S_{j}(y,t)}{C(y,t)}\approx 1+ \ln \left \{ Y_{j} \right \} + \ln \left \{ S_{j}(x,t) \right \}-\ln \left \{ C(x,t) \right \}$. Thus,

$$ \begin{array}{@{}rcl@{}} &&Y_{j} \sum\limits_{y=1}^{L^{2}} \frac{d_{j}(|y-x|) S_{j}(y,t)}{C(y,t)} \approx\\ &&L^{2} \mathbf{E}\left[{d_{j}(|y-x|)}\right] \left( 1+ \ln \left\{ Y_{j} \right\} \right.\\ &+&\left.\mathbf{E}\left[{\ln \left\{ S_{j}(y,t) \right\}}\right]-\mathbf{E}\left[{\ln \left\{ C(y,t) \right\}} \right)\right]\\ &+& L^{2} \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\}-\ln \left\{ C(y,t) \right\} \right). \end{array} $$

(24)

Finally, note that ${\sum }_{y=1}^{L^{2}}d_{j}(|y-x|)=1$, as it is a dispersal kernel. Thus, $\mathbf {E}\left [{d_{j}(|y-x|)}\right ]=1/L^{2}$, and

$$ \begin{array}{@{}rcl@{}} &&Y_{j} \underset{y}{\sum}\frac{d_{j}(|y-x|) S_{j}(y,t)}{C(y,t)} \approx \\ &&1+ \ln \left\{ Y_{j} \right\} + \mathbf{E}\left[{\ln \left\{ S_{j}(y,t) \right\}}\right]-\mathbf{E}\left[{\ln \left\{ C(y,t) \right\}}\right] \\ &+& L^{2} \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right)-\\ &&L^{2}\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right). \end{array} $$

(25)

Putting these back into Eq. 22,

$$ \begin{array}{@{}rcl@{}} &&{\lambda}_{j}^{\prime}(x,t) \approx\\ &&\ln \left\{ Y_{j} \right\} + \mathbf{E}\left[{\ln \left\{ S_{j}(y,t) \right\}}\right]-\mathbf{E}\left[{\ln \left\{ C(y,t) \right\}}\right]\\ &+& L^{2} \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right)\\ &-&L^{2}\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right). \end{array} $$

(26)

This equation represents the fitness of an adult at site x. To calculate the finite rate of increase of a population, ${\tilde {\lambda }}_{j}^{\prime }(t)$, we average this over all sites containing adults of species j,

$$ \begin{array}{@{}rcl@{}} &&{\tilde{\lambda}}_{j}^{\prime}(t) = \mathbf{E}\left[ \left.{\lambda}_{j}^{\prime}(x,t)\right| N_{j}(x,t)=1 \right] \\ &=& \mathbf{E}[{{\lambda}_{j}^{\prime}(x,t) N_{j}(x,t)}] \left( \mathbf{E}[{N_{j}(t)}]\right)^{-1} \\ &=& \mathbf{E}\left[{{\lambda}_{j}^{\prime}(x,t)}\right] + \text{\textbf{cov}}\left( {\lambda}_{j}^{\prime}(y,t), \frac{N_{j}(y,t)}{\mathbf{E}[{N_{j}(t)}]} \right). \end{array} $$

(27)

The first term represents the average fitness that an individual would have if it were placed in random location, and the second term describes the impact of spatial structure (Chesson 2000). For example, if species j is rare but highly aggregated, then $\mathbf {E}[{{\lambda }_{j}^{\prime }(x,t)}]$ will be positive (i.e., the species is rare, so a random adult is unlikely to be in a high-mortality area), but $\text {\textbf {cov}}\left ({\lambda }_{j}^{\prime }(y,t), N_{j}(y,t)/\mathbf {E}[{N_{j}(t)}] \right )$ will be negative (i.e., most of the individuals are aggregated in high-mortality areas). We then substitute Eq. 26 into this. The term $\ln \left \{ Y_{j} \right \}$ and the spatial average terms are already constants, and they will not contribute to the covariance. Thus,

$$ \begin{array}{@{}rcl@{}} &&{\tilde{\lambda}}_{j}^{\prime}(t) \approx \ln \left\{ Y_{j} \right\} + \mathbf{E}[{\ln \left\{ S_{j}(y,t) \right\}}]-\mathbf{E}[{\ln \left\{ C(y,t) \right\}]}\\ &+&\! L^{2} \mathbf{E}\left[{\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right)}\right]\\ &-&\! L^{2}\mathbf{E}\left[\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right)\right]\\ &+& \! L^{2} \text{\textbf{cov}}\left( \text{\textbf{cov}}\left( {d_{j}(|y-x|)}, {\ln \left\{ S_{j}(y,t) \right\}} \right), \frac{N_{j}(x,t)}{\mathbf{E}[{N_{j}(t)}]} \right)\\ &-&\! L^{2} \text{\textbf{cov}}\left( \text{\textbf{cov}}\left( {d_{j}(|y-x|)}, {\ln \left\{ C(y,t) \right\}} \right), \frac{N_{j}(x,t)}{\mathbf{E}[{N_{j}(t)}]} \right) \end{array} $$

(28)

(where the outer covariances are over all sites x). This can be simplified further by approximating survival, S_j(x, t). By Eq. 5,

$$ \begin{array}{@{}rcl@{}} &&\ln \left\{ S_{j}(x,t) \right\} = \ln \left\{ I_{j} \right\} +\\ &&\ln \left\{ 1- \sum\limits_{y=1}^{L^{2}}p_{j}d_{\text{{p}}-j}\left( |y-x|\right) N_{j}(y,t) \right\} . \end{array} $$

(29)

p_j is O(𝜖); therefore,

$$ \ln \left\{ S_{j}(x,t) \right\} \approx \ln \left\{ I_{j} \right\} - \sum\limits_{y=1}^{L^{2}}p_{j}d_{\text{{p}}-j}\left( |y-x|\right) N_{j}(y,t). $$

(30)

The above equation is the reason that our results will hold if S_j(x, t) is an approximately linear function of conspecific density (as stated in the main text); this is why we felt justified modeling S_j(x, t) using a log-likelihood function in our simulations (described below in Appendix A3, Eq. 65). Thus,

$$ \begin{array}{@{}rcl@{}} &&{\tilde{\lambda}}_{j}^{\prime}(t) \approx \ln \left\{ Y_{j} \right\} -\mathbf{E}\left[{\ln \left\{ C(y,t) \right\}}\right]\\ &+& \! \ln \left\{ I_{j} \right\} - \mathbf{E}\left[{\sum\limits_{y=1}^{L^{2}}p_{j}d_{\text{{p}}- j}\left( |y-x|\right) N_{j}(y,t)} \right]\\ &+& \! L^{2} \mathbf{E}\left[\text{\textbf{cov}}\left( d_{j}(|y- x|), \ln \left\{ S_{j}(y,t ) \right\} \right)\right]\\ &-&\! L^{2}\mathbf{E}\left[{\text{\textbf{cov}}\left( d_{j}(|y- x), \ln \left\{ C(y,t) \right\} \right)}\right]\\ &+& \! L^{2} \text{\textbf{cov}}\left( \text{\textbf{cov}}\left( {d_{j}(|y- x|)}, {\ln \left\{ S_{j}(y,t) \right\}} \right), \frac{N_{j}(x,t)}{\mathbf{E}[{N_{j}(t)}]} \right)\\ &-& \! L^{2} \text{\textbf{cov}}\left( \text{\textbf{cov}}\left( {d_{j}(|y- x|)}, {\ln \left\{ C(y,t) \right\}} \right), \frac{N_{j}(x,t)}{\mathbf{E}[{N_{j}(t)}]} \right) \end{array} $$

(31)

A1.3 Invasion analysis and partitioning terms

Here, we perform an invasion analysis. We assume that species i is rare, and the other R − 1 resident species are at equilibrium in the absence of i. We assume that all species are at a long-term steady state and stable spatial distribution. Because the residents are at equilibrium, $\tilde {\lambda }_{r}^{\prime }(t)=0$ for all r≠i. Thus,

$$ {\tilde{\lambda}}_{i}^{\prime}(t) = {\tilde{\lambda}}_{i}^{\prime}(t) - \frac{1}{R-1} \underset{r\neq i}{\sum} \tilde{\lambda}_{r}^{\prime}(t) $$

(32)

(where the summation is over all resident species). Substituting approximations Eqs. 30 and 31 into the above equation, we can partition the invader growth rate into a series of additive terms.

We combine the terms $\ln \left \{ Y_{j} \right \}$ and $\ln \left \{ I_{j} \right \}$ (giving us $\ln \left \{ Y_{j} I_{j} \right \}$). Taking the invader-resident difference in this, we are left with ΔY_i in Eq. 7. The invader-resident difference in the fourth term gives us ΔP_i, Eq. 8. The invader-resident difference in the fifth term gives us Δκ_Pi, Eq. 10. The invader-resident difference in the sixth term gives Δκ_Ci, Eq. 13. The invader-resident difference in the seventh term gives ${\varDelta } {\kappa }_{\text {{P}}i}^{\prime }$, Eq. 13. Finally, the invader-resident difference in the last term gives ${\varDelta } {\kappa }_{\text {{C}}i}^{\prime }$, Eq. 14. Together, this gives us Eq. 6.

Note that Eq. 6 does not have a term comparing the $\mathbf {E}[{\ln \left \{ C(y,t) \right \}}]$ values, as these are the same for all species, and therefore $\mathbf {E}[{\ln \left \{ C(y,t) \right \}}] - \frac {1}{R-1}\underset {r\neq i}{\sum } \mathbf {E}[{\ln \left \{ C(y,t) \right \}}] = 0$.

Appendix A2: Simplifying the Δ terms

In this section, we simplify the Δ terms, to give the approximations in Table 2. Additionally, we break them into their stabilizing components (their average over all species), and their mean fitness-difference components (the value for a given species minus the average).

To simplify notation, we define the following:

$$ \begin{array}{@{}rcl@{}} &&\widetilde{{d}_{j}^2 } = \mathbf{E}\left[{\sum\limits_{y=1}^{L^{2}} \left( d_{j}(|x-y|)\right)^{2} }\right]\\ &&\widetilde{d_{j}d_{\text{{P}},j}} = \mathbf{E}\left[{\sum\limits_{y=1}^{L^{2}} d_{j}(|x-y|)d_{\text{{p}}-j}\left( |x-y|\right) }\right]\\ &&\widetilde{d_{j} {d}_{\text{{P}},j}^{2}} = \mathbf{E}\left[{\sum\limits_{y=1}^{L^{2}} d_{j}(|x-y|)\left( d_{\text{{p}}-j}\left( |x-y|\right)\right)^{2} }\right] . \end{array} $$

(33)

Thus, $\widetilde {{d}_{j}^2 }$ is the sum of d_j(|x − y|)² over all sites (i.e., the dot product of the dispersal kernel with itself), $\widetilde {d_{j}d_{\text {{P}},j}}$ is the sum of the product of d_j(|x − y|) and $d_{\text {{p}}-j}\left (|x-y|\right )$ over all sites (i.e., the dot product of the dispersal kernel and the enemy distribution kernel), and $\widetilde {d_{j} {d}_{\text {{P}},j}^{2}}$ is the sum of the product of d_j(|x − y|) and $(d_{\text {{p}}-j}\left (|x-y|\right ))^{2}$ (i.e., the dot product of the dispersal kernel with the square of the enemy distribution kernel). These terms simplify notation below. For example, we show below that the probability that a seedling is killed by its parent’s natural enemies is related to $\widetilde {d_{j}d_{\text {{P}},j}}$.

A2.1 Appendix: Simplifying ΔY_j

We first simplify ΔY_i, Eq. 7, which is the invader-resident difference in the term $\ln \left \{ Y_{j}I_{j} \right \}$. To do this, we rewrite it such that the summation is over all species,

$$ \begin{array}{@{}rcl@{}} {\varDelta} Y_{i} &=& \ln \left\{ Y_{i}I_{i} \right\} - \frac{1}{R-1} \underset{r\neq i}{\sum} \ln \left\{ Y_{r} I_{r} \right\}\\ &=&\ln \left\{ Y_{i}I_{i} \right\} + \frac{\ln \left\{ Y_{i}I_{i} \right\}}{R-1} \\ &-& \frac{\ln \left\{ Y_{i}I_{i} \right\}}{R-1} - \frac{1}{R-1} \underset{r\neq i}{\sum} \left( \ln \left\{ Y_{r}I_{r} \right\} \right)\\ &=& \frac{(R-1)\left( \ln \left\{ Y_{i}I_{i} \right\}\right) + \left( \ln \left\{ Y_{i}I_{i} \right\}\right)}{R-1}\\ &-& \frac{1}{R-1} \sum\limits_{j=1}^{R} \left( \ln \left\{ Y_{j}I_{j} \right\} \right) \\ &=& \frac{R\left( \ln \left\{ Y_{i}I_{i} \right\}\right)}{R-1} - \frac{R}{R-1} \frac{1}{R}\sum\limits_{j=1}^{R} \left( \ln \left\{ Y_{j}I_{j} \right\} \right) \\ &=& \frac{R}{R-1} \left( \ln \left\{ Y_{i}I_{i} \right\} - \overline{ \ln \left\{ Y_{j}I_{j} \right\} } \right). \end{array} $$

(34)

When we average ΔY_i over all species, the $\ln \left \{ Y_{j}I_{j} \right \}$ terms will cancel with $\overline {\ln \left \{ Y_{j}I_{j} \right \} }$, leaving

$$ \overline{{\varDelta} Y} =0. $$

(35)

The fitness effect for species j, ${\varDelta } Y_{j} - \overline {{\varDelta } Y}$, can thus be written with either Eq. 7 or 34.

A2.2 Simplifying ΔP_j

We next simplify ΔP_i, Eq. 8, which comes from differences in $\mathbf {E}\left [{\sum \limits _{y=1}^{L^{2}}p_{j}d_{\text {{p}}-j}\left (|y-x|\right ) N_{j}(y,t) }\right ]$. The only terms that vary spatially are the N_j(y, t) terms. Thus,

$$ \begin{array}{@{}rcl@{}} {\varDelta} P_{i} &=& - \sum\limits_{y=1}^{L^{2}}p_{i} d_{\text{{p}}-i}\left( |y-x|\right) \mathbf{E}\left[{N_{i}(x,t)}\right]\\ &+& \frac{1}{R-1} \underset{r\neq i}{\sum} \sum\limits_{y=1}^{L^{2}}p_{r} d_{\text{{p}}-r}\left( |y-x|\right) \mathbf{E}\left[{N_{r}(x,t)}\right]\\ &=& - p_{i} \mathbf{E}\left[{N_{i}(x,t)}\right] \sum\limits_{y=1}^{L^{2}} d_{\text{{p}}-i}\left( |y-x|\right) \\ &+& \frac{1}{R-1} \underset{r\neq i}{\sum}p_{r} \mathbf{E}\left[{N_{r}(x,t)}\right] \sum\limits_{y=1}^{L^{2}} d_{\text{{p}}-r}(|y-x|). \end{array} $$

(36)

$\mathbf {E}\left [{N_{i}(t)}\right ]=0$, and d_j(|y − x|) sums to 1, therefore

$$ {\varDelta} P_{i} = \frac{1}{R-1} \underset{r\neq i}{\sum} p_{r} \mathbf{E}\left[{N_{r}(t)}\right]. $$

(37)

Because $\underset {r\neq i}{\sum } \mathbf {E}\left [{N_{r}(t)}\right ]=1$, this term is simply $\frac {1}{R-1}$ times the frequency-weighted mean of p_r. This can rewritten as

$$ {\varDelta} P_{i} \approx \frac{1}{(R-1)^{2}} \underset{r\neq i}{\sum} p_{r} + {\textbf{cov}}^{\text{{S}}}\left( \mathbf{E}\left[{N_{r}(t)}\right], p_{r} \right)^{\{-i\}}, $$

(38)

where the covariance is over all species except species i. Thus, it tends to be smaller if the most abundant species are the least affected by their natural enemies.

Next, we calculate the stabilizing mechanism. The average of $\frac {1}{R-1} \underset {r\neq i}{\sum } p_{r}$, when averaged over all species, is simply the mean of p_r across species, $\overline {p_{j}}$. Previous work has suggested that the average of $\text {\textbf {cov}}^{\text {{S}}}\left (\mathbf {E}[{N_{r}(t)}], p_{r} \right )^{\{-i\}}$ can be approximated as $\text {\textbf {cov}}^{\text {{S}}}\left (\mathbf {E}[{N_{j}(t)}], p_{j} \right )$, calculated over all species when they are at equilibrium (Stump and Comita 2018). Thus,

$$ \overline{{\varDelta} P} \approx \frac{\overline{p_{j}}}{R-1} + \text{\textbf{cov}}^{\text{{S}}}\left( \mathbf{E}[{N_{j}(t)}], p_{j} \right). $$

(39)

Additionally, previous work has suggested that the difference in the $\text {\textbf {cov}}^{\text {{S}}}\left (\mathbf {E}[{N_{j}(t)}], p_{j} \right )$ terms will have a negligible effect on fitness-differences (Stump and Comita 2018). Thus, the fitness effect of ΔP_j is

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} P_{j} - \overline{{\varDelta} P} \approx \frac{1}{(R-1)^{2}}\underset{r\neq i}{\sum} p_{r} - \frac{\overline{p_{j}}}{R-1}\\ &=&- \frac{p_{i}}{(R-1)^{2}} + \frac{p_{i}}{(R-1)^{2}}\\ &+& \frac{1}{(R-1)^{2}}\underset{r\neq i}{\sum} p_{r} - \frac{\overline{p_{j}}}{R-1} \\ &=&- \frac{p_{i}}{(R-1)^{2}} + \frac{1}{(R-1)^{2}}\sum\limits_{j=1}^{R} p_{j} - \frac{\overline{p_{j}}}{R-1} \\ &=&- \frac{p_{i}}{(R-1)^{2}} + \frac{\overline{p_{j}} R }{(R-1)^{2}} - \frac{\overline{p_{j}}}{R-1} \\ &=&- \frac{p_{i}}{(R-1)^{2}} + \frac{\overline{p_{j}} R - \overline{p_{j}} (R-1)}{(R-1)^{2}}\\ &=& -\frac{p_{i} - \overline{p_{j}}}{(R-1)^{2}}. \end{array} $$

(40)

As in previous work (Stump and Comita 2018), this effect is generally minor.

A2.3 Simplifying Δκ_Pj

We next simplify Δκ_Pi, Eq. 10, which begins as

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} \kappa_{\text{P}i} = L^{2} \mathbf{E}\left[\text{\textbf{cov}}\left( d_{i}(|y-x|), \ln \left\{ S_{i}(y,t) \right\} \right)\right]\\ &-& \frac{1}{R-1} \underset{r\neq i}{\sum} L^{2} \mathbf{E}\left[\text{\textbf{cov}}\left( d_{r}(|y-x|), \ln \left\{ S_{r}(y,t) \right\} \right)\right] \end{array} $$

(41)

We first calculate the covariance $\text {\textbf {cov}}\left (d_{j}(|y-x|), \ln \left \{ S_{j}(y,t) \right \} \right )$. By the law of total covariance,

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right) = \\&& \mathbf{E}[{\text{\textbf{cov}}\left( \left. d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right\vert Z\right)}]\\ &+& \text{\textbf{cov}}\left( \mathbf{E}\left[ \left.d_{j}(|y-x|)\right| Z \right], \mathbf{E}\left[ \left.\ln \left\{ S_{j}(y,t) \right\}\right| Z \right] \right), \end{array} $$

(42)

where $\text {\textbf {cov}}\left (\left . X, Y \right \vert Z\right )$ is the covariance of X and Y conditioned on some condition Z (the expectation is then over all possible Z’s). We condition our covariances on |y − x| = z (i.e., on the assumption that a seed dispersed a distance z); for notational simplicity, we write this as conditioning on z. The first term is 0, because d_j(|y − x|) will be constant when z is constant. Similarly, $\mathbf {E}\left [ \left .d_{j}(|x-y|)\right | z \right ] = d_{j}(z)$. Finally, we use our approximation in Eq. 30 to write $\ln \left \{ S_{j}(y,t) \right \}$ as $\ln \left \{ Y_{j} \right \}-{\sum }_{y=0}^{L^{2}}p_{j}d_{\text {{p}}-j}\left (|y-x|\right ) N_{j}(|y-x|,t)$. The term $\ln \left \{ Y_{j} \right \}$ is constant and thus will not affect the covariance. Thus,

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right) =\\ &-&\text{\textbf{cov}}\left( d_{j}(z), \mathbf{E}\left[ \left.\sum\limits_{y=0}^{L^{2}}p_{j}d_{\text{{p}}-j}\left( |y-x|\right) N_{j}(|y-x|,t)\right| z \right] \right). \end{array} $$

(43)

We next find a formula for the mean impact of natural enemies averaged over space, $\mathbf {E}\left [ \left .{\sum }_{y=1}^{L^{2}}p_{j}d_{\text {{p}}-j}\left (|y-x|\right ) N_{j}(y,t)\right | z \right ]$. We are averaging over all possible sites x, so spatial structure does not matter (this is taken into account with ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$). Let us define ψ_z as the number of sites a distance z from x. Let us first begin with the z that we are conditioning on, z = |x − y|. We know that at least one of those sites must contain a conspecific adult (i.e., the parent of the seed at x). The remaining ψ_z − 1 sites will contain a conspecific adult with probability E[N_j(t)]. Thus, there will on average be ψ_zE[N_j(t)] + 1 −E[N_j(t)] adults a distance z away (and they will reduce survival by $p_{j}d_{\text {{p}}-j}\left (z\right )$). At any other distance w≠|x − y|, there will be ψ_w sites, each will contain a conspecific adult with probability E[N_j(t)], and they will each reduce survival by $p_{j}d_{\text {{p}}-j}\left (w\right )$. Thus,

$$ \begin{array}{@{}rcl@{}} &&\mathbf{E}\left[ \left.\sum\limits_{y=1}^{L^{2}}p_{j}d_{\text{{p}}-j}\left( |y-x|\right) N_{j}(y,t)\right| z \right]\\ &=& p_{j} d_{\text{{p}}-j}\left( z\right) (1-\mathbf{E}[{N_{j}(t)}]) + \sum\limits_{y=1}^{L^{2}} p_{j} d_{\text{{p}}-j}\left( y-x\right) \mathbf{E}[{N_{j}(t)}]\\ &=& p_{j} d_{\text{{p}}-j}\left( z\right) (1-\mathbf{E}[{N_{j}(t)}]) + p_{j} \mathbf{E}[{N_{j}(t)}] \end{array} $$

(44)

Only the first term varies with z. Therefore,

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( \mathbf{E}\left[ \left.d_{j}(|y-x|)\right| z \right], \mathbf{E}\left[ \left.\ln \left\{ S_{j}(y,t) \right\}\right| z \right] \right) \\ &=& \text{\textbf{cov}}\left( d_{j}(z), p_{j} d_{\text{{p}}-j}\left( z\right) (1-\mathbf{E}[{N_{j}(t)}]) + p_{j} \mathbf{E}[{N_{j}(t)}] \right)\\ &=& \left( 1-\mathbf{E}[{N_{j}(t)}]\right)p_{j} \text{\textbf{cov}}\left( d_{j}(z), d_{\text{{p}}-j}\left( z\right) \right) \end{array} $$

(45)

The covariance term $\text {\textbf {cov}}\left (d_{j}(z), d_{\text {{p}}-j}\left (z\right ) \right )$ is $p_{j} \mathbf {E}[{d_{j}(z)d_{\text {{p}}-j}\left (z\right )}]$ − $p_{j} \mathbf {E}[{d_{j}{z}}]\mathbf {E}[{d_{\text {{p}}-j}\left (z\right )}]$. The mean of the product is $\frac {1}{L^{2}}\widetilde {d_{j}d_{\text {{P}},j}}$ (by our definition in Eq. 33). The mean of d_j(z) and $d_{\text {{p}}-j}\left (z\right )$ are each L², as they are dispersal kernels and therefore must sum to 1 over all L² sites. Thus,

$$ \begin{array}{@{}rcl@{}} &&\mathbf{E}[{\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right)}]\\ &=& \frac{1}{L^{2}} \left( 1-\mathbf{E}[{N_{j}(t)}]\right) \left( p_{j}\widetilde{d_{j}d_{\text{{P}},j}} - \frac{1}{L^{2}} p_{j} \right). \end{array} $$

(46)

If L is large, we can ignore the $-\frac {1}{L^{2}}p_{j}$ term.

Substituting this into Eq. 10

$$ \begin{array}{@{}rcl@{}} {\varDelta} \kappa_{\text{{P}}i} &\approx& -p_{i} \left( 1-\mathbf{E}\left[{N_{i}(t)}\right]\right) \widetilde{d_{i}d_{\text{{P}},i}}\\ &&+ \frac{1}{R-1}\underset{r\neq i}{\sum} p_{r}\left( 1-\mathbf{E}\left[{N_{r}(t)}\right]\right) \widetilde{d_{r}d_{\text{{P}},r}}\\ &\approx& \frac{-1}{R-1}\underset{r\neq i}{\sum}p_{r} \mathbf{E}\left[{N_{r}(t)}\right] \widetilde{d_{r}d_{\text{{P}},r}} \\ &&-\left( p_{i} \widetilde{d_{i}d_{\text{{P}},i}} - \frac{1}{R-1}\underset{r\neq i}{\sum} p_{r}\widetilde{d_{r}d_{\text{{P}},r}} \right). \end{array} $$

(47)

Using a procedure similar to Eq. 34, this becomes

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} \kappa_{\text{{P}}i} \approx \frac{-1}{R-1}\underset{r\neq i}{\sum}p_{r} \mathbf{E}\left[{N_{r}(t)}\right] \widetilde{d_{r}d_{\text{{P}},r}} \\ &-&\frac{R}{R-1}\left( p_{i} \widetilde{d_{i}d_{\text{{P}},i}} - \overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}} \right). \end{array} $$

(48)

Thus, averaging across species, this becomes a frequency-weighted mean of $\overline {p_{j}\widetilde {d_{j}d_{\text {{P}},j}}}$. We believe that the approximation for $\overline {{\varDelta } P}$ will hold (i.e., that we can approximate the mean of the covariance over invaders as the covariance when all species are at equilibrium); in this case,

$$ \overline{{\varDelta} \kappa_{\text{{P}} }} \approx \frac{-\overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}}}{R-1} - \text{\textbf{cov}}^{\text{{S}}}\left( \mathbf{E}[{N_{j}(t)}], p_{j} \widetilde{d_{j}d_{\text{{P}},j}} \right). $$

(49)

Additionally, if the covariances cancel out of the mean fitness differences, then

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} \kappa_{\text{{P}}j} - \overline{{\varDelta} \kappa_{\text{{P}} }} \approx \left( \frac{-1}{R-1}\underset{r\neq i}{\sum}p_{r} \widetilde{d_{r}d_{\text{{P}},r}}\mathbf{E}\left[{N_{r}(t)}\right] \right.\\ &&\left. -\frac{R}{R-1}\left( p_{j} \widetilde{d_{j}d_{\text{{P}},j}} - \overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}} \right) \right) + \frac{\overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}}}{R-1} . \end{array} $$

(50)

E[N_i(t)] = 0, and therefore the summation $\underset {r\neq i}{\sum }p_{r} \widetilde {d_{r}d_{\text {{P}},r}}\mathbf {E}[{N_{r}(t)}]$ over all residents is equal to the summation $\sum \limits _{j=1}^{R}p_{j} \widetilde {d_{j}d_{\text {{P}},j}}\mathbf {E}[{N_{j}(t)}]$ over all species. This term is $\overline {p_{j} \widetilde {d_{j}d_{\text {{P}},j}}}$ $+R\text {\textbf {cov}}^{\text {{S}}}\left (p_{j} \widetilde {d_{j}d_{\text {{P}},j}}, \mathbf {E}[{N_{j}(t)}] \right )$, although if we can cancel out the covariance term, then

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} \kappa_{\text{{P}}j} - \overline{{\varDelta} \kappa_{\text{{P}} }} \approx \left( -\frac{\overline{p_{j} \widetilde{d_{j}d_{\text{{P}},j}}}}{R-1} -\frac{R}{R-1}\left( p_{j} \widetilde{d_{j}d_{\text{{P}},j}} - \overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}} \right) \right) \\ &+& \frac{\overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}}}{R-1} \\ &=&-\frac{R}{R-1}\left( p_{j} \widetilde{d_{j}d_{\text{{P}},j}} - \overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}}\right). \end{array} $$

(51)

A2.4 Simplifying Δκ_Cj

We next simplify Δκ_Ci, Eq. 13. Using the law of total covariance and conditioning on |y − x| = z (as before),

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right) = \\ &&\text{\textbf{cov}}\left( d_{j}(z), \mathbf{E}\left[ \left.\ln \left\{ C(y,t) \right\}\right| z \right] \right). \end{array} $$

(52)

We next find a formula for $\mathbf {E}\left [ \left .\ln \left \{ C(y,t) \right \}\right | z \right ]$. We define C(j, y, t) as the number of seedlings at y of species j; it is

$$ C(j,y,t) = \sum\limits_{x=1}^{L^{2}} d_{j}(|x-y|) Y_{j} S_{j}(y,t) N_{j}(x,t). $$

(53)

Thus, $C(y,t)=\sum \limits _{j=1}^{R} C(j,y,t)$. C(j, y, t) can be rewritten as

$$ \begin{array}{@{}rcl@{}} &&C(j,y,t) = Y_{j} \left( I_{j} \left( 1 - p_{j}\sum\limits_{x=1}^{L^{2}} d_{\text{{p}}-j}\left( |x-y|\right) N_{j}(x,t) \right) \right) \times\\ &&\left( \sum\limits_{x=1}^{L^{2}} d_{j}(|x-y|) N_{j}(x,t) \right). \end{array} $$

(54)

This will be easier to work with if its value is close to E[N_j(t)]. Thus, we next divide this divide both sides by the across-species mean of Y_jI_j, which we call (Y I)^∗. Then, we define $\widehat {YI}_{j} = Y_{j}I_{j}/(YI)^{*} - 1$, and

$$ \begin{array}{@{}rcl@{}} &&\frac{C(j,y,t)}{(YI)^{*}} = \left( 1+\widehat{YI}_{j} \right)\left( 1 - p_{j}\sum\limits_{x=1}^{L^{2}} d_{\text{{p}}-j}\left( |x-y|\right) N_{j}(x,t) \right) \times\\ &&\left( \sum\limits_{x=1}^{L^{2}} d_{j}(|x-y|) N_{j}(x,t) \right). \end{array} $$

(55)

We next approximate the conditional average of C(j, y, t) given an N_j(x, t) for some location x. N_j(x, t) it could be 0 or 1, we will define this later. We approximate the mean of C(j, y, t) as the value of C(j, y, t) when N_j(x, t) = E[N_j(t)]; the true mean would include an effect of the variance of N_j(x, t); however, we find that the approximations that this simplification produces fit well (see Appendix A4). Thus,

$$ \begin{array}{@{}rcl@{}} &&\frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t) \right]}{(YI)^{*}} \approx \left( 1+\widehat{YI}_{j} \right) \times\\ &&\left( 1 - p_{j} \mathbf{E}[{N_{j}(t)}] - \left( N_{j}(x,t)-\mathbf{E}[{N_{j}(t)}] \right)p_{j} d_{\text{{p}}-j}\left( |x-y|\right) \right) \times \\ &&\left( \mathbf{E}[{N_{j}(t)}] + \left( N_{j}(x,t)-\mathbf{E}[{N_{j}(t)}] \right) d_{j}(|x-y|) \right). \end{array} $$

(56)

Next, we multiply these terms. The term $\widehat {YI}_{j}$ is O(𝜖); thus, $\widehat {YI}_{j} p_{j}$ will be small, and we will ignore it. If N_j(x, t) = 0, this becomes

$$ \begin{array}{@{}rcl@{}} &&\frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t)=0 \right]}{(YI)^{*}} \approx \\ &&\mathbf{E}[{N_{j}(t)}] \left[1 - d_{j}(|x-y|) + \widehat{YI}_{j} \left( 1 - d_{j}(|x-y|) \right) \right.\\ &&\left. -p_{j} \mathbf{E}[{N_{j}(t)}] \left( 1 + d_{\text{{p}}-j}\left( |x-y|\right) \right) \left( 1 - d_{j}(|x-y|)\right) \right]. \end{array} $$

(57)

If N_j(x, t) = 1,

$$ \begin{array}{@{}rcl@{}} &&\frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t)=1 \right]}{(YI)^{*}} \approx \\ &&\frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t)=0 \right]}{(YI)^{*}} + d_{j}(|x-y|) \\ &+& \widehat{YI}_{j} d_{j}(|x-y|) - p_{j} d_{j}(|x-y|) d_{\text{{p}}-j}\left( |x-y|\right) \\ &-&p_{j}d_{j}(|x-y|) \mathbf{E}[{N_{j}(t)}] \left( 1 - d_{\text{{p}}-j}\left( |x-y|\right) \right) \\ &-&p_{j} d_{\text{{p}}-j}\left( |x-y|\right) \mathbf{E}[{N_{j}(t)}] \left( 1 - d_{j}(|x-y|) \right). \end{array} $$

(58)

For notational simplicity, we will define ξ_j(|x − y|, t) as the amount that a seed of species j would change seedling competition at y by being at x,

$$ \begin{array}{@{}rcl@{}} &&\xi_{j}(|x-y|,t) = \frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t)=1 \right]}{(YI)^{*}}\\ &-& \frac{\mathbf{E}\left[ \left.C(j,y,t)\right| N_{j}(x,t)=0 \right]}{(YI)^{*}}\\ &\approx& - p_{j} \left( d_{j}(|x-y|) d_{\text{{p}}-j}\left( |x-y|\right) \right.\\ &-& 2 \mathbf{E}[{N_{j}(t)}] d_{j}(|x-y|) d_{\text{{p}}-j}\left( |x-y|\right) \\ &&\left. + \mathbf{E}[{N_{j}(t)}] d_{j}(|x-y|) + \mathbf{E}[{N_{j}(t)}] d_{\text{{p}}-j}\left( |x-y|\right)\right) \\ &+& d_{j}(|x-y|) + \widehat{YI}_{j} d_{j}(|x-y|). \end{array} $$

(59)

We can write the mean of C(y, t) using the C(j, y, t) terms

$$ \frac{\mathbf{E}\left[{C(y,t)}\right]}{(YI)^{*}} = \sum\limits_{j=1}^{R} \frac{\mathbf{E}\left[{C(k,y,t)}\right]}{(YI)^{*}} $$

(60)

(summed over all species k). If we focus on a particular site x, then this could also be written as

$$ \begin{aligned} \frac{\mathbf{E}{C(y,t)}}{(YI)^{*}} =& \sum\limits_{k=1}^{R} \frac{\mathbf{E}\left[ \left.C(k,y,t)\right| N_{k}(x,t)=0 \right]}{(YI)^{*}} \\ & + \sum\limits_{k=1}^{R} \mathbf{E}\left[{N_{k}(t)}\right]\xi_{j}(|x-y|,t), \end{aligned} $$

(61)

i.e., the first term is the impact of all sites other than x, and the second term is the specific contribution that site x has on seedling competition at y. If we knew that species j occupied x, then

$$ \begin{array}{@{}rcl@{}} \frac{\mathbf{E}\left[ \left.C(y,t)\right| N_{j}(x,t)=1 \right]}{(YI)^{*}} &=& \sum\limits_{k=1}^{R} \frac{\mathbf{E}\left[ \left.C(k,y,t)\right| N_{j}(x,t)=0 \right]}{(YI)^{*}}\\ &\quad+& \xi_{j}(|x-y|,t). \end{array} $$

(62)

Thus, combining the above two equations,

$$ \begin{aligned} \frac{\mathbf{E}\left[ \left.C(y,t)\right| N_{j}(x,t)=1 \right]}{(YI)^{*}} =&\frac{\mathbf{E}[{C(y,t)}]}{(YI)^{*}} + \xi_{j}(|x-y|,t) \\ &- \sum\limits_{k=1}^{R} \mathbf{E}[{N_{k}(t)}]\xi_{j}(|x-y|,t). \end{aligned} $$

(63)

The terms E[C(y, t)] and (Y I)^∗ are within O(𝜖), and thus their ratio will be close to 1. Additionally, the ξ_j(|x − y|, t) terms be within O(𝜖). Thus, if we take the natural log of both sides and cancel the $\ln \left \{ (YI)^{*} \right \}$ terms, then we are left with

$$ \begin{array}{@{}rcl@{}} &&\ln \left\{ \mathbf{E}\left[ \left.C(y,t)\right| N_{j}(x,t)=1 \right] \right\} \approx \ln \left\{ \mathbf{E}[{C(y,t)}] \right\}\\ &&\qquad+ \xi_{j}(|x-y|,t) - \sum\limits_{j=1}^{R} \mathbf{E}[{N_{j}(t)}]\xi_{j}(|x-y|,t). \end{array} $$

(64)

Plugging (64) into our formula for the covariance, Eq. 52

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right) \approx \\ &&\text{\textbf{cov}}\left( d_{j}(z), \xi_{j}(d,t) - \sum\limits_{k=1}^{R} \mathbf{E}[{N_{j}(t)}]\xi_{k}(d,t) \right). \end{array} $$

(65)

(Note that the $\ln \left \{ \mathbf {E}[{C(y,t)} \right \}]$ terms cancel, as they are constant). Therefore,

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} \kappa_{\text{{C}}i} \approx - \text{\textbf{cov}}\left( d_{i}(z), \xi_{i}(z,t) - \sum\limits_{j=1}^{R} \mathbf{E}[{N_{j}(t)}]\xi_{j}(z,t) \right)\\ &+& \frac{1}{R-1} \underset{r\neq i}{\sum} \text{\textbf{cov}}\left( d_{r}(z), \xi_{r}(z,t) - \sum\limits_{j=1}^{R} \mathbf{E}[{N_{j}(t)}]\xi_{j}(z,t) \right). \end{array} $$

(66)

Our attempts at simplifying this equation further resulted in a poor approximation. Thus, we leave it as this. The stabilizing effect is thus the above equation averaged over species, and the fitness effect is the above equation minus the across-species average.

A2.5 Simplifying the ${\varDelta } \kappa ^{\prime }$ terms

The terms ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ account for the effects of spatial aggregation. In order to find a simple approximation, we would need an approximation for the spatial structure of each species. We were unable to derive this from first principals. Thus, we need to study its behavior using computer simulations. Here, we outline how these terms can be calculated from simulations.

We begin with ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$. The average probability that a seed would survive can be determined by calculating the mean survival of seeds produced by each adult, and is

$$ \frac{1}{Y\mathbf{E}[{N_{j}(t)}]L^{2}} \sum\limits_{x=1}^{L^{2}} Y_{j}N_{j}(x,t) \sum\limits_{y=1}^{L^{2}}d_{j}(|y-x|)S_{j}(y,t). $$

(67)

Thus, the per-seed average of $\ln \left \{ S_{j}(y,t) \right \}$ can be simplified to

$$ \begin{array}{@{}rcl@{}} &&\left( \text{per-seed mean} \ln \left\{ S_{j}(y,t) \right\} \right) \\ &=& \frac{1}{Y\mathbf{E}[{N_{j}(t)}]L^{2}} \sum\limits_{x=1}^{L^{2}} Y_{j}N_{j}(x,t) \sum\limits_{y=1}^{L^{2}}d_{j}(|y-x|)\ln \left\{ S_{j}(y,t) \right\}\\ &=& \mathbf{E}[{\ln \left\{ S_{j}(y,t) \right\}}] + \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right)\\ &+& \text{\textbf{cov}}\left( d_{j}(|y-x|)\ln \left\{ S_{j}(y,t) \right\}, N_{j}(x,t) \right). \end{array} $$

(68)

Thus,

$$ \begin{array}{@{}rcl@{}} &&\text{\textbf{cov}}\left( d_{j}(|y-x|)\ln \left\{ S_{j}(y,t) \right\}, N_{j}(x,t) \right)\\ &=& \left( \text{per-seed mean} \ln \left\{ S_{j}(y,t) \right\} \right) \\ &-& \mathbf{E}[{\ln \left\{ S_{j}(y,t) \right\}}] - \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(y,t) \right\} \right) \end{array} $$

(69)

The effect of I_j will cancel out of the mean log seed survival and $\mathbf {E}[{\ln \left \{ S_{j}(y,t) \right \}}]$. Putting this into our resident-invader comparison, we find that

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} {\kappa}_{\text{{P}}i}^{\prime} \approx \left[ \left( \text{per-seed mean}~ \ln \left\{ S_{i}(y,t) \right\} \right) \right.\\ &-&\left. \frac{1}{R-1}\underset{r\neq i}{\sum} \left( \text{per-seed mean}~\ln \left\{ S_{r}(x,t) \right\} \right) \right]\\ &-& {\varDelta} P_{i} - {\varDelta} \kappa_{\text{{P}}i} . \end{array} $$

(70)

This result is somewhat tautological—we wanted to quantify the species-specific differences in the per-seed average log chance of seed survival, so we defined the third term as the per-seed average log chance minus the other two terms. We did this because it was mathematically correct (by definition), and because each term represents a qualitatively different effect. And, intuitively, the covariance in ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ is any impact on the per-seed mean log survival that could not be explained by the other two effects. This equation would be useless if ΔP_i and Δκ_Pi could not be independently verified; however, as we show in the next section, we were able to run simulations where spatial aggregation was removed (i.e., so that ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ would be 0 on average). We found that, indeed, the community dynamics could be predicted fairly accurately by ΔP_j, Δκ_Pj, Δκ_Cj, and ΔY_j.

We can show similarly that the per-seed mean of $\ln \left \{ C(y,t) \right \}$ (i.e., the average amount of competition experienced by a species) is

$$ \begin{array}{@{}rcl@{}} &&\left( \text{per-seed mean}\ln \left\{ C(y,t) \right\} \right)_{j}\\ &=& \frac{1}{Y\mathbf{E}[{N_{j}(t)}]L^{2}} \sum\limits_{x=1}^{L^{2}} Y_{j}N_{j}(x,t) \sum\limits_{y=1}^{L^{2}}d_{j}(|y-x|)\ln \left\{ C(y,t) \right\}\\ &=& \mathbf{E}\left[{\ln\left\{ C(y,t) \right\}}\right] + \text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right)\\ &+& \text{\textbf{cov}}\left( d_{j}(|y-x|)\ln \left\{ C(y,t) \right\}, N_{j}(x,t) \right). \end{array} $$

(71)

Therefore,

$$ \begin{array}{@{}rcl@{}} &&{\varDelta} {\kappa}_{\text{{C}}i}^{\prime} \approx \left( \left( \text{per-seed mean}~\! \ln \left\{ C(y,t) \right\} \right)_{i} \right.\\ &&\left. - \frac{1}{R-1}\underset{r\neq i}{\sum} \left( \text{per-seed mean}~\!\ln \left\{ C(x,t) \right\} \right)_{r} \right)\! -\! {\varDelta} \kappa_{\text{{C}}i} . \end{array} $$

(72)

In Appendix A3, we run several simulations and check how ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ vary with a species mean dispersal distance and resident-invader state.

Appendix A3: Computer simulations

In this section, we describe how we implemented our computer simulations.

To speed computation, our model worked with the Manhattan distance function (i.e., the distance between (x₁, y₁) and (x₂, y₂) is |x₁ − x₂| + |y₁ − y₂|), rather than the Euclidean distance function (i.e., the distance between (x₁, y₁) and (x₂, y₂) is $\sqrt {(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}$). The area of a circle is different under each distance function, which affects how many trees are a given distance away. For example, imagine we placed point in a grid, with each point separated by 1 m; then, say we selected a focal point. Under the Manhattan distance function, there would be 20 points a distance 5 m away from the focal point; however, under the Euclidean distance function, there would be 28 points between 4.5 and 5.5 m away from the focal point. But, given that none of our analysis depends on the distance function used (and indeed, d_A,j(z) and d_j(z) are arbitrary), we do not believe that this qualitatively affected our results.

We model dispersal using a discrete approximation of a 2-dimensional t (2Dt) distribution. The probability distribution function of a 2Dt distribution with mean 0 and 2 degrees of freedom (i.e., the pdf that a seed will disperse a distance z if it is equally likely to disperse in every direction) is

$$ D(z)=\frac{2\pi z}{\pi \exp\{\alpha_{j}\} (1+z^{2} \exp\{-\alpha_{j}\})^{2}} $$

(73)

where α_j is a species-specific parameter. Note that the 2πz in the numerator is used because D(z) calculates the chance that it lands in anywhere a distance z away; this is generally left off to model the chance that it lands at a particular location a distance D(z) away. As stated above, Euclidean and Manhattan distances are different; thus, when we were generating seed dispersal kernels, we decided to preserve the probability that a seed would disperse a given distance, rather than the probability it would disperse to a given point.

In our simulations, we calculate d_j(z) as follows: A seed that does not disperse beyond l/2 (i.e., half a site) remains in its parent’s site, and a seed that disperses between l(z − 0.5) and l(z + 0.5) is dispersed to one of the 4z sites a distance z away with equal likelihood. To reduce computational load, we assume that if a seed dispersed more than 5.5l, then it was placed in a random site (this could include sites within 5 sites, though the impact was trivial if L was large). Thus,

$$ d_{j}(z) = \begin{cases} {\int}_{0}^{l/2} D(y)dy + {\int}_{5.5l}^{\infty}\frac{D(y)}{L^{2}}dy & z=0\\ \frac{1}{4z} {\int}_{l(z-0.5)}^{l(z+0.5)} D(y)dy + {\int}_{5.5l}^{\infty}\frac{D(y)}{L^{2}}dy & 1\leq z\leq 5\\ {\int}_{5.5l}^{\infty}\frac{D(y)}{L^{2}}dy & z>5 \end{cases} $$

(74)

We modeled the community as a torus (i.e., a seed that moved off one side was warped to the other side), so that there were no edge effects.

In simulation, we modeled survival such that enemies affected tho log-odds probability of survival, rather than having a linear impact on survival. We did this so that survival would never become negative, and because this connected better with empirical literature (Comita et al. 2010). However, when p_A,j and p_S,j were small, the log-odds function is close enough to linear to be approximated by Eq. 2. We defined the terms $\tilde {I}_{j}$ and $\tilde {p}_{j}$, such that survival on the log-odds scale was

$$ S_{j}(x,t)=\text{logit}\left( \tilde{I}_{j} - \sum\limits_{y=1}^{L^{2}} \tilde{p}_{j} d_{\text{{p}}-j}\left( |x-y|\right) N_{j}(y,t)\right). $$

(75)

To relate these terms to our previous definitions, $\tilde {I}_{j}$ was defined as the inverse logit of I_j (i.e., so that when there are no natural enemies, survival is I_j) and $\tilde {p}_{j}$ was defined as $p_{j}(1+\exp \{\tilde {I}_{j}\})$ (i.e., so that a small increase 𝜖 in the number of nearby conspecifics reduces survival by approximately p_j𝜖).

We modeled the impact of distance-responsive enemies as declining exponentially according an exponential distance function (Comita et al. 2010). We calculate survival at the midpoint of a site (which was 0.25 sites for the parent’s site) and ignore adults more than 5 sites away. Thus,

$$ d_{\text{{A}},j}(z)= \begin{cases} c_{1} \exp\{-0.25 c_{2}\} & z=0\\ c_{1} \exp\{-c_{2} z\} & 1\leq z\leq 5\\ 0 & z> 5 \end{cases} $$

(76)

where c₁ is a distance constant and c₂ is a constant that allows d_A,j(z) to sum to 1. In all cases below, we chose c₂ = 0.5 (to match empirical literature of 0.2 per m, Comita et al. (2010)), and c₁ = 0.752.

A3.1 Computational details for running simulations

We modeled sites as l = 10 m squares, and usually simulated a 200 × 200 site (400 ha) community. The simulation was usually initiated by filling each site with a random individual with equal likelihood. In cases where one species’ density was held constant, we randomly filled that number of sites with the focal species and filled the remainder by randomly selecting the R − 1 other species with equal probability.

Each time step, we examined each adult to determine if it died. If it survived, then there was no change; if the adult died, then we calculated the number of seeds that were expected to fall in each site. Thus, an individual of species j at site x would contribute Y_jd_j(|y − x|) seeds to site y. We then calculated the impact of distance- and density-responsive enemies. If there were s_j(y) seeds of species j at y, then every seed’s log-odds survival was reduced by $s_{j}(y)\tilde {p}_{j}^{\text {{S}}}$; similarly, an adult of species j at x reduced the log-odds survival to species j seedlings at y by ${\tilde {p}}_{j}^{\text {{A}}}{d}_{j}^{\text {{A}}}(|y-x|)$. We used these to calculate S_j(y, t) and used this to calculate the expected number of seedlings at each site (it was rarely a whole number). We then chose a species randomly in proportion to the relative fraction of seedlings at the site. This procedure was repeated each time step.

In some cases, we held the frequency of one species constant. To do this, we began each time step by recording where each individual of that species was located. We then calculated the births and deaths across the community as described above. At the end of the time step, we calculated the amount that the focal population changed in order to estimate ${\tilde {\lambda }}_{j}^{\prime }(t)$, and the values needed to estimate our Δ terms. If the population had increased by n individuals, then we randomly killed n individuals of the focal population. If the population had decreased by n, then we randomly selected n individuals that had died (i.e., using the record from the start of the time step) and restored them.

We performed invasion analyses using the following steps. First, we simulated a community for 2000 time steps, to eliminate any species that was not coexisting. We also recorded the mean frequency of each remaining species. Then, each species was selected in turn to be an invader. We removed the invader and increased every other species’ density proportionately. We then simulated the community dynamics for 100 time steps, so that it could reach equilibrium and a stable spatial structure. We then introduced a small number of invader individuals (typically 50 adults). The invaders were placed in random locations, replacing the adult that was there. We ran community dynamics for several time steps (typically 13 time steps, or 5 generations), so that the invader could build up a spatial distribution, but not so long that they became common. We then continued to run the simulation for 37 time steps, and each time step, recorded data on the invader. We calculated the expected value of ${\tilde {\lambda }}_{j}^{\prime }(t)$ by calculating the probability that species j captures each site (i.e., the number of seedlings divided by the total number of seedlings, both after accounting for natural enemies), dividing that value by E[N_j(t)], multiplying that value by δ, and adding (1 − δ). We also calculated ΔY_j, ΔP_j, Δκ_Pj, Δκ_Cj, ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$, and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$, using approximations derived in Appendix ??. We averaged our results over 37 time steps and ran 120 replicate simulations. It was rare for an invader to die off stochastically, but when it did, we only recorded growth and the Δ values when it was alive.

A3.2 Parameter sets

The six parameter sets used for figures are listed in Table 3. In all cases, we assumed that $\tilde {I}_{j}=-1$ for all species (i.e., I_j ≈ 0.27), that, δ = 0.4, and that $\tilde {p}_{\text {{A}},j}=\tilde {p}_{\text {{S}},j}$. We generated parameters for 8 species, though in some cases not all persisted. The “trade-off” parameter set assumed that there was a strict trade-off between α_j and $\tilde {p}_{j}$, such that increasing $\tilde {p}_{j}$ by 0.1 increased α_j by 0.5. The “equal sensitivity” parameter set assumed that α_j varied, but all species had the same $\tilde {p}_{j}$ value; this was used to isolate the particular effects of differences in dispersal. Similarly, the “equal dispersal” parameter set assumed that $\tilde {p}_{j}$ varied, but α_j was the same between species and was used to isolate the impact of differences in natural enemy susceptibility. In the “yield differences” parameter set, species differ in only their yield values; we also assumed that $\tilde {p}_{j}$ values were higher than in other parameter sets. We also included two “random parameters” sets, in which $\tilde {p}_{j}$, α_j, and Y_j were selected randomly (without replacement) from the same distribution; these were used to test what would occur if there was not a strict trade-off.

Table 3 Parameters in displayed figures. If only one number is given, then all species have that value for that parameter

We tested several additional parameter sets during our analysis. We included a community where species differed in $\tilde {p}_{j}$ and had near-universal dispersal, we tested two parameter sets that lacked specialist natural enemies, and we tested an additional random parameter set. We chose not to include these data, because we felt they did not add anything. We also ran similar parameter sets but changed the range over which α_j and $\tilde {p}_{j}$ could vary, and parameters where $\tilde {p}_{\text {{A}},j}\neq \tilde {p}^{\text {{S}}}_{j}$ (including $\tilde {p}^{\text {{S}}}_{j}=0$). Our results were basically similar.

Appendix A4: Checking our approximation

In this section, we describe how we used our computer simulations to check our analytical approximations. We simulated a community without spatial structure. We did this by randomly rearranging the community at the end of each time step (a method previously referred to as the “shaken” method, (Stump et al. 2018b)). This allowed us to remove the effects of spatial structure generated by clumped adult distributions (i.e., the ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ terms), without removing the impact than adult had on its offspring (i.e., the Δκ_Pj and Δκ_Cj terms).

We studied our model in two ways. First, we performed invasion analysis on our model. We calculated the ${\varDelta } {\kappa }_{\text {{P}}i}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}i}^{\prime }$ values for each invader i, and averaged them over the 120 simulations. The approximations for ${\varDelta } {\kappa }_{\text {{P}}i}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}i}^{\prime }$ calculates the actual impact of natural enemies and competition, and then subtract out the expectation in a model without adult spatial structure. Thus, if our approximation of the non-spatial structure case is accurate, then we predict that ${\varDelta } {\kappa }_{\text {{P}}i}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}i}^{\prime }$ should be approximately 0 on average.

Second, we attempted to study our model in a non-invasion context. We could no longer assume that ${\tilde {\lambda }}_{j}^{\prime }(t)=0$ for all species; however, because the community had a fixed number of sites, $\mathbf {E}[{N_{j}(t)}]\tilde {\lambda }_{j}$ needed to sum to 1 across all species. Thus, for all species

$$ {\tilde{\lambda}}_{j}^{\prime}(t) = {\tilde{\lambda}}_{j}^{\prime}(t) - {\sum}_{r=1}^{R}\mathbf{E}\left[{N_{r}(t)}\right] \tilde{\lambda}_{r}^{\prime}(t) $$

(77)

where the summation is over all species. We then defined the following terms:

$$ \begin{array}{@{}rcl@{}} \widetilde{{\varDelta} Y_{j}} &=& \ln \left\{ Y_{j} \right\} + \ln \left\{ I_{j} \right\} - {\sum}_{r=1}^{R}\mathbf{E}\left[{N_{r}(t)}\right]\left( \ln \left\{ Y_{r} \right\} + \ln \left\{ I_{r} \right\} \right).\\ \widetilde{{\varDelta} P_{j}} &=& - \mathbf{E}\left[{ \sum\limits_{y=1}^{L^{2}}p_{j}d_{\text{{p}}-j}\left( |y-x|\right) N_{j}(|y-x|,t)}]\right] \\ &+& {\sum}_{r=1}^{R}\mathbf{E}\left[{N_{r}(t)}\right]\mathbf{E}\left[{ \sum\limits_{y=1}^{L^{2}}p_{r}d_{\text{{p}}-r}\left( |y-x|\right) N_{r}(|y-x|,t) }\right]\\ \widetilde{{\varDelta} \kappa_{\text{{P}}j}}& =& -L^{2} \mathbf{E}\left[{\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ S_{j}(x,t) \right\} \right)}\right] \\ &+& {\sum}_{r=1}^{R}\mathbf{E}\left[{N_{r}(t)}\right] L^{2}\mathbf{E}\left[{ \text{\textbf{cov}}\left( d_{r}(|y-x|), \ln \left\{ S_{j}(x,t) \right\} \right)}\right] \\ \widetilde{{\varDelta} \kappa_{\text{{C}}j}} &\approx& -L^{2}\mathbf{E}[{\text{\textbf{cov}}\left( d_{j}(|y-x|), \ln \left\{ C(y,t) \right\} \right)}] \\ &+& {\sum}_{r=1}^{R}\mathbf{E}\left[{N_{r}(t)}\right] L^{2} \mathbf{E}[{\text{\textbf{cov}}\left( d_{r}(|y-x|), \ln \left\{ C(y,t) \right\} \right)}] \end{array} $$

(78)

Therefore, if spatial structure in the adult community had no effect,

$$ {\tilde{\lambda}}_{j}^{\prime}(t) \approx \widetilde{{\varDelta} Y_{j}} + \widetilde{{\varDelta} P_{j}} + \widetilde{{\varDelta} \kappa_{\text{{P}}j}} + \widetilde{{\varDelta} \kappa_{\text{{C}}j}}. $$

(79)

To test the quality of our approximation, we simulated a shaken community with a large number of sites. We compared this to a community which was governed by the above dynamical system. All of the parameters and initial conditions were identical; thus, if our approximation perfectly represented the simulation, their dynamics would be identical to within demographic stochasticity. Additionally, at each time step, we calculated the covariance between seed density at a site x and both $\ln \left \{ S_{j}(x,t) \right \}$ and $\ln \left \{ C(x,t) \right \}$, along with our estimates for the covariance between those terms at d_j(|y − x|), Eqs. 46 and 65.

Our approximations performed extremely well. Our approximation tended to slightly overestimate the invader growth rate of each species (Fig. 7). The measured values of ${\varDelta } {\kappa }_{\text {{P}}i}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}i}^{\prime }$ tended to be close to 0. For example, in the six communities in Fig. 7, $\overline {{\varDelta } \kappa ^{\prime }_{\text {{P}}}}$ was measured to be between − 0.0016 and − 0.0033, and $\overline {{\varDelta } \kappa ^{\prime }_{\text {{C}}}}$ was between − 0.0020 and − 0.0036 (i.e., indicating that there was this much error). In the non-equilibrium setting, our predictions from Eqs. 78 and 79 tended to match the simulations fairly well (Fig. 8). Additionally, during the simulations, the measured the covariance between seeds at x and both $\ln \left \{ S_{j}(x,t) \right \}$ and $\ln \left \{ C(x,t) \right \}$ closely matched our estimates for those terms (Fig. 8).

Fig. 7

We ran invasion analyses on communities without spatial structure (i.e., using the shaken model) and compared its outcome to invasion analyses run using using Eqs. 78 and 79. Each dot represents the invader growth rate of a particular species. If the approximations were perfect, the dots would be on the one-to-one line. Each graph lists the parameter set used in the title

Fig. 8

a We simulated a community without spatial structure (i.e., using the shaken model) in a large community (L = 300) and present the populations as the solid lines. We also simulated the growth of the community using Eqs. 78 and 79 and present the populations as the dotted lines (corresponding colors were the same species). We used “equal sensitivity” parameters. b At each time step of the simulation, we recorded $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \sum \limits _{y=1}^{L^{2}}p_{i}d_{\text {{p}}-j}\left (|y-x|\right ) N_{j}(|y-x|,t) \right )$, and also used Eq. 46 to estimate what it should be. Each dot represents an estimate at one time point for one species (different colors refer to different species). If our approximation was exact, all dots would be exactly over the one-to-one line. c This graph is similar to b, except comparing $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ C(x,t) \right \} \right )$ to approximation (65). d through f are the same as a through c, except with the “equal dispersal” parameters. g through i are the same as a through c, except with the “random” parameters

Appendix A5: Studying our model with simulations

In this section, we describe how we used computer simulations to analyze our model.

First, simulations suggest that spatial structure made communities less stable (similar results have been shown by Bolker and Pacala (1999) and Durrett and Levin (1994)). In several of our parameter sets, fewer species persisted when there was spatial structure than when there was no spatial structure. For example, in the “equal sensitivity” parameter set, 8 species persisted without spatial structure, but only 6 persisted with spatial structure. In the “equal dispersal” parameter set, 8 species persisted without spatial structure, but only 7 persisted with spatial structure. In the “Random parameter” set, 7 species persisted without spatial structure, and 5 persisted without spatial structure. Additionally, when we ran an invasion analysis on a community with spatial structure, and then compared this to an identical community without spatial structure (i.e., a shaken model), we found that invader growth rates were lower with spatial structure (Fig. 9). In general, differences in invader growth rates can mostly be attributed to the impacts of spatial structure, ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$. For example, in the communities shown in Fig. 9, the ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ terms explained between 67% and 98% of the difference between the simulated community and the non-spatial approximation. The remaining differences are likely caused in part by demographic stochasticity in simulations, but could also be caused by shifts in the relative abundances (which would impact $\overline {{\varDelta } P_{j}}$).

Fig. 9

We ran invasion analyses on communities with spatial structure, and compared its outcome to invasion analyses run using Eqs. 78 and 79 (which does not account for spatial structure). Each dot represents the invader growth rate of a particular species. If a dot is below the line, it indicates that spatial structure reduced that species’ invader growth rates. Each graph lists the parameter set used in the title

The impact of adult aggregation on mortality risk, ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$, appeared to undermine diversity. In all of the our parameter sets, $\overline {{\varDelta } {\kappa }_{\text {{P}}j}^{\prime }}$ was negative (i.e. it was destabilizing). It was always lower than $\overline {{\varDelta } \kappa _{\text {{P}}j}}$, and often about twice the magnitude. ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ also produced fitness differences, which tended to be smaller in magnitude than the fitness differences produced by Δκ_Pj. We also found that ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ was negative for all species in nearly every simulation we ran, i.e. that even the biggest fitness advantage granted by ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ could not make up for its destabilizing effect.

The impact of adult aggregation on competition, ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$, also undermined diversity, though the impact was less strong. In all cases, $\overline {{\varDelta } {\kappa }_{\text {{C}}j}^{\prime }}$ was negative; however, its actual effect was typically small, and usually much smaller than $\overline {{\varDelta } \kappa _{\text {{C}}j}}$ and $\overline {{\varDelta } {\kappa }_{\text {{P}}j}^{\prime }}$. It also tended to produce fitness effects; the fitness effects were always smaller than those of Δκ_Cj (the impact of the parent on competition), but could be bigger than the fitness effects of natural enemies.

We found a strong and consistent correlation between a species’ abundance and its aggregation. Measuring aggregation as $\ln \left \{ g(3)-1 \right \}$ (i.e. the natural log of the pair correlation at 3 sites, or 30m, minus 1), we found that across our parameter sets, frequency and aggregation was approximately

$$ \ln \left\{ g(3)-1 \right\} \approx a + b \mathbf{E}[{N_{j}(t)}] $$

(80)

for some a and b. In most cases, the slope of the line for a particular species was about − 1.2. This varied slightly but was always between − 1.1 and − 1.25 in the parameter range we tested. We are not sure why this occurred, and how much it was driven by our specific parameters. The intercept value, a differed between species. It appeared that dispersal differences had the biggest impact on a, with species with high dispersal being the least aggregated. Additionally, species with strong potential NDD had a lower a, i.e. were less aggregated, though the impact was slightly weaker than the impact of dispersal.

We found that the terms the covariance between seed density and both $\ln \left \{ S_{j}(x,t) \right \}$ and $\ln \left \{ C(x,t) \right \}$, the covariance which contribute to ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$, tend to increase with pair correlation function. For example, in Figs. 10 and 11 we show how they change with the pair correlation at 30m for the three different datasets. The $\ln \left \{ S_{j}(x,t) \right \}$ covariance increases with clustering in all cases, and the $\ln \left \{ C(x,t) \right \}$ covariance increases with clustering in all but 1 case. Both of these covariances contribute negatively to population growth, thus, the fact that they increase shows why $\overline {{\varDelta } \kappa ^{\prime }_{\text {{P}}}}$ and $\overline {{\varDelta } \kappa ^{\prime }_{\text {{C}}}}$ both destabilize coexistence.

Fig. 10

We simulated community dynamics, holding one species’ abundance constant. We tracked how clustered the adult population was at each abundance, and also what $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \sum \limits _{y=1}^{L^{2}}p_{i}d_{\text {{p}}-i}\left (|y-x|\right ) N_{i}(|y-x|,t) \right )$ was at that abundance. Each line represents a different species

Fig. 11

We simulated community dynamics, holding one species’ abundance constant. We tracked how clustered the adult population was at each abundance, and also what $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ C(x,t) \right \} \right )$ was at that abundance. Each line represents a different species

The covariance terms $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ S_{j}(x,t) \right \} \right )$ and $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ C(x,t) \right \} \right )$ differ significantly between species, even if they are at the same clustering (Figs. 10 and 11). This shows how $\overline {{\varDelta } \kappa ^{\prime }_{\text {{P}}}}$ and $\overline {{\varDelta } \kappa ^{\prime }_{\text {{C}}}}$ produce fitness differences. The differences occur because the covariances are driven by other factors as well. For example, species that are more susceptible to natural enemies will have a stronger value of $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ S_{j}(x,t) \right \} \right )$. Similarly, species with low seed dispersal will have a higher value of the covariance $\text {\textbf {cov}}\left (\text {seed density at \textit {x}}, \ln \left \{ S_{j}(x,t) \right \} \right )$, because their seeds are more likely to stay in high-mortality areas.

Finally, we tested if trade-offs would emerge naturally due to community assembly. To do this, we randomly assembled communities of 15 species. Each species was given a potential NDD for distance-responsive enemies (p_A,j) and a dispersal distance parameter (α_j); these values were assigned randomly without replacement. Additionally, potential NDD for distance-responsive enemieS,p_S,j, was assigned as 0.78p_A,j. We calculated the correlation between p_A,j and α_j among the 15 species. We then simulated community dynamics for 1500 time steps, to determine which species persisted. We then calculated the correlation between p_A,j and α_j among species that persisted. We repeated this for 150 communities.

Community assembly tended to create a positive correlation between $p_{j}^{\text {{A}}}$ and α_j, such that species that dispersed farther were more susceptible to enemies (Fig. 12a). Our initial communities began with 0 correlation on average, though this varied with a standard deviation of 0.26. On average, 9.3 species survived in a given random community. The average correlation between p_A,j and α_j among survivors was 0.43 (with a standard deviation of 0.29), it was positive in 92% of simulations, and was higher than the initial community in 97% of simulations.

Fig. 12

a We assembled communities out of 15 species. Each individual had its p_A,j and α_j values assigned randomly from a distribution (without replacement). We then allowed the species to compete, and determine which persisted. We calculated the correlation between p_A,j and α_j for all 15 species, and for those that persisted; each dot represents one such comparison. A positive correlation indicates that species that dispersed farther were more sensitive to their natural enemies (i.e. a trade-off was was occurring). p_A,j values were uniformly distributed from 0.1 to 0.5. α_j values were uniformly distributed from 5.5 to 7. p_S,j values were set to 0.78p_A,j. For other parameters, we used: I_j = 0.27, δ = 0.4, L = 120, l = 10, and Y_j = 7.8. b Here we show one such randomly assembled community. Each circle represents a species that persisted, and each x represents a species that was excluded. In the initial community, the correlation between p_A,j and α_j was − 0.04; in the final community, the correlation was 0.39

Additionally, these results show that while a trade-off would help species coexist, a strict trade-off was not necessary, as it is for the competition-colonization trade-off and many other models of dispersal-mediated coexistence. In fact, a strict trade-off did not emerge in any of our 150 communities. Fig. 12b shows an example of this. Here, 9 of the original 15 species persisted. Among survivors, the species with the third highest dispersal has the third lowest enemy susceptibility, making that species strictly better than 5 of its 8 competitors. However, filtering still produced a correlation of 0.39; this occurred because the species that were susceptible to natural enemies could not survive unless they had reasonably high dispersal, and species that had low dispersal could not survive unless they were tolerant of their natural enemies.

Appendix A6: Measuring terms empirically

In this section, we discuss how various model parameters could be measured empirically.

A6.1 Measuring ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$

The terms ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ quantify how the spatial distribution of adults affects the amount of seedling mortality (due to natural enemies) and seedling competition that species experience. We believe that elements of them can be measured in an equilibrium setting. To do this, one would need to map the distribution of seedlings in many locations and estimate the impact of natural enemies at each location. This was done, for example, by Comita et al. (2010), who tracked seedling survival across 20,000 plots in Barro Colorado Island, Panama. Alternatively, one could grow seedlings in soil collected from many locations (using standard methods for plant-soil feedback studies; (Mangan et al. 2010)). First, the covariance that would produce ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }+{\varDelta } \kappa _{\text {{P}}j}$ can be quantified as the covariance between seedling density and seedling mortality. For species j, one would need to calculate the expected chance that a seedling of species j would die in each of the 20,000 plots, and calculate how this covaries with the actual number of seedlings in each plot, and then divide this by the mean number of seedlings found per plot. It is unclear what would be the best way to estimate survival in a plot without seedlings (such sites could be left out, or perhaps estimated using a regression model). Similarly, the covariance that would produce ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }+{\varDelta } \kappa _{\text {{C}}j}$ is the covariance between the number of species j seedlings in a plot and the total number of seedlings in that plot, divided by the mean number of species j seedlings in a plot.

The above method has a few drawbacks. First, both of these effects calculate the total impact of spatial structure (i.e. including the impact of the seedling’s parent); this is probably preferable, though if the impact of the parent needs to be removed, one could subtract what we expect Δκ_Pj and Δκ_Cj to be (based on our model). Second, the method for calculating ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }+{\varDelta } \kappa _{\text {{P}}j}$ technically measures the spatial structure of seedling survival, rather enemies specifically. It could potentially be flawed if habitat partitioning is occurring, and species are clustered in areas where they have high density-independent survival (Stump and Chesson 2015). Similarly, the method for calculating ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }+{\varDelta } \kappa _{\text {{C}}j}$ is based on the assumption that inter- and intraspecific competition are similar (i.e. that crowding matters, but not which species are causing the crowding); our method will give an inaccurate picture if this assumption is far from true. Last, both measure the covariances in an equilibrium condition, whereas ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }+{\varDelta } \kappa _{\text {{P}}j}$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }+{\varDelta } \kappa _{\text {{C}}j}$ are for invasion situations. This could still be useful, as it would likely indicate the scope of fitness differences. It could also be used to test some of our models predictions, such as that the covariances contributing to ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }+{\varDelta } \kappa _{\text {{P}}j}$ and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }+{\varDelta } \kappa _{\text {{C}}j}$ will be higher for species with low seed dispersal, and for rarer species.

A6.2 Measuring how limited seed dispersal reduces the stabilizing mechanism

In this section, we derive (19) of the main text, which shows how limited seed dispersal reduces the stabilizing effect of natural enemies. We then give an example of how one could measure this, using parameter estimates from Muller-Landau et al. (2008).

The basic principal of Eq. 19 is that we calculate what the stabilizing effect from $\overline {{\varDelta } P_{j}}$ and $\overline {{\varDelta } \kappa _{\text {{P}}j}}$ are, and divided this by what it would be if seeds had 100% dispersal (i.e. $\overline {{\varDelta } P_{j}}$ alone). The term $\overline {{\varDelta } P_{j}}$ is the mean of p_jE[N_j(t)], and $\overline {{\varDelta } \kappa _{\text {{P}}j}}$ is the mean of $-p_{j}\widetilde {d_{j}d_{\text {{P}},j}}\mathbf {E}[{N_{j}(t)}]$. Thus,

$$ \begin{array}{@{}rcl@{}} &&\frac{\overline{{\varDelta} P_{j}} +\overline{{\varDelta} \kappa_{\text{{P}}j}}}{\overline{{\varDelta} P_{j}} } \approx \frac{\overline{p_{j} \mathbf{E}[{N_{j}(t)}]} - \overline{p_{j}\widetilde{d_{j}d_{\text{{P}},j}}\mathbf{E}[{N_{j}(t)}]}}{\overline{p_{j} \mathbf{E}[{N_{j}(t)}]}}\\ &=& 1 - \frac{\overline{p_{j} \mathbf{E}[{N_{j}(t)}]} \overline{\widetilde{d_{j}d_{\text{{P}},j}}} + \text{\textbf{cov}}\left( p_{j} \mathbf{E}[{N_{j}(t)}], \widetilde{d_{j}d_{\text{{P}},j}} \right)}{\overline{p_{j} \mathbf{E}[{N_{j}(t)}]}}\\ &=& 1 - \overline{\widetilde{d_{j}d_{\text{{P}},j}}} - \frac{\text{\textbf{cov}}\left( p_{j} \mathbf{E}[{N_{j}(t)}], \widetilde{d_{j}d_{\text{{P}},j}} \right)}{\overline{p_{j} \mathbf{E}[{N_{j}(t)}]}} \end{array} $$

(81)

$\widetilde {d_{j}d_{\text {{P}},j}}$ is the chance a seed encounters its parent’s natural enemies, and the second term is the mean of that term across species. Of course, this effect is not the same for all species, and the third term accounts for this variation. That third term is not the same as $\frac {1}{R-1}\text {\textbf {cov}}\left ({p_{j}}, {\widetilde {d_{j}d_{\text {{P}},j}}} \right ) + \overline {p_{j}}\text {\textbf {cov}}\left ({\mathbf {E}[{N_{j}(t)}]}, {\widetilde {d_{j}d_{\text {{P}},j}}} \right )$, though it is a good enough approximation for intuition. First, if there is a dispersal-sensitivity trade-off, then the species with the highest p_j (those most susceptible to natural enemies) will have the lowest $\widetilde {d_{j}d_{\text {{P}},j}}$ (i.e. most seeds will escape their parents natural enemies); in this case, the covariance will be negative, and thus the stabilizing mechanism will be strengthened. However, if the species that are the most common have the highest dispersal (which we might expect given the fitness advantage), then this will contribute positively to the covariance, and therefore reduce the stabilizing mechanism.

The third term requires significantly more information to assess. As such, simply quantifying $\overline {\widetilde {d_{j}d_{\text {{P}},j}}}$ will give a first rough estimate, as it will show how much limited seed dispersal weakens the mechanism in the absence of a dispersal-sensitivity trade-off. If a dispersal-sensitivity trade-off is occurring, then $\overline {\widetilde {d_{j}d_{\text {{P}},j}}}$ is an upper-bound of how harmful limited seed dispersal is.

To give an example of how this can be done, we used data from Muller-Landau et al. (2008). They estimated dispersal kernels for 41 species on Barro Colorado Island, Panama. For our work, we focused on the 26 canopy tree species they used (“canopy species” defined in (Comita et al. 2007)). Those values are in Table 4. In a few cases, Muller-Landau et al. (2008) calculated two estimates for α (the dispersal distance parameter in a 2Dt distribution); in these cases, we used data that was fit to seed-fall data rather than fruit.

Table 4 Here we display data from Muller-Landau et al. (2008). We display only canopy trees, as classified by Comita et al. (2007). The parameter α_j is used for calculating a 2Dt distribution

Our model predicts that the dispersal kernel of density-responsive enemies is equal to the seed dispersal kernel. Thus, for density-responsive enemies, $\widetilde {d_{j}d_{\text {{P}},j}}$ is $\widetilde {d_{j}^2}$. So, we first calculated d_j(z) for each species, for a distance between 0 and 8 (we still assumed that trees would take up approximately a 10×10m area). Instead of using Eq. 75, which relied on Manhattan distances, we used the following formula:

$$ d_{j}(z) = \begin{cases} \frac{4}{\pi}{\int}_{0}^{l/2} D(y)dy & z=0\\ \frac{1}{(2\pi z)} {\int}_{l(z-0.5)}^{l(z+0.5)} D(y)dy & z\geq 1\\ \end{cases} $$

(82)

(with D(y) in Eq. 74). For example, Anacardium excelsum has a dispersal parameter α_j = 4.9. Thus,

$$ \begin{array}{@{}rcl@{}} d_{j}(0) &=& 0.2002 \\ d_{j}(1) &=& 0.0747 \\ d_{j}(2) &=& 0.0157 \\ d_{j}(3) &=& 0.0041 \\ d_{j}(4) &=& 0.0015 \\ d_{j}(5) &=& 0.0006 \\ d_{j}(6) &=& 0.0003 \\ d_{j}(7) &=& 0.0002 \\ d_{j}(8) &=& 0.0001 . \end{array} $$

(83)

Next, we squared each term. Finally, we multiplied d_j(z)² by the area between 10(z − 0.5) and 10(z + 0.5) (i.e., 0.79, 6.28, 12.57, 18.85, 25.13, 31.42, 37.70, 43.98, and 50.27m²), and summed those values. Thus, for A. excelsum, it is 0.07.

We found that across canopy trees, $\widetilde {d_{j}d_{\text {{P}},j}}$ varies from 0.002 (Zanthoxylum ekmanii) to 0.193 (Beilschmiedia pendula), with a mean of $\overline {\widetilde {d_{j}d_{\text {{P}},j}}} = 0.057$. This suggests that in the absence of a dispersal-sensitivity trade-off, limited seed dispersal weakens the stabilizing effect of density-responsive enemies by about 6%. Future studies could improve this estimate by quantifying how p_j covaries with $\widetilde {d_{j}d_{\text {{P}},j}}$, and whether more common species have higher $\widetilde {d_{j}d_{\text {{P}},j}}$.

A6.3 Measuring dispersal-susceptibility trade-offs

In this subsection, we discuss how our model could be used to test for dispersal-susceptibility trade-offs.

At a basic level, a dispersal-susceptibility trade-off indicates that species that disperse their seeds farther are more susceptible to natural enemies. This could be done, by calculating each species’ mean dispersal distance, or some similar parameter (e.g. α in a 2Dt distribution, Muller-Landau et al. (2008)), and correlating this with their potential CNDD (measured using techniques discussed in the main text). If the correlation is positive (i.e. species with higher dispersal are more susceptible), it indicates a trade-off is occurring. Indeed, this is essentially what we did in our simulations in Fig. 12. We are not currently aware of anyone who has done this.

Ideally, one would show that this trade-off is equalizing. This is somewhat complicated, as there is no good definition of “equalizing” in a multispecies community (Barabás et al. 2018). A mechanism is called equalizing if it reduces mean fitness-differences between species. In a two-species community, this is simple to calculate, as sp. 1’s fitness is equal to negative sp. 2’s fitness; thus, if a trade-off causes the absolute value of both decrease, it is equalizing. In a community with three or more species, however, this result is more complex.

With that caveat aside, we would suggest that the best way to see if this trade-off is equalizing is to check whether the trade-off affects the variance in mean fitness differences. The term ΔP_j will not contribute much to fitness-differenceS,Δκ_Cj requires a huge number of parameters (and we worry the quantitative values in Eq. 71 may be model-specific), and ${\varDelta } {\kappa }_{\text {{C}}j}^{\prime }$ and ${\varDelta } {\kappa }_{\text {{P}}j}^{\prime }$ cannot be calculated reliably from parameters alone; therefore, we suggest focusing on Δκ_Pj. To do this, one would calculate $p_{j} \widetilde {d_{j}d_{\text {{P}},j}}$ for each species, using methods described above, and determine its variance. Then, one would calculate what the variance of $p_{j} \widetilde {d_{j}d_{\text {{P}},j}}$ is if all species had identical $\widetilde {d_{j}d_{\text {{P}},j}}$. We would recommend setting $\widetilde {d_{j}d_{\text {{P}},j}}$ to what $\widetilde {d_{j}d_{\text {{P}},j}}$ would be if every species had the same mean seed dispersal, although an argument could be made for setting $\widetilde {d_{j}d_{\text {{P}},j}}$ to the community-average value of $\widetilde {d_{j}d_{\text {{P}},j}}$. In either case, if the variance is higher when species have identical $\widetilde {d_{j}d_{\text {{P}},j}}$ values, then the trade-off is equalizing; otherwise it increases fitness differences.