Consequences of spatial patterns for coexistence in species-rich plant communities

Study areas

Nine large forest dynamics plots of areas between 20 and 50 ha were used in the present study (Supplementary Table 1). The forest plots are part of the ForestGEO network⁴ and are situated in Asia and the Americas at locations ranging in latitude from 9.15° N to 45.55° N. Tree species richness among the plots ranges from 36 to 468. All free-standing individuals with diameter at breast height (dbh) ≥1 cm were mapped, size measured and identified. We focused our analysis here on individuals with dbh ≥ 10 cm (resulting in a sample size of 131,582 individuals) and focal species with more than 50 individuals (resulting in 289 species). The 10 cm size threshold excludes most of the saplings and enables comparisons with previous spatial analyses^20,35,47,48. Shrub species were also excluded.

Some of our analyses require estimation of the ratio β_fi/β_ff that describes the relative individual-level competitive effect¹⁸ of individuals of species i on an individual of the focal species f. We used for this purpose phylogenetic distances⁴⁹ based on molecular data, given in Myr, that assume that functional traits are phylogenetically conserved^19,26,27. In this case, close relatives are predicted to compete more strongly or to share more pests than distant relatives²⁶. To obtain consistent measures among forest plots, phylogenetic similarities were scaled between 0 and 1, with conspecifics set to 1, and a similarity of 0 was assumed for a phylogenetic distance of 1,200 Myr, which was somewhat larger than the maximal observed distance (1,059 Myr). This was necessary to avoid discounting crowding effects from the most distantly related neighbours²⁶.

Observed spatial patterns at species-rich forests

Figure 1 and Supplementary Data Table 1 show the intraspecific variation in our three crowding indices n_kff, n_kfh and n_kfβ that can be approximated by gamma distributions. To assess how well the gamma distribution described the observed distribution, we used an error index defined as the sum of the absolute differences of the two cumulative distributions divided by the number of bins (spanning the two distributions). The maximal value of the error index is one, and a smaller value indicates a better fit.

Equations (6, 8 and 9) relate the measures of the emerging spatial patterns (that is, k_ff, k_fh and B_f) to macroscale properties and conditions for species coexistence. Even though our multiscale model (equation (7)) is simplified, it allows for a direct comparison with the emerging patterns in our nine fully stem-mapped forest plots. We estimate the key quantities of equations (8) and (9) directly from the forest plot data (Fig. 4), with the exception of the carrying capacities K_f, which were indirectly estimated from the observed species abundances (assuming approximate equilibrium; equation (8) and Supplementary Data Table 1). This allowed us to estimate the feasibility index µ_f (equation (9)). Because statistical analyses with individual-based neighbourhood models^19,26 based on neighbourhood crowding indices have shown that the performance of trees depends on their neighbours for R between 10 and 15 m, we estimate all measures of spatial neighbourhood patterns with a neighbourhood radius of R = 10 m. Analyses with R = 15 or R = 20 gave similar results.

The spatial multispecies model and equilibrium

We use a general macroscale model to describe the dynamics of a community of S species:

where r_f is the mean number of recruits per adult of species f within time step Δt, s_f is a density-independent background survival rate of species f and the α_fi are the population-level interaction coefficients, yielding α_ff = c γ_ff k_ff β_ff and α_fi = c γ_fβ k_fh β_ff B_f (equation (6)). The β_fi are the assumed individual-level interaction coefficients between individuals of species i and f; k_ff = K_ff(R) / π R² and k_fh = K_fh(R) / π R² measure intraspecific clustering and interspecific segregation, respectively, with K_ff(R) being the univariate K function for species f and K_fh(R) the bivariate K function describing the pattern of all heterospecifics ‘h’ around individuals of species f. A is the area of the observation window.

Following equation (5), B_f can be estimated as

and is the weighted average of the relative individual-level interaction coefficients β_fi/β_ff between species i and the focal species f, weighted by the mean number of individuals of species i in the neighbourhoods of the individuals of the focal species (that is, c k_fi N_i(t)). For competitive interactions, B_f ranges between zero and one; B_f = 1 indicates that heterospecific and conspecific neighbours compete equally, and smaller values of B_f indicate reduced competition with heterospecific neighbours. The denominator can be rewritten in terms of segregation k_fh to all heterospecifics and the total number of heterospecifics ∑_i≠f N_i(t).

