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

Ecology cannot yet fully explain why so many tree species coexist in natural communities such as tropical forests. A major difficulty is linking individual-level processes to community dynamics. We propose a combination of tree spatial data, spatial statistics and dynamical theory to reveal the relationship between spatial patterns and population-level interaction coefficients and their consequences for multispecies dynamics and coexistence. Here we show that the emerging population-level interaction coefficients have, for a broad range of circumstances, a simpler structure than their individual-level counterparts, which allows for an analytical treatment of equilibrium and stability conditions. Mechanisms such as animal seed dispersal, which result in clustering of recruits that is decoupled from parent locations, lead to a rare-species advantage and coexistence of otherwise neutral competitors. Linking spatial statistics with theories of community dynamics offers new avenues for explaining species coexistence and calls for rethinking community ecology through a spatial lens.

U nderstanding the mechanisms that maintain high species diversity in plant communities such as tropical forests has long challenged ecologists 1 and has stimulated major efforts in field and theoretical ecology [2][3][4][5] . However, despite a multitude of coexistence mechanisms that have been proposed 6 and recent advances in coexistence theories [7][8][9][10][11][12][13][14][15][16][17][18] , this fundamental question has not been fully resolved 8,9,14 . For example, theoretical models indicate that stable coexistence is difficult to reach in large communities 11,13 . We argue that consideration of spatial patterns of plant individuals, such as intraspecific clustering and interspecific segregation, may allow for a better understanding of mechanisms of coexistence in species-rich communities 17 .
Although many studies suggest that spatial patterns and neighbourhood effects may play an important role in diversity maintenance [17][18][19][20][21] , the integration of spatial patterns into coexistence theories of species-rich communities is difficult. A major difficulty is linking spatial processes at the individual level to community dynamics. One reason for this is a scale mismatch. The analytical models that form the basis of most coexistence theories 7,8,[11][12][13][14][15][16]22 have state variables that operate at the macroscale (that is, the population or community-level abundances), use parameters that describe average 'mesoscale' properties of the individuals (such as population-level interaction coefficients and demographic rates) and often rely on 'mean-field' approximations 18,23 where spatial patterns are neglected.
However, spatial patterns and population-level interaction coefficients emerge at the mesoscale from the microscale behaviour of individuals and their interactions with other individuals and the environment. Therefore, studying the impact of spatial patterns on species coexistence requires multiscale approaches such as spatial moment equations 18,23,24 that incorporate pattern-forming processes operating at the level of individuals and translate these into population and community dynamics.
We propose here such a multiscale approach. To this end, we first derive population-level interaction coefficients α fi from individual-level interaction coefficients β fi and neighbourhood crowding indices 19,21 that are commonly used to describe interactions among tree individuals at the microscale, and then incorporate the emerging coefficients α fi into analytical macroscale models. Our approach is based on separation of timescales (adiabatic approximation 25 ), given that mesoscale spatial patterns usually build up quickly and approach a quasi-steady state whereas the macroscale state variables (for example, abundances) change slowly 23 . Therefore, we do not need to describe the dynamics of the quick mesoscale patterns explicitly (as, for example, is done in approaches based on moment equations 18,23,24 ) but concentrate instead on spatial patterns that transport the critical information from the microscale into macroscale models. This approach requires information on mesoscale spatial patterns that can be obtained from fully mapped forest plots

Consequences of spatial patterns for coexistence in species-rich plant communities
such as those of the Forest Global Earth Observatory (ForestGEO) network 4 .
More specifically, we (1) derive species-level interaction coefficients from individual-level interactions using empirical information on spatial patterns in nine ForestGEO megaplots 4 , (2) integrate the resulting species-level interaction coefficients into analytical macroscale multispecies models and (3) study their consequences for multispecies dynamics and coexistence.

Results and discussion
Species-level interaction coefficients. We first derive species-level interaction coefficients from individual-level neighbourhood crowding indices 19,21,26 that quantify how the performance of a focal individual depends on interactions with its neighbours. To this end, we describe the survival rate of a focal individual k of species f in dependence on the local number of neighbours as where s f is a density-independent background survival rate of species f; the crowding indices n kff , n kfi and n kfh are the number of conspecifics, neighbours of species i and heterospecific neighbours within distance R of a focal individual k, respectively; the subscript 'h' indicates all heterospecifics together (that is, n kfh = ∑ i≠f n kfi ) and the crowding index n kfβ weights each heterospecific neighbour by its relative competition strength β fi /β ff (equation (1a)), with β fi being the individual-level interaction coefficients between species f and i ( Fig. 1a-c). The corresponding population-level survival rate is given by where N i (t) is the abundance of species i at time step t and α fi is the population-level interaction coefficient between species f and i.
To estimate surv f we average the survival rates s kf (equation (1a)) of all individuals k of species f: where p x and p y are the distributions of the crowding indices x = n kff and y = n kfβ for individuals of species f, respectively.
To determine the distributions p x and p y , we analysed forest inventory data from nine 20-50 ha forest dynamics plots (Supplementary  Table 1) in the ForestGEO network 4 . We used phylogenetic similarity between tree species as a surrogate for the relative competition strength β fi /β ff because it is available for the species in the nine plots (Methods). This is an established approach in species-rich communities 19,26,27 to approximate niche differences in the absence of other data.
The number of con-and heterospecific neighbours and the heterospecific interaction index n kfβ vary widely among conspecifics and can be described by gamma distributions (Fig. 1 and Extended Data Figs. 1 and 2). Detailed analysis of the empirical crowding indices reveals additional relationships that are relevant for our subsequent analysis. First, we find that the crowding indices n kff and n kfβ are not, or are only weakly, correlated for a given species f (Extended Data Fig. 3a). Second, we find for trees of a given species f high correlations between the two crowding indices n kfh and n kfβ (Extended Data Fig. 3b) with a common factor B f (that is, n kfβ ≈ B f n kfh ). This result suggests operation of diffuse neighbourhood competition, in which the competition strength of heterospecifics is on average a factor B f lower than that of conspecifics.
The integral of equation (2) can be solved analytically for independent gamma distributions p x and p y and yields where n ff and n fβ are the average values of the crowding indices n kff and n kfβ , respectively, and γ ff and γ fβ contain the variance-to-mean ratios of the gamma distributions p x and p y , respectively, but in our case have values close to one (Methods). The last step in deriving pairwise population-level interaction coefficients is to relate the averages of the different crowding indices to the macroscale population abundances N f (t). We accomplish this by taking advantage of connections between crowding indices and the summary functions of spatial point process theory 21 . The mean of the crowding index n kff (that is, the mean number of further conspecific neighbours within distance R) is proportional to Ripley's K, a well-known quantity in point process theory 28,29 : where K ff (R) is the univariate K function for species f and A is the area of the observation window. The K function describes the spatial pattern of conspecifics within a neighbourhood distance R, indicating clustering if K ff (R) > πR 2 , a random pattern if K ff (R) = πR 2 and regularity if K ff (R) < πR 2 . In the following, we use the normalized K function k ff (R) = K ff (R) /πR 2 to quantify the spatial neighbourhood patterns, and therefore n ff = k ff (R) πR 2 A N f (t). Analogously, the mean number of heterospecific neighbours is given byn where the bivariate normalized K function k fh (R) indicates segregation to heterospecifics (subscript 'h') within distance R if k fh (R) < 1. Independent placement occurs if k fh (R) = 1, and attraction if k fh (R) > 1. Motivated by the finding n kfß ≈ B f n kfh (Extended Data Fig. 3b) we rewrite the mean crowding index n fβ as where the point process summary function B f indicates how much the competition strength of one heterospecific neighbour differs on average from that of one conspecific neighbour. The values of B f depend mainly on the individual-level interaction coefficients β fi but also on the spatial pattern of the different species and their relative abundances (equation (12)). Inserting the expressions for n ff and n fβ into equation (3) and comparing with equation (1b) leads to our first main result, the analytical expressions of the population-level interaction coefficients: with scaling constant c = πR 2 /A. Notably, equation (6b) indicates that the emerging population-level interaction coefficients α fi are the same for all heterospecifics (that is, α fi = α fh for i ≠ f). Thus, even if the individual-level interactions coefficients β fi differ among species pairs, the emerging population-level interaction coefficients α fi have a substantially simpler structure. This phenomenon is an example of simplicity emerging from complex species interactions 30 and is likely to occur only in species-rich communities 9,12,31 . The population-level interaction coefficients α fi depend on several factors that can influence the macroscale balance between intra-and interspecific competition: (1) intraspecific clustering k ff and interspecific segregation k fh (Fig. 2a,b), (2) the relative competition strength B f of one heterospecific neighbour (Fig. 2c) and (3) the shape of the response of survival to crowding (γ ff and γ fβ , which contain the variance-to-mean ratios of the distribution of the crowding indices). Note that absence of spatial patterns (that is, k ff = 1 and k fh = 1) and a linear approximation of equation (3) lead to γ ff = γ fβ = 1 and direct proportionality α fi = cβ fi of the individualand population-level interaction coefficients, as assumed by Lotka-Volterra models.
The information on spatial patterns extracted from the inventory data of our nine forests allows us to estimate the relative population-level interaction coefficients α fi /α ff for all pairs of species i and f ( Fig. 3 and Extended Data Fig. 4). For example, at BCI, the values of α fi /α ff differ substantially from the corresponding individual-level coefficients β fi /β ff (Fig. 3a,c), and for 83% of all species pairs, we find α fi /α ff < β fi /β ff . Thus, the mesoscale spatial patterns can reduce, at the population level, the strength of heterospecific interactions relative to conspecific interactions by 'diluting' encounters with heterospecific neighbours relative to conspecific neighbours. Spatial patterns therefore have a strong potential to alter the outcome of deterministic individual-level interactions.
Conditions for coexistence in the multiscale model. To study the consequences of the emerging spatial patterns for community dynamics and coexistence we insert the population-level interaction coefficients α fi (equation (6)) into a simple macroscale model In this model we assume that survival is governed by neighbourhood competition with α fi = α fh , and the number of recruits of species f during a time step Δt is given by r f N f (t), where r f is the per capita reproduction rate of species f.
The carrying capacity of species f (that is, the equilibrium of equation (7) with N i (t) = 0 for all i ≠ f) is given by . Note that our theory also applies, after redefinition of the carrying capacity, to alternative macroscale models (Supplementary Table 3). From equation (7) we find abundance of species f in equilibrium and J * the equilibrium community size (that is, J * = ∑ i N * i ; see also equation (14)). This leads, under the assumption that the population-level interaction coeffi- that is positive if denominator and numerator are both positive or both negative. However, the invasion criterion (equation (18)) is only fulfilled if both are positive. In this case, equation (8) suggests two different ways a species can go extinct. First, the denominator indicates that a species with strong clustering k ff will show a small equilibrium abundance since in this case α ff ≫ α fh (equation (6)). Large values of k ff can be expected for species of low abundance under dispersal limitation, where recruitment happens close to conspecific adults. Second, the numerator of equation (8) indicates a positive abundance of species f if α ff K f > α fh J * and α fh /α ff < 1. Therefore, we introduce a new feasibility index that indicates a positive abundance if µ f < 1 given that heterospecific interactions at the population level are weaker than conspecific interactions (that is, α fh /α ff < 1). The invasion criterion 7,8 that tests whether a species with low abundance can invade the equilibrium community of all other species turned out to be basically the same as the feasibility condition (equation (9)) if the invading species does not show strong clustering (Methods and equation (18)). Note that we did not assume Allee effects 8 . Further analysis that considers the dependency of J * on the values of K f and α fh /α ff shows that the values of µ f must be similar for all species f to fulfil the condition µ f < 1, and that µ f can show larger interspecific variability if the species richness S is smaller and/or if the mean of (14) and (15)).
The feasibility index µ f therefore governs species assembly by determining the subset of species of a larger species pool that can persist 13 , but any addition of a new species changes µ f and may lead to reassembly of the community. Using the observed abundances in the forest plots (and assuming equilibrium) allows us to test our theory. We can estimate from the observed abundances the carrying capacities K f and therefore also the indices µ f for all focal species of our nine plots (Fig. 4). For 282 of our 289 focal species, we found µ f < 1 and α fh /α ff < 1, which means that the two conditions for stable coexistence are indeed satisfied for nearly all species. However, this is not a given, as shown by the seven species from BCI with µ f > 1 and α fh /α ff > 1. Thus, our theory is compatible with the observed coexistence of most species at our nine forest plots if the assumption of approximately constant population-level interaction coefficients holds.
In addition, we get information on the typical values of µ f and α fh /α ff that allows for insight into the stability of the communities. In agreement with the predictions of our theory, we find that µ f tends to be smaller if α fh /α ff is smaller (Fig. 4c). Furthermore, the values of µ f were, for most species, larger than the expectation of µ f for the corresponding communities without interspecific variability in µ f (equation (10a)) but with the same number of species and the same mean values of α fh /α ff (Fig. 4c), but all of the values were relatively close to the critical value of 1 (the median of all 289 species was 0.938; Fig. 4a).
Consequences of spatial patterns for coexistence. Our theory predicts that coexistence requires, in the limit of high species richness, that species approach functional identity with respect to the feasibility index µ f (Methods). This resembles neutral theory 2,32 , but our theory allows for trade-offs among demographic parameters and emerging spatial patterns to reach this equivalence (equation (9)). To study the consequences of spatial patterns for coexistence, we analysed a symmetric 33 version of our model where all species have the same parameters and follow the same stochastic rules and where all individuals compete identically (leading to B f = 1). Thus, we eliminate any potential coexistence mechanism other than that resulting from spatial patterns.
If the mesoscale patterns k ff and k fh converge to a stochastic equilibrium, we find that the feasibility and invasion criteria are always fulfilled if α fh /α ff < 1 since Equations (10a) and (10b) follow from equations (16) and (18), respectively, if α fh /α ff is the same for all species f. Thus, the spatial patterns that emerge at the mesoscale from the individual-level interactions can stabilize if α fh /α ff < 1. The underlying mechanism is a positive fitness-density covariance 34  Intraspecific clustering is indicated by k ff > 1 and interspecific segregation by k fh < 1, and the less a heterospecific neighbour competes on average relative to a conspecific neighbour, the more B f decreases. The neighbourhood radius used was R = 10 m. For the analysis, we used all individuals with dbh ≥ 10 cm and included focal species with more than 50 individuals. For forest plot names, see Supplementary Data Table 1.
To reveal the conditions that can lead to coexistence in a spatially explicit context, we use a spatially explicit and individual-based 35 implementation of the symmetric version of the multiscale model (equation (7) and Methods). Models of this type are able to produce realistic spatial patterns consistent with mapped species distributions of large forest plots 17,35 . While our analytical approach in equation (7) only allows us to make simplified assumptions about the spatial component of the recruitment process, the simulation model allows us to explore the role of the spatial component of recruitment in more detail.
Indeed, the way recruits were placed was critical for coexistence. Randomly placed recruits produced unstable dynamics (Extended Data Fig. 5a) characterized by regularity (mean of k ff = 0.92) and segregation (mean of k fh = 0.92), both caused by competition 18 , and the instability was caused by con-and heterospecifics competing equally at the population level (that is, α fh /α ff ≈ 1) (Extended Data Fig. 5d,g). When we followed the common approach of placing recruits with a kernel around conspecific adults 17,18,[35][36][37][38][39] to mimic dispersal limitation, we again found unstable dynamics (Extended Data Fig. 5b), despite intraspecific clustering and interspecific segregation (that is, α fh /α ff < 1; Extended Data Fig. 5n). The reason for the instability was high clustering of rare species 24,35 (Extended Data Fig. 6) that completely negated the potentially positive effects of α fh /α ff < 1.
In contrast, community dynamics can be stabilized if recruits are placed in small clusters but independent of the location of conspecific adults (Extended Data Fig. 5c). With this mechanism we mimic canopy gaps 40 , animal seed dispersal 41 or other mechanisms that can generate clustering independent of parent locations, as found at BCI 42 . Decoupling clustering from the parent locations does not lead to the negative relationship between clustering and abundance, and all measures of spatial patterns converged quickly into quasi-equilibrium (Extended Data Fig. 5f,i,o).
The simulation data reveal the spatial coexistence mechanism underlying the positive fitness-density covariance 34 (Extended Data Figs. 7 and 8). We find that the emerging spatial patterns lead to a situation where individuals of a common species are more likely to be near more neighbours and tend therefore to experience stronger competition (Extended Data Fig. 7). While the number of heterospecific neighbours remains approximately constant, the number of conspecific neighbours decreases with decreasing abundance if clustering does not change with abundance (equation (4a)). However,  Table 1.  Table 1.
if clustering increases with decreasing abundance, the rare-species advantage is weakened and the dynamics become unstable 24 . The data of several ForestGEO forest plots were compatible with a positive fitness-density covariance (Extended Data Fig. 8f-i) as they show that, when a species becomes rare, areas of higher conspecific crowding tend to have fewer total competitors. Comparison of the results of the stable versus unstable simulation showed that even relatively weak tendencies in this relationship are sufficient to stabilize the dynamics (Extended Data Figs. 8a,b). However, this was not the case for three temperate forest plots where the power-law clustering-abundance relationship showed exponents of b < −0.5 (that is, species tended to have high clustering at low abundances), but the other plots showed b > −0.5 (Extended Data Fig. 8n-p).
The apparent contradiction with previous theoretical studies 18,24,35,37 where intraspecific clustering and interspecific segregation could not stabilize community dynamics thus arises as a consequence of the assumption of placing recruits close to their parents. This finding has important consequences for ecological theory because it shows, in contrast to the prevalent view 36,43,44 , that spatial patterns alone can lead to coexistence of multiple species. This is even more important since the specific spatial patterns required for this coexistence mechanism also exist in real forests.