The analytical expression of the equilibrium (equation (8)) relies on the assumption that the values of B_f are approximately constant in time. This assumption may not apply in our model during the initial burn-in phase of the simulations if the β_fi/β_ff show large intraspecific variability (Supplementary Text and Figs. 1–5). The underlying mechanism is the central niche effect introduced by Stump⁴⁵ where a species has reduced average fitness if it has high niche overlap with competitors.

Finally, the factors γ_ff = ln(1 + b_ff β_ff) (b_ff β_ff)⁻¹ and γ_fβ = ln(1+ b_fβ β_ff) (b_fβ β_ff)⁻¹ describe the influence of the variance-to-mean ratios b_ff and b_fβ of the gamma distribution of the crowding indices n_kff and n_kfβ, respectively. For high survival rates during one time step (for example, >85%), the values of γ_ff and γ_fβ are close to one; in this case the exponential function in equation (1a) can be approximated by its linear expansion and γ_ff = γ_fβ = 1.

In equilibrium we have (N_f(t + Δt) ‒ N_f(t))/Δt = 0, which leads, with equation (7), to:

with (K_f = – {mathrm{ln}}left( {frac{{1 – r_f}}{{s_f}}} right) left( alpha _{ff} right)^{-1}) and the total number of individuals being (J^{ast} = sum _iN_i^{ast}). Rewriting equation (13) yields (frac{K_f}{J^{ast}} = left( frac{N_f^{ast}}{J^{ast}} right) left(1- frac{alpha_{f{mathrm{h}}}}{alpha_{ff}}right) + frac{alpha_{f{mathrm{h}}}}{alpha_{ff}}). For α_fh/α_ff < 1 we therefore find K_f < J^*, which indicates that a multispecies forest would host more individuals than a monoculture. To estimate J^* we sum equation (13) over all species i and find (J^{ast} = mathop {sum }limits_i frac{{K_i}}{{1 – alpha _{i{mathrm{h}}}/alpha _{ii}}} – J^{ast} mathop {sum }limits_i frac{{alpha _{i{mathrm{h}}}/alpha _{ii}}}{{1 – alpha _{i{mathrm{h}}}/alpha _{ii}}}). Therefore, we obtain

with (m_K = frac{1}{S}mathop {sum }limits_i^{,} frac{{K_i}}{{1 – alpha _{i{mathrm{h}}}/alpha _{ii}}}) and (m_alpha = frac{1}{S}mathop {sum }limits_i^{,} frac{{alpha _{i{mathrm{h}}}/alpha _{ii}}}{{1 – alpha _{i{mathrm{h}}}/alpha _{ii}}}) being averages over the S species of the community.

All species have positive abundances at equilibrium if the two conditions µ_f = α_fhJ^*/α_ffK_f < 1 and α_fh/α_ff < 1 are met (see equation (9)). We now show that the chance that these conditions are satisfied for all species is larger if the values of µ_f show little intraspecific variability. To understand this, we assume that the µ_f can be approximated by their mean (bar mu). In this case (J^{ast} /bar mu) is also approximately constant and we can replace K_i in equation (14) by ((J^{ast} /bar mu ) (alpha _{i{mathrm{h}}}/alpha _{ii})^{-1}) and obtain

where S is the number of species in the community, and therefore

Thus, in the case of a perfect interspecific balance in µ_f we always have a feasible equilibrium if α_fh/α_fh < 1, and species can go extinct only if the intraspecific variability in µ_f becomes too large. The smaller the mean value of µ_f, the more variability in µ_f is allowed. Equation (16) shows that (bar mu) is smaller if the number S of species in the community is smaller and/or if the mean value of m_α is smaller.

Equation (16) also suggests that communities with more species need to show stronger species equivalence in µ_f because the term S m_α(1 + S m_α)⁻¹ approaches a value of one for a large number of species S. This finding mirrors the results of analyses of Lotka–Volterra models with random interaction matrices¹¹ that showed that the larger the number of species S, the more difficult it becomes to generate a feasible community.

Invasion criterion

Using the population-level interaction coefficients (equation (6)) in the macroscale model, we now derive conditions for coexistence based on the invasion criterion^7,8 for a species m. The growth rate of an invading species m with low density M(t) into the equilibrium community of all other S – 1 species N_i(t) should be positive; thus, with equation (7), we have

Considering that (J_m^{ast} = mathop {sum}nolimits_{i = 1}^{S – 1} {N_i^{ast}}) and (alpha_{mm} M(t) ll alpha_{f{mathrm{h}}} J_m^{ast}) (that is, the invading species m is at low abundance and does not show strong clustering) we find (- ln left( frac{1-r_m}{s_m} right) > alpha_{m{mathrm{h}}} J_m^ ast), and by dividing by α_mm we obtain the invasion condition

which is basically identical to the condition for feasibility (equation (9)), but here the community size (J_m^ ast) of the reduced community appears instead of the equilibrium community size J^* of all species, including species m. Thus, a new species m is more likely to invade if it has a high value of the carrying capacity K_m and if it more strongly reduces heterospecific interactions relative to conspecific interaction (that is, α_mh/α_mm is smaller). However, if the species is too efficient (that is, has too large a capacity K_m and/or too low an α_mh/α_mm) it may increase the value of J^* too much (equation (14)), thereby causing the extinction of the weakest species with the highest values of µ_f (that is, a too-low value of K_m and a too-high value of α_fh/α_ff). Equation (18) also suggest that an equilibrium with µ_f > 1 and α_mh/α_mm > 1 will be unstable.

Fitness–density covariance

To place our new spatial coexistence mechanism in the context of existing coexistence theory, we apply scale transition theory³⁴ to our model version where spatial effects are the only potential coexistence mechanism (that is, all species have the same parameters and all individuals compete equally; β_fi/β_ff = 1, B_f = 1).

Following equation (1a), the expected fitness of an individual k of a focal species f (that is, its expected contribution to the population after some defined interval of time Δt) in the macroscale model (equation (7)) is

where W_k = n_kff + ∑_i≠fn_kfi is the fitness factor of individual k, f(W_k) = exp(−β_ff  W_k) is the fitness function and n_kff and ∑_i≠f n_kfi are the number of conspecific and heterospecific neighbours, respectively, of individual k within distance R. The spatial average of the fitness factor over the entire plot is

where c = πR² / A, and J(t) = ∑_iN_i(t) is the total number of individuals in the plot. Given that J(t) converges very quickly into equilibrium J^* (Extended Data Fig. 5 and Supplementary Figs. 1 and 2), we find for the spatial average fitness (bar lambda _f = 1).

The average individual fitness (tilde lambda _f(t)) of a focal species f is the average of λ_k,f over all individuals k of species f and can be estimated for the macroscale model (equations (6 and 7)) as (tilde lambda _fleft( t right) = N_fleft( {t + {Delta}t} right)/N_fleft( t right)). A key ingredient of scale transition theory³⁴ is that the fitness–density covariance is given by ({mathrm{cov}} = tilde lambda _fleft( t right) – bar lambda _f). With equation (3) and γ_ff ≈ γ_fβ ≈ 1 and (bar n_{fbeta } = bar n_{f{mathrm{h}}}) we find

where the mean of the crowding indices is given by (bar n_{ff}(N_f) = ck_{ff}(N_f)N_f) and (bar n_{f{mathrm{h}}}(N_f) = ck_{f{mathrm{h}}}(J^{ast} – N_f)) (equation (4)). Therefore, if clustering k_ff and segregation k_fh are independent from abundance N_f, more abundant species have more neighbours, since

Thus, a positive fitness–density covariance in our model means that individuals of a common species are more likely to be near more trees in total.

Extended Data Fig. 7 shows the quantities (bar n_{ff} + bar n_{f{mathrm{h}}}), (bar n_{ff}), (bar n_{f{mathrm{h}}}) and (tilde lambda _f – bar lambda _f) plotted over abundance N_t for data generated by our spatially explicit simulation model for the scenarios of stable and unstable dynamics (Extended Data Fig. 5b,c). Indeed, the stable simulations show a positive fitness–density covariance, however, there is no such trend for the dynamics of the unstable community (Extended Data Fig. 7g,h).

Spatial patterns will act as positive fitness–density covariance if, when a species becomes rare, areas of high conspecific crowding have fewer competitors. We tested this for the data generated by our simulation model and for the nine forest plots (Extended Data Fig. 8). We could estimate for each focal species f the covariance between the number of conspecific neighbours (that is, n_kff) and the total number of neighbours (that is, n_kff + n_kfh) and demand that the covariance should be mostly positive and larger for more abundant species. However, since the quantity n_kff appears in this test on both sides, a positive covariance can be expected. To compensate for this artefact, we instead used the covariance between the local dominance of conspecifics in the neighbourhood of individuals k (that is, d_kff = n_kff (n_kff + n_kfh)⁻¹) and total number of neighbours (that is, n_kff + n_kfh) (Extended Data Fig. 8).