Conclusions
Understanding the mechanisms that maintain high species diversity in communities such as tropical forests is at the core of ecological theory, but these mechanisms are not yet fully resolved.
Here, we argue that spatial patterns may play an important role in species coexistence of high diversity plant communities 17 . To test this hypothesis, we introduced a multiscale framework that reveals how pattern-forming processes operating at the level of individuals translate into mesoscale spatial patterns and how those patterns influence macroscale community dynamics.
We showed that the population-level interaction coefficients α fi can have, for a broad range of common circumstances, a simpler structure than the underlying individual-level interaction coefficients β fi . This simplicity, which emerged from spatially explicit species interactions 30 , allowed for an analytical treatment of equilibrium, feasibility and invasion conditions of the corresponding macroscale models (equations (8)(9)(10)). Inserting the emerging α fi coefficients into macroscale community models (for example, equation (7); Supplementary Table 3) should, in principle, allow us to take advantage of macroscale theory 9,[11][12][13]45 . However, our results also indicate that the population-level interaction coefficients may not be temporally constant as commonly assumed but depend on spatial patterns that may change with abundance. This is especially likely if recruitment is mainly located close to the parents.
It is also possible to expand our framework to take into account more detailed neighbourhood crowding indices that consider not only the number of neighboured trees of a given species but also their distance and size 19,21,26 . This requires redefinition of the quantities k ff and k fh that describe intraspecific clustering and interspecific segregation, respectively, but does not change the overall structure of our equations. A special strength of our approach is that the population-level interaction coefficients contain measures of spatial neighbourhood patterns that can be directly estimated from fully mapped forest plots 4 . Together with additional information, this may allow for estimating network structures as well as stability of the whole community.
Our analysis revealed that communities of competing species can show a stable mode where the mesoscale patterns converge quickly into quasi-equilibrium and an unstable mode where negative relationships between species clustering and abundance emerge (Extended Data Figs. 5 and 6). The two modes are governed by the way species clustering is generated: the well-known unstable mode is related to clustering of recruits around their parents 17,18,[35][36][37][38][39] whereas the stable mode is related to clustering in locations that are independent from the parent locations, due, for example, to animal seed dispersal 41 or canopy gaps 40 . This result calls for a closer examination of the spatial relationship between the recruits and adults. Indeed, independent placement of recruits from conspecific large trees may not be unusual. For example, Getzin et al. 42 found in detailed analyses of the BCI forest that recruits were for most species spatially independent of large conspecific trees. For the stable mode we could identify conditions for coexistence, and forthcoming work may extend to quantifying the ability of additional mechanisms such as niche differences 7,8 , habitat associations 46 , spatial and temporal relative nonlinearity 7,8 and storage effects 7,8 to alleviate the destabilizing increase of clustering if species become rare.
This study explicitly incorporates spatial patterns in theoretical models of plant communities and combines analytical theory with spatial simulations and field data analysis. Our finding that species with similar attributes may show stable coexistence has profound implications for ecological theory. Furthermore, the multiscale framework we propose here opens exciting new avenues to explain species coexistence through a spatial lens.

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 4 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 18 of individuals of species i on an individual of the focal species f. We used for this purpose phylogenetic distances 49 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 26 . 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 26 .
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 2 and k fh = K fh (R) / π R 2 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) Finally, the factors γ ff = ln(1 + b ff β ff ) (b ff β ff ) −1 and γ fβ = ln(1+ b fβ β ff ) (b fβ β ff ) −1 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: −1 and the total number of individuals being 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) with mK = 1 S ∑ i Ki 1−αih/αii and mα = 1 S ∑ i αih/αii 1−αih/αii being averages over the S species of the community.
All species have positive abundances at equilibrium if the two conditions µ f = α fh J * /α ff K 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 μ. In this case J * /μ is also approximately constant and we can replace K i in equation (14) by (J * /μ)(α ih /αii) −1 and obtain where S is the number of species in the community, and thereforē 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 μ 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 α ) −1 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 11 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 = ∑ S−1 i=1 N * i and αmmM(t) ≪ α fh J * m (that is, the invading species m is at low abundance and does not show strong clustering) we find > α mh J * m , 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 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 34 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≠f n 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 2 / A, and J(t) = ∑ i N 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 λ f = 1.
The average individual fitness λ 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 λ f (t) = N f (t + Δt) /N f (t). A key ingredient of scale transition theory 34 is that the fitness-density covariance is given by cov =λ f (t) −λ f . With equation (3) and γ ff ≈ γ fβ ≈ 1 and n fβ =n fh we find where the mean of the crowding indices is given by n ff (N f ) = ck ff (N f )N f and n fh (N f ) = ck fh (J * − 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, sincen ff +n fh = c(k ff − k fh )N f + ck fh J * 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 n ff +n fh , n ff , n fh and λ f −λ 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 ) −1 ) 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 17 , 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 (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 28 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 σ 2 ). 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 42 in a simple way, such as contagious seed dispersal by animals 50 or forest gaps that may imprint clumped distributions of recruits of pioneer species 40 . 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 −1 ) and an overall 5 yr mortality rate (10%) similar to those of trees with dbh ≥ 10 cm in the BCI plot 51 .
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.

Data availability
The data that support the findings in this manuscript (and the raw data for Figs

Code availability
The source code of the simulation model is provided in the Supplementary Information. Fig. 1 | examples for intraspecific variability in the crowding indices. a, b, c, Malus baccata, Baihua plot. d, e, ), middle: number of heterospecific neighbours (n kfh ), and right: heterospecifics neighbours weighted by their relative competitive effect β fi /β ff (n kfβ ). We used a 10 m plant neighbourhood 19,26 and phylogenetic similarity as surrogate for pairwise interaction strength 19,26 . Solid blue lines show gamma distributions with the same mean and variance-to-mean as the observed distributions, and the vertical red lines indicates the mean values. The 95% percentiles for the error indices quantifying the departures from a Gamma distribution for all 289 focal species were 0.051, 0.027, and 0.022 for n kff , n kfh , n kfβ , respectively, indicating a reasonable fit.

NAturE EcoLogy & EvoLutioN
Extended Data Fig. 2 | Characteristics of the neighbourhood crowding indices. Distribution of the mean and the variance-to-mean ratio of the crowding indices for the different species at each forest plot with boxplots indicating 10th, 25th, 50th, 75th, 90th percentiles and outliers. a. The mean n ff of the conspecific neighbourhood crowding index n kff over species. b. The mean n fh of the heterospecific neighbourhood crowding index n kfh over species. c. The mean n fβ of the interaction neighbourhood crowding index n kfβ over species. d. The variance-to-mean ratio b kf of the conspecific neighbourhood crowding index n kff over species. e. The variance-to-mean ratio b kh of the heterospecific neighbourhood crowding index n kfh over species. f. The variance-to-mean ratio b kβ of the interaction neighbourhood crowding index n kfβ over species. The neighbourhood radius was R = 10 m. We used for the analysis all individuals with dbh ≥ 10 cm and included focal species with more than 50 individuals. For plot names see Supplementary Table 1. Fig. 3 | Correlation between different crowding indices. We estimated for all individuals of a species f the correlation between their crowding indices n kff and n kfβ (a) and n kfh and n kfβ (b). The crowding index n kff counts the conspecific neighbours of individual k within distance R = 10 m, n kfh counts the corresponding number of heterospecific neighbours, and n kfβ weights each heterospecific neighbour by its relative competition strength β fi /β ff . The boxplots show the distribution of the Pearson correlation coefficients for each focal species, separately for the nine forest plots, indicating 10th, 25th, 50th, 75th, 90th percentiles and outliers. We used for the analysis all individuals with dbh ≥ 10 cm and included focal species with more than 50 individuals. For plot names see Supplementary Table 1.