Spatially explicit simulation model

The model is a spatially explicit and stochastic implementation of the spatial multispecies model (equation (7)), similar to that of May et al.^35,37 and Detto and Muller-Landau¹⁷, and simulates the dynamics of a community of S tree species in a given plot of a homogeneous environment (for example, 50 ha) in 5 yr time steps adapted to the ForestGEO census interval (Extended Data Fig. 5 and Supplementary Figs. 1 and 2). Only reproductive (adult) trees are considered, but size differences between them are not considered. During a given time step the model first simulates stochastic recruitment of reproductive trees and placement of recruits, and second, stochastic survival of adults that depends on the neighbourhood crowding indices for conspecifics (n_kff) and heterospecifics (n_kfβ) (but excluding recruits). In the next time step, the recruits count as reproductive adults and are subject to mortality. No immigration from a metacommunity is considered. To avoid edge effects, torus geometry is assumed.

The survival probability of an adult k of species f is given by (s_f exp left(-beta_{ff} left( n_{kff}+n_{kfbeta}right)right)) (equation (1a)). The two neighbourhood indices n_kff and n_kfβ describe the competitive neighbourhood of the focal individual k and sum up all conspecific and heterospecific neighbours, respectively, within distance R, but weight them with the relative individual-level interaction coefficients β_fi/β_ff (refs. ^19,21,26).

Each individual produces on average r_f recruits, and their locations are determined by a type of Thomas process²⁸ to obtain clustering. To this end, the spatial position of the recruits is determined by two independent mechanisms. First, a proportion 1 – p_d of recruits is placed stochastically around randomly selected conspecific adults by using a two-dimensional kernel function (here a Gaussian with variance σ²). This is the most common way to generate species clustering in spatially explicit models^{17,18,35,36,37,38,39}. Specifically, we first randomly select one parent for each of these recruits among the conspecific adults and then determine the position of the recruit by sampling from the kernel. Second, the remaining proportion p_d of recruits is distributed in the same way around randomly placed cluster centres that are located independently of conspecific adults. This mode mimics spatial clustering of recruits independent of the parent locations⁴² in a simple way, such as contagious seed dispersal by animals⁵⁰ or forest gaps that may imprint clumped distributions of recruits of pioneer species⁴⁰. For each species we assume a density λ_fc of randomly distributed cluster centres, which have, at each time step, a probability p_fp of changing location. For each of these recruits, we first randomly select one cluster centre among the cluster centres of the corresponding species and then determine the position of the recruit by sampling from the kernel. For the simulation shown in Extended Data Fig. 5a, the recruits were located at random positions within the plot.

Parameterization of the simulation model

Extended Data Fig. 5 shows simulations of the individual-based model conducted in a 200 ha area containing approximately 83,000 trees with, initially, 80 species. There was no immigration. The model parameters were the same for all species, and all species followed exactly the same model rules. We selected β_fi = β_ff to obtain no differences in con- and heterospecific interactions and s_f = 1 (no background mortality), and we adjusted the parameters β_ff = 0.0075 and r_f = 0.1 to yield tree densities (415 ha⁻¹) and an overall 5 yr mortality rate (10%) similar to those of trees with dbh ≥ 10 cm in the BCI plot⁵¹.

The Gaussian kernel used to place recruits around conspecific adults or around random cluster centres had a parameter σ = 10 m. There were 40 random cluster centres in total for each species that had a probability of p_fp = 0.3 of changing location within one census interval. The only difference between the simulation shown in Extended Data Fig. 5b and the one shown in Extended Data Fig. 5c is that in the former, we used a proportion p_d = 0.05 of recruits to be placed around randomly distributed cluster centres (that is, 95% of the recruits were placed close to their parents), but in the latter, we selected p_d = 0.95 (that is, 95% of the recruits were placed around randomly distributed cluster centres). In our simulations, on average, one of these cluster centres received four recruits per time step, which were scattered within a radius of approximately 30 m, and received approximately 13 recruits during its lifetime (at each time step it had a probability of 0.3 of changing location). In contrast, in Extended Data Fig. 5a recruits were placed at random locations within the plot.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com