nature research | reporting summary
April 2020 Corresponding author(s): Thorsten Wiegand, Xugao Wang Last updated by author(s): Feb 26, 2021 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted  Table. This table contains also the raw data (including sample sizes) and estimation of the p-values shown in Extended Data Figure 8.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability The data that support the findings in this manuscript (and the raw data for figures 2, 3 and 4 and Extended Data Figures 2, 3, 4 and 8) can be found in the Supplementary Information as Data Table. To generate this data, we used the raw census data of the ForestGEO network that can only be shared on request because most PI's did not make them publicly available. For data request see https://forestgeo.si.edu/sites-all .

nature research | reporting summary
April 2020 Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. All studies must disclose on these points even when the disclosure is negative.

Study description
This study develops a theoretical multiscale framework to reveal how pattern-forming processes operating at the level of individual trees translate into mesoscale spatial patterns and how they influence macroscale population and community dynamics. The approach includes analysis of mathematical and simulation models and comparison of the simulation results with data of fullymapped forest plots of the Forest Global Earth Observatory (ForestGEO).

Research sample
We used ForestGeo ( https://forestgeo.si.edu/explore-data ) data sets of nine large forest dynamics plots of areas between 20 and 50 ha.

Sampling strategy
The nine forest plots included in the study have been completely censused for trees at least 10 cm in diameter, so there is no sampling within each forest.

Data collection
The study did not involve data collection. The data from model simulations was collected as described in the methods.
Timing and spatial scale The study did not involve data collection. We used for each forest dynamics plot (size ranging between 20 and 50ha) data of one census.

Data exclusions
We used only trees with sizes >= 10cm dbh (diameter at breast height), which is an established approach in the field. This size threshold excludes most of the saplings and enables comparisons with previous spatial analyses. Tree species with a minimum of 50 individuals (dbh >= 10 cm) were used in the study to ensure that estimates of spatial patterns were reliable.

Reproducibility
The estimation of the measures of spatial patterns and the simulation model is specified exactly in the publication such that it could be reproduced. Additionally the simulation code (that includes estimation of measures of spatial patterns) is added as Supplementary File.