Articles
https://doi.org/10.1038/s41559-021-01440-0
Consequences of spatial patterns for coexistence 
in species-rich plant communities
Thorsten Wiegand   1,2 ✉, Xugao Wang   3 ✉, Kristina J. Anderson-Teixeira   4,5, Norman A. Bourg4, 
Min Cao6, Xiuqin Ci7,8, Stuart J. Davies5, Zhanqing Hao3,9, Robert W. Howe10, W. John Kress11, 
Juyu Lian12, Jie Li7,8, Luxiang Lin6, Yiching Lin13, Keping Ma   14, William McShea4, Xiangcheng Mi   14, 
Sheng-Hsin Su   15, I-Fang Sun16, Amy Wolf10, Wanhui Ye12 and Andreas Huth1,2,17
Ecology cannot yet fully explain why so many tree species coexist in natural communities such as tropical forests. A major diffi-
culty is linking individual-level processes to community dynamics. We propose a combination of tree spatial data, spatial statis-
tics 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.
Understanding the mechanisms that maintain high species However, spatial patterns and population-level interaction coeffi-diversity in plant communities such as tropical forests has cients emerge at the mesoscale from the microscale behaviour of long challenged ecologists1 and has stimulated major efforts individuals and their interactions with other individuals and the 
in field and theoretical ecology2–5. However, despite a multitude environment. Therefore, studying the impact of spatial patterns on 
of coexistence mechanisms that have been proposed6 and recent species coexistence requires multiscale approaches such as spatial 
advances in coexistence theories7–18, this fundamental question has moment equations18,23,24 that incorporate pattern-forming processes 
not been fully resolved8,9,14. For example, theoretical models indicate operating at the level of individuals and translate these into popula-
that stable coexistence is difficult to reach in large communities11,13. tion and community dynamics.
We argue that consideration of spatial patterns of plant individuals, We propose here such a multiscale approach. To this end, 
such as intraspecific clustering and interspecific segregation, may we first derive population-level interaction coefficients αfi from 
allow for a better understanding of mechanisms of coexistence in individual-level interaction coefficients βfi and neighbourhood 
species-rich communities17. crowding indices19,21 that are commonly used to describe interac-
Although many studies suggest that spatial patterns and neigh- tions among tree individuals at the microscale, and then incorporate 
bourhood effects may play an important role in diversity mainte- the emerging coefficients αfi into analytical macroscale models. Our 
nance17–21, the integration of spatial patterns into coexistence theories approach is based on separation of timescales (adiabatic approxima-
of species-rich communities is difficult. A major difficulty is link- tion25), given that mesoscale spatial patterns usually build up quickly 
ing spatial processes at the individual level to community dynam- and approach a quasi-steady state whereas the macroscale state vari-
ics. One reason for this is a scale mismatch. The analytical models ables (for example, abundances) change slowly23. Therefore, we do 
that form the basis of most coexistence theories7,8,11–16,22 have state not need to describe the dynamics of the quick mesoscale patterns 
variables that operate at the macroscale (that is, the population or explicitly (as, for example, is done in approaches based on moment 
community-level abundances), use parameters that describe average equations18,23,24) but concentrate instead on spatial patterns that 
‘mesoscale’ properties of the individuals (such as population-level transport the critical information from the microscale into mac-
interaction coefficients and demographic rates) and often rely on roscale models. This approach requires information on mesoscale 
‘mean-field’ approximations18,23 where spatial patterns are neglected. spatial patterns that can be obtained from fully mapped forest plots 
1Department of Ecological Modelling, Helmholtz Centre for Environmental Research – UFZ, Leipzig, Germany. 2German Centre for Integrative Biodiversity 
Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany. 3CAS Key Laboratory of Forest Ecology and Management, Institute of Applied Ecology, Chinese 
Academy of Sciences. 4Conservation Ecology Center, Smithsonian Conservation Biology Institute, Front Royal, VA, USA. 5Forest Global Earth Observatory 
(ForestGEO), Smithsonian Tropical Research Institute, Washington, DC, USA. 6CAS Key Laboratory of Tropical Forest Ecology, Xishuangbanna Tropical 
Botanical Garden, Chinese Academy of Sciences. 7Center of Conservation Biology, Core Botanical Gardens, Chinese Academy of Sciences. 8Centre for 
Integrative Conservation, Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences. 9School of Ecology and Environment, Northwestern 
Polytechnical University. 10Department of Natural and Applied Sciences, University of Wisconsin–Green Bay, Green Bay, WI, USA. 11Department of Botany, 
National Museum of Natural History, Smithsonian Institution, Washington, DC, USA. 12Key Laboratory of Vegetation Restoration and Management 
of Degraded Ecosystems, South China Botanical Garden, Chinese Academy of Sciences. 13Department of Life Science, Tunghai University. 14State Key 
Laboratory of Vegetation and Environmental Change, Institute of Botany, Chinese Academy of Sciences. 15Taiwan Forestry Research Institute. 16Center 
for Interdisciplinary Research on Ecology and Sustainability, National Dong Hwa University. 17Institute of Environmental Systems Research, University of 
Osnabrück, Osnabrück, Germany. ✉e-mail: thorsten.wiegand@ufz.de; wangxg@iae.ac.cn
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles NAturE EcoLogy & EvoLutioN
such as those of the Forest Global Earth Observatory (ForestGEO) Data Fig. 3b) with a common factor Bf (that is, nkfβ ≈ Bfnkfh). This 
network4. result suggests operation of diffuse neighbourhood competition, in 
More specifically, we (1) derive species-level interaction coeffi- which the competition strength of heterospecifics is on average a 
cients from individual-level interactions using empirical informa- factor Bf lower than that of conspecifics.
tion on spatial patterns in nine ForestGEO megaplots4, (2) integrate The integral of equation (2) can be solved analytically for inde-
the resulting species-level interaction coefficients into analytical pendent gamma distributions px and py and yields
macroscale multispecies models and (3) study their consequences 
for multispecies dynamics and coexistence. ( ( ))
survf = sf exp −βff γffn̄ff + γfβn̄fβ (3)
Results and discussion
Species-level interaction coefficients. We first derive species-level where n̄ff  and n̄fβ are the average values of the crowding indices nkff 
interaction coefficients from individual-level neighbourhood and nkfβ, respectively, and γff and γfβ contain the variance-to-mean 
crowding indices19,21,26 that quantify how the performance of a focal ratios of the gamma distributions px and py, respectively, but in our 
individual depends on interactions with its neighbours. To this end, case have values close to one (Methods).
we describe the survival rate of a focal individual k of species f in The last step in deriving pairwise population-level interaction 
dependence on the local number of neighbours as coefficients is to relate the averages of the different crowding indices 
to the macroscale population abundances Nf(t). We accomplish this 
   by taking advantage of connections between crowding indices and 
   the summary functions of spatial point process theory21. The mean 
�
s s exp kf = f −βff nkff + (βfi/βff)n

kfi , (1a) of the crowding index nkff (that is, the mean number of further con-  i̸=f
   specific neighbours within distance R) is proportional to Ripley’s K, 
� �� �
n a well-known quantity in point process theory28,29kfβ :
where sf is a density-independent background survival rate of spe- n̄ff = Kff (R)Nf(t)/A (4a)
cies f; the crowding indices nkff, nkfi and nkfh are the number of con-
specifics, neighbours of species i and heterospecific neighbours where Kff(R) is the univariate K function for species f and A is the 
within distance R of a focal individual k, respectively; the subscript area of the observation window. The K function describes the spatial 
‘h’ indicates all heterospecifics together (that is, nkfh = ∑i≠f nkfi) and pattern of conspecifics within a neighbourhood distance R, indicat-
the crowding index nkfβ weights each heterospecific neighbour by ing clustering if Kff(R) > πR2, a random pattern if K (R) = πR2ff  and 
its relative competition strength βfi/βff (equation (1a)), with βfi being regularity if Kff(R) < πR2. In the following, we use the normalized 
the individual-level interaction coefficients between species f and K function kff (R) = Kff (R) / R2π  to quantify the spatial neighbour-2
i (Fig. 1a–c). The corresponding population-level survival rate is hood patterns, and therefore n̄ = k (R) πRff ff A Nf(t).
given by Analogously, the mean number of heterospecific neighbours is 
given by
  
∑
survf = sf exp αff Nf(t) + (afi/αff)N 2i(t) ∑− , (1b) n̄fh = k
πR
fh (R)
i̸=f A Ni(t), (4b)
i≠ f
where Ni(t) is the abundance of species i at time step t and αfi is the where the bivariate normalized K function kfh(R) indicates segrega-
population-level interaction coefficient between species f and i. To tion to heterospecifics (subscript ‘h’) within distance R if kfh(R) < 1. 
estimate survf we average the survival rates skf (equation (1a)) of all Independent placement occurs if kfh(R) = 1, and attraction if 
individuals k of species f: kfh(R) > 1.
Motivated by the finding nkfß ≈ Bf nkfh (Extended Data Fig. 3b) we N (t)
∫ ∫
1 f ( )∑ rewrite the mean crowding index n̄fβ assurvf = Nf (t
skf = sf exp −βff (x+ y) pxpydxdy (2))
k=1 x y n̄fβ = Bf n̄fh, (5)
︸︷︷︸
n̄fβ
where px and py are the distributions of the crowding indices x = n n̄kff fh
and y = nkfβ for individuals of species f, respectively.
To determine the distributions px and py, we analysed forest inven- where the point process summary function Bf indicates how much 
tory data from nine 20–50 ha forest dynamics plots (Supplementary the competition strength of one heterospecific neighbour differs on 
Table 1) in the ForestGEO network4. We used phylogenetic similar- average from that of one conspecific neighbour. The values of Bf 
ity between tree species as a surrogate for the relative competition depend mainly on the individual-level interaction coefficients βfi 
strength βfi/βff because it is available for the species in the nine plots but also on the spatial pattern of the different species and their rela-
(Methods). This is an established approach in species-rich commu- tive abundances (equation (12)).
nities19,26,27 to approximate niche differences in the absence of other Inserting the expressions for n̄ff  and n̄fβ into equation (3) and 
data. comparing with equation (1b) leads to our first main result, the ana-
The number of con- and heterospecific neighbours and the het- lytical expressions of the population-level interaction coefficients:
erospecific interaction index nkfβ vary widely among conspecifics 
and can be described by gamma distributions (Fig. 1 and Extended αff = c γff kffβff intraspecific interactions of species f (6a)
Data Figs. 1 and 2). Detailed analysis of the empirical crowding 
indices reveals additional relationships that are relevant for our sub-
sequent analysis. First, we find that the crowding indices n α = cγ k β B interspecific interactionswith species i (6b)kff and nkfβ fi fβ fh ff f
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 cor- with scaling constant c = πR2/A. Notably, equation (6b) indicates 
relations between the two crowding indices nkfh and nkfβ (Extended that the emerging population-level interaction coefficients αfi are 
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
a b c Heterospecific neighbours
Conspecific neighbours Heterospecific neighbours weighted by βfi /βfi
k
d Conspecifics e Heterospecifics f Interactions
300 200 160
140
250
150 120
200 100
150 100 80
60
100
50 40
50 20
0 0 0
0 5 10 15 20 0 10 20 30 40 50 60 70 0 10 20 30 40
nkff nkfh nkfβ
Fig. 1 | Neighbourhood crowding indices describe individual-level interactions and their intraspecific variability. a, The conspecific crowding index nkff is 
the number of conspecific neighbours (filled red circles) within distance R (black circle) of the focal individual k (filled red square). b, The heterospecific 
crowding index nkfh is the number of heterospecific neighbours (filled grey circles) within distance R of the focal individual k. c, The heterospecific 
interaction crowding index nkfβ additionally weights heterospecifics by their relative competitive effect βfi/βff, symbolized by the arrows. Different colours 
indicate different species. d, Distribution of the number nkff of conspecific neighbours with diameter at breast height (dbh) ≥ 10 cm of the species 
Castanopsis cuspidata of the 25 ha Fushan plot. e, Corresponding distribution of heterospecific neighbours nkfh. f, Corresponding distribution of the crowding 
index nkfβ. In d–f, blue lines show gamma distributions with the same means and variance-to-mean ratios as the observed distributions, and the vertical 
black line indicates the mean value. See Supplementary Data Table 1 for additional examples.
the same for all heterospecifics (that is, αfi = αfh for i ≠ f). Thus, even neighbours. Spatial patterns therefore have a strong potential to 
if the individual-level interactions coefficients βfi differ among spe- alter the outcome of deterministic individual-level interactions.
cies pairs, the emerging population-level interaction coefficients 
αfi have a substantially simpler structure. This phenomenon is an Conditions for coexistence in the multiscale model. To study 
example of simplicity emerging from complex species interactions30 the consequences of the emerging spatial patterns for community 
and is likely to occur only in species-rich communities9,12,31. dynamics and coexistence we insert the population-level interaction 
The population-level interaction coefficients αfi depend on coefficients αfi (equation (6)) into a simple macroscale model
several factors that can influence the macroscale balance between 
intra- and interspecific competition: (1) intraspecific clustering kff Nf(t+Δt)−Nf(t) = N (t)
and interspecific segregation kfh (Fig. 2a,b), (2) the relative com- Δt
f
[ ( )]
petition strength Bf of one heterospecific neighbour (Fig. 2c) and (7)( ) ∑
(3) the shape of the response of survival to crowding (γ  and γ , rf − 1 + sf exp −αffNf (t)− αff fβ fh Ni(t)
which contain the variance-to-mean ratios of the distribution of i≠ f
the crowding indices). Note that absence of spatial patterns (that is, In this model we assume that survival is governed by neighbour-
kff = 1 and kfh = 1) and a linear approximation of equation (3) lead hood competition with αfi = αfh, and the number of recruits of spe-
to γff = γfβ = 1 and direct proportionality αfi = cβfi of the individual- cies f during a time step Δt is given by rfNf(t), where rf is the per 
and population-level interaction coefficients, as assumed by Lotka– capita reproduction rate of species f.
Volterra models. The carrying capacity of species f (that is, the equilib-
The information on spatial patterns extracted from the inven- rium of equation (7) with Ni(t) = 0 for all i ≠ f) is given by 
tory data of our nine forests allows us to estimate the relative ( )1−rf −1
population-level interaction coefficients α /α  for all pairs of spe- Kf = − ln s αff . Note that our theory also applies, after fi ff f
cies i and f (Fig. 3 and Extended Data Fig. 4). For example, at redefinition of the carrying capacity, to alternative macroscale 
BCI, the values of αfi/αff differ substantially from the correspond- models (Supplementary Table 3). From equation (7) we find 
( )
ing individual-level coefficients βfi/βff (Fig. 3a,c), and for 83% of K N∗ α= + fh
∑
N∗ = N∗ 1 αfh α− + fh J∗, where N∗f  is the 
all species pairs, we find αfi/αff < βfi/βff. Thus, the mesoscale spatial 
f f αff i f α αi̸=f ff ff
patterns can reduce, at the population level, the strength of hetero- abundance of species f in equilibrium and J* the equilibrium com-
∑
specific interactions relative to conspecific interactions by ‘diluting’ munity size (that is, J∗ = N∗i i ; see also equation (14)). This leads, 
encounters with heterospecific neighbours relative to conspecific under the assumption that the population-level interaction coeffi-
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Frequency
Articles NAturE EcoLogy & EvoLutioN
a b c
1,000 2.0
1.2
100 1.0
1.5 0.8
10 0.6
1.0 0.4
1
0.2
0.1 0.5 0
AB
a
BS ihu CB
I S S IFS N B S a
I
H T B S S
S N I B S a I S I
C D G S B
C A B ihu CB H T F B BCC G AD S CB ih
u CB DH
S TS F BN CW Ba S X W Ba S X W Ba S
G XS
B
Forest plot Forest plot Forest plot
Fig. 2 | emergent spatial patterns in the nine forest dynamics plots. a–c, Distribution of the values of the different measures of spatial patterns, taken over 
the focal species of the different forest plots. Boxplots show the 10th, 25th, 50th, 75th and 90th percentiles; outliers are indicated by filled black points. 
Intraspecific clustering is indicated by kff > 1 and interspecific segregation by kfh < 1, and the less a heterospecific neighbour competes on average relative to 
a conspecific neighbour, the more Bf 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.
cients αfh are constant (Supplementary text), to a single equilibrium 282 of our 289 focal species, we found µf < 1 and αfh/αff < 1, which 
of the macroscale model for species f means that the two conditions for stable coexistence are indeed sat-
isfied for nearly all species. However, this is not a given, as shown 
( )
1 r by the seven species from BCI with µf > 1 and αfh/αff > 1. Thus, our −
− ln fsf theory is compatible with the observed coexistence of most species 
︷︸︸︷ at our nine forest plots if the assumption of approximately constant 
∗
N∗
αffKf −αfhJ
= (8) population-level interaction coefficients holds.f αff − αfh 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 
that is positive if denominator and numerator are both positive agreement with the predictions of our theory, we find that µf tends to 
or both negative. However, the invasion criterion (equation (18)) be smaller if αfh/αff is smaller (Fig. 4c). Furthermore, the values of µf 
is only fulfilled if both are positive. In this case, equation (8) sug- were, for most species, larger than the expectation of µf for the cor-
gests two different ways a species can go extinct. First, the denomi- responding communities without interspecific variability in µf (equa-
nator indicates that a species with strong clustering kff will show a tion (10a)) but with the same number of species and the same mean 
small equilibrium abundance since in this case αff ≫ αfh (equation values of αfh/αff (Fig. 4c), but all of the values were relatively close to 
(6)). Large values of kff can be expected for species of low abun- the critical value of 1 (the median of all 289 species was 0.938; Fig. 4a).
dance under dispersal limitation, where recruitment happens close 
to conspecific adults. Consequences of spatial patterns for coexistence. Our theory 
Second, the numerator of equation (8) indicates a positive abun- predicts that coexistence requires, in the limit of high species rich-
dance of species f if αffKf > α *fhJ  and αfh/αff < 1. Therefore, we intro- ness, that species approach functional identity with respect to the 
duce a new feasibility index feasibility index µ  (Methods). This resembles neutral theory2,32f , but 
our theory allows for trade-offs among demographic parameters 
αfhJ∗μ = (9) and emerging spatial patterns to reach this equivalence (equation f αffKf (9)). To study the consequences of spatial patterns for coexistence, 
we analysed a symmetric33 version of our model where all species 
that indicates a positive abundance if µf < 1 given that heterospe- have the same parameters and follow the same stochastic rules and 
cific interactions at the population level are weaker than conspecific where all individuals compete identically (leading to Bf = 1). Thus, 
interactions (that is, αfh/αff < 1). The invasion criterion7,8 that tests we eliminate any potential coexistence mechanism other than that 
whether a species with low abundance can invade the equilibrium resulting from spatial patterns.
community of all other species turned out to be basically the same If the mesoscale patterns kff and kfh converge to a stochastic equi-
as the feasibility condition (equation (9)) if the invading species librium, we find that the feasibility and invasion criteria are always 
does not show strong clustering (Methods and equation (18)). Note fulfilled if αfh/αff < 1 since
that we did not assume Allee effects8.
Further analysis that considers the dependency of J* on the val- S αfhα
ues of Kf and αfh/αff shows that the values of µf must be similar for μ
ff
f = αfh < 1 (10a)
all species f to fulfil the condition µf < 1, and that µf can show larger 1+ α (S− 1)ff
interspecific variability if the species richness S is smaller and/or if 
( )−1
the mean of αfh 1 αfhα − α  is smaller (equations (14) and (15)). αff ff (S− 1) fhα
The feasibility index µf therefore governs species assembly by deter- μ
i ff
f = αfh < 1 (10b)
mining the subset of species of a larger species pool that can per- 1+ α (S− 2)ff
sist13, but any addition of a new species changes µf and may lead to 
reassembly of the community. Equations (10a) and (10b) follow from equations (16) and (18), 
Using the observed abundances in the forest plots (and assum- respectively, if αfh/αff is the same for all species f. Thus, the spatial 
ing equilibrium) allows us to test our theory. We can estimate from patterns that emerge at the mesoscale from the individual-level 
the observed abundances the carrying capacities Kf and therefore interactions can stabilize if αfh/αff < 1. The underlying mechanism is 
also the indices µf for all focal species of our nine plots (Fig. 4). For a positive fitness–density covariance34 (Methods).
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Intraspecific pattern kff
Interspecific pattern kkh
Bf
NAturE EcoLogy & EvoLutioN Articles
a b c
2.0
0.6
1.5 αfi /αff
0.4
1.0 βfi /βff
0.5 0.2
0 0
0 0.5 1.0 1.5 2.0
AB BS
a I S S S N I 0 0.5 1.0 1.5
Individual-level β /β W C aih
u
SC
B T F C
DH G B B
fi ff B X
S
Relative interaction coefficients
Forest plot
Fig. 3 | The emergent population-level interaction coefficients. a, Relationship between the relative population-level interaction coefficients αfi/αff and 
the corresponding relative individual-level interaction coefficients βfi/βff, (equation (6)) for the 75 focal species of the BCI plot. The βfi/βff were based 
on phylogenetic dissimilarity. b, Distribution of the values of αfi/αff, taken over the 289 focal species of the different forest plot. Boxplots show the 10th, 
25th, 50th, 75th and 90th percentiles; outliers are indicated by filled black points. For full, separate distributions for tropical, subtropical and temperate 
forests, see Extended Data Fig. 4. c, Example for the distribution of the relative individual- and population-level interaction coefficients for the BCI plot. 
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.
a b c
0.20
No feasibility Unstable
1.0
0.15
0.8
0.10
0.6
0.05 0.4 Stable No
coexistence feasibility
0 0.2
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 B S ua IA S S S
I 0 0.5 1.0 1.5
CBW aih C
B H GT
F
D SB
N
S B
C
Feasibility index µ B X Stabilization α /αf fh ff
Forest plot
Fig. 4 | The feasibility index µf for the 289 focal species of different forest plots. a, The distribution of the index µf (equation (9)) for the 289 analysed 
tree species. The vertical dashed line indicates the median of the distribution. b, Distribution of µf, taken over the 289 analysed species of the different 
forest plot. Boxplots show the 10th, 25th, 50th, 75th and 90th percentiles; outliers are represented by black points. c, Relationship between the feasibility 
index µf and the ratio αfh/αff of heterospecific to conspecific population-level interaction coefficients. Values of µf < 1 indicate a positive abundance. The 
red lines show the dependence of µf on αfh/αff expected for communities without interspecific variability in µf (equation (10a)) with 18 focal species (CBS 
plot, lower line) and 75 focal species (BCI plot, upper line). The neighbourhood radius used was R = 10 m. For the analysis, we used all individuals with 
dbh ≥ 10 cm and included tree species with more than 50 individuals. For plot names, see Supplementary Data Table 1.
To reveal the conditions that can lead to coexistence in a spatially the instability was high clustering of rare species24,35 (Extended Data 
explicit context, we use a spatially explicit and individual-based35 Fig. 6) that completely negated the potentially positive effects of 
implementation of the symmetric version of the multiscale model αfh/αff < 1.
(equation (7) and Methods). Models of this type are able to produce In contrast, community dynamics can be stabilized if recruits 
realistic spatial patterns consistent with mapped species distribu- are placed in small clusters but independent of the location of con-
tions of large forest plots17,35. While our analytical approach in equa- specific adults (Extended Data Fig. 5c). With this mechanism we 
tion (7) only allows us to make simplified assumptions about the mimic canopy gaps40, animal seed dispersal41 or other mechanisms 
spatial component of the recruitment process, the simulation model that can generate clustering independent of parent locations, as 
allows us to explore the role of the spatial component of recruitment found at BCI42. Decoupling clustering from the parent locations 
in more detail. does not lead to the negative relationship between clustering and 
Indeed, the way recruits were placed was critical for coexistence. abundance, and all measures of spatial patterns converged quickly 
Randomly placed recruits produced unstable dynamics (Extended into quasi-equilibrium (Extended Data Fig. 5f,i,o).
Data Fig. 5a) characterized by regularity (mean of kff = 0.92) and The simulation data reveal the spatial coexistence mechanism 
segregation (mean of k  = 0.92), both caused by competition18, and underlying the positive fitness–density covariance34fh  (Extended Data 
the instability was caused by con- and heterospecifics competing Figs. 7 and 8). We find that the emerging spatial patterns lead to a 
equally at the population level (that is, αfh/αff ≈ 1) (Extended Data situation where individuals of a common species are more likely to 
Fig. 5d,g). When we followed the common approach of placing be near more neighbours and tend therefore to experience stronger 
recruits with a kernel around conspecific adults17,18,35–39 to mimic competition (Extended Data Fig. 7). While the number of hetero-
dispersal limitation, we again found unstable dynamics (Extended specific neighbours remains approximately constant, the number of 
Data Fig. 5b), despite intraspecific clustering and interspecific seg- conspecific neighbours decreases with decreasing abundance if clus-
regation (that is, αfh/αff < 1; Extended Data Fig. 5n). The reason for tering does not change with abundance (equation (4a)). However, 
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Proportion Population-level αfi /αff
Feasibility index µf
Proportion
Articles NAturE EcoLogy & EvoLutioN
if clustering increases with decreasing abundance, the rare-species independent from the parent locations, due, for example, to ani-
advantage is weakened and the dynamics become unstable24. mal seed dispersal41 or canopy gaps40. This result calls for a closer 
The data of several ForestGEO forest plots were compatible with examination of the spatial relationship between the recruits and 
a positive fitness–density covariance (Extended Data Fig. 8f–i) as adults. Indeed, independent placement of recruits from conspecific 
they show that, when a species becomes rare, areas of higher con- large trees may not be unusual. For example, Getzin et al.42 found in 
specific crowding tend to have fewer total competitors. Comparison detailed analyses of the BCI forest that recruits were for most spe-
of the results of the stable versus unstable simulation showed that cies spatially independent of large conspecific trees. For the stable 
even relatively weak tendencies in this relationship are sufficient to mode we could identify conditions for coexistence, and forthcom-
stabilize the dynamics (Extended Data Figs. 8a,b). However, this was ing work may extend to quantifying the ability of additional mecha-
not the case for three temperate forest plots where the power-law nisms such as niche differences7,8, habitat associations46, spatial and 
clustering–abundance relationship showed exponents of b < −0.5 temporal relative nonlinearity7,8 and storage effects7,8 to alleviate the 
(that is, species tended to have high clustering at low abundances), destabilizing increase of clustering if species become rare.
but the other plots showed b > −0.5 (Extended Data Fig. 8n–p). This study explicitly incorporates spatial patterns in theoretical 
The apparent contradiction with previous theoretical stud- models of plant communities and combines analytical theory with 
ies18,24,35,37 where intraspecific clustering and interspecific segre- spatial simulations and field data analysis. Our finding that species 
gation could not stabilize community dynamics thus arises as a with similar attributes may show stable coexistence has profound 
consequence of the assumption of placing recruits close to their par- implications for ecological theory. Furthermore, the multiscale 
ents. This finding has important consequences for ecological theory framework we propose here opens exciting new avenues to explain 
because it shows, in contrast to the prevalent view36,43,44, that spatial species coexistence through a spatial lens.
patterns alone can lead to coexistence of multiple species. This is 
even more important since the specific spatial patterns required for Methods
this coexistence mechanism also exist in real forests. 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 
4
Conclusions 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 
Understanding the mechanisms that maintain high species diver- from 36 to 468. All free-standing individuals with diameter at breast height (dbh) 
sity in communities such as tropical forests is at the core of eco- ≥1 cm were mapped, size measured and identified. We focused our analysis here on 
logical theory, but these mechanisms are not yet fully resolved. individuals with dbh ≥ 10 cm (resulting in a sample size of 131,582 individuals) and 
Here, we argue that spatial patterns may play an important role in focal species with more than 50 individuals (resulting in 289 species). The 10 cm 
species coexistence of high diversity plant communities17. To test size threshold excludes most of the saplings and enables comparisons with previous spatial analyses20,35,47,48. Shrub species were also excluded.
this hypothesis, we introduced a multiscale framework that reveals Some of our analyses require estimation of the ratio βfi/βff that describes 
how pattern-forming processes operating at the level of individu- the relative individual-level competitive effect18 of individuals of species i on an 
als translate into mesoscale spatial patterns and how those patterns individual of the focal species f. We used for this purpose phylogenetic distances49 
influence macroscale community dynamics. based on molecular data, given in Myr, that assume that functional traits are phylogenetically conserved19,26,27We showed that the population-level interaction coefficients α  . In this case, close relatives are predicted to fi compete more strongly or to share more pests than distant relatives26. To obtain 
can have, for a broad range of common circumstances, a simpler consistent measures among forest plots, phylogenetic similarities were scaled 
structure than the underlying individual-level interaction coef- between 0 and 1, with conspecifics set to 1, and a similarity of 0 was assumed for a 
ficients β . This simplicity, which emerged from spatially explicit phylogenetic distance of 1,200 Myr, which was somewhat larger than the maximal fi
species interactions30, allowed for an analytical treatment of equi- observed distance (1,059 Myr). This was necessary to avoid discounting crowding effects from the most distantly related neighbours26.
librium, feasibility and invasion conditions of the corresponding 
macroscale models (equations (8–10)). Inserting the emerging αfi Observed spatial patterns at species-rich forests. Figure 1 and Supplementary 
coefficients into macroscale community models (for example, equa- Data Table 1 show the intraspecific variation in our three crowding indices nkff, nkfh 
tion (7); Supplementary Table 3) should, in principle, allow us to and nkfβ that can be approximated by gamma distributions. To assess how well the 
take advantage of macroscale theory9,11–13,45. However, our results gamma distribution described the observed distribution, we used an error index defined as the sum of the absolute differences of the two cumulative distributions 
also indicate that the population-level interaction coefficients may divided by the number of bins (spanning the two distributions). The maximal 
not be temporally constant as commonly assumed but depend on value of the error index is one, and a smaller value indicates a better fit.
spatial patterns that may change with abundance. This is especially Equations (6, 8 and 9) relate the measures of the emerging spatial patterns (that 
likely if recruitment is mainly located close to the parents. is, kff, kfh and Bf) to macroscale properties and conditions for species coexistence. 
It is also possible to expand our framework to take into account 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. 
more detailed neighbourhood crowding indices that consider not only We estimate the key quantities of equations (8) and (9) directly from the forest 
the number of neighboured trees of a given species but also their dis- plot data (Fig. 4), with the exception of the carrying capacities Kf, which were 
tance and size19,21,26. This requires redefinition of the quantities k  and indirectly estimated from the observed species abundances (assuming approximate ff
kfh that describe intraspecific clustering and interspecific segregation, equilibrium; equation (8) and Supplementary Data Table 1). This allowed us to 
respectively, but does not change the overall structure of our equa- estimate the feasibility index µf (equation (9)). Because statistical analyses with individual-based neighbourhood models19,26 based on neighbourhood crowding 
tions. A special strength of our approach is that the population-level indices have shown that the performance of trees depends on their neighbours 
interaction coefficients contain measures of spatial neighbourhood for R between 10 and 15 m, we estimate all measures of spatial neighbourhood 
patterns that can be directly estimated from fully mapped forest plots4. patterns with a neighbourhood radius of R = 10 m. Analyses with R = 15 or R = 20 
Together with additional information, this may allow for estimating gave similar results.
network structures as well as stability of the whole community. The spatial multispecies model and equilibrium. We use a general macroscale 
Our analysis revealed that communities of competing species model to describe the dynamics of a community of S species:
can show a stable mode where the mesoscale patterns converge   
quickly into quasi-equilibrium and an unstable mode where nega- Nf (t + Δt) Nf (t) ∑− ( )
= N (t) r − 1 + s exp−α N (t) − α N (t)
tive relationships between species clustering and abundance emerge Δt f f f ff f fi ii≠ f
(Extended Data Figs. 5 and 6). The two modes are governed by (11)
the way species clustering is generated: the well-known unstable 
mode is related to clustering of recruits around their parents17,18,35–39 where rf is the mean number of recruits per adult of species f within time step Δt, 
whereas the stable mode is related to clustering in locations that are sf is a density-independent background survival rate of species f and the αfi are the 
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
population-level interaction coefficients, yielding αff = c γff kff βff and αfi = c γfβ kfh βff Bf Invasion criterion. Using the population-level interaction coefficients (equation 
(equation (6)). The βfi are the assumed individual-level interaction coefficients (6)) in the macroscale model, we now derive conditions for coexistence based on 
between individuals of species i and f; k 2ff = Kff(R) / π R  and kfh = Kfh(R) / π R2 measure the invasion criterion7,8 for a species m. The growth rate of an invading species m 
intraspecific clustering and interspecific segregation, respectively, with Kff(R) with low density M(t) into the equilibrium community of all other S – 1 species 
being the univariate K function for species f and Kfh(R) the bivariate K function Ni(t) should be positive; thus, with equation (7), we have
describing the pattern of all heterospecifics ‘h’ around individuals of species f. A is 
( )
the area of the observation window. S−1∑
Following equation (5), B  can be estimated as (rm − 1) + sm exp −αmmM (t) − αmh N
∗
i > 0. (17)
f i=1
∑ [ ] ∑
n̄f i f ckfiNi (t
β
) fi k N β(t) fi ∑
β β i ∗ S−1 ∗ ∗≠B ff
i=f fi̸ βff
f = = ∑ [ ] = ∑n̄ k N t , (12)
Considering that Jm = i 1 Ni  and αmmM(t) ≪ αfhJm (that is, the invading =
fh ckfiNi (t) fh i f i ( ) species m is at low abundance and does not show strong clustering) we find i f ( )≠ ≠ ln 1−rm > α J∗− s mh m, and by dividing by αmm we obtain the invasion conditionm
and is the weighted average of the relative individual-level interaction coefficients αmhJ∗
βfi/βff between species i and the focal species f, weighted by the mean number of μ
i
m =
m < 1 (18)
individuals of species i in the neighbourhoods of the individuals of the focal species αmmKm
(that is, c kfi Ni(t)). For competitive interactions, Bf ranges between zero and one; which is basically identical to the condition for feasibility (equation (9)), but 
Bf = 1 indicates that heterospecific and conspecific neighbours compete equally, and here the community size J∗ of the reduced community appears instead of the 
smaller values of Bf indicate reduced competition with heterospecific neighbours. 
m
equilibrium community size J* of all species, including species m. Thus, a new 
The denominator can be rewritten in terms of segregation kfh to all heterospecifics species m is more likely to invade if it has a high value of the carrying capacity K  
and the total number of heterospecifics ∑ mi≠f Ni(t). and if it more strongly reduces heterospecific interactions relative to conspecific 
The analytical expression of the equilibrium (equation (8)) relies on the interaction (that is, αmh/αmm is smaller). However, if the species is too efficient (that assumption that the values of Bf are approximately constant in time. This is, has too large a capacity K  and/or too low an α /α ) it may increase the value 
assumption may not apply in our model during the initial burn-in phase of the m mh mmof J* too much (equation (14)), thereby causing the extinction of the weakest species 
simulations if the βfi/βff show large intraspecific variability (Supplementary Text with the highest values of µ  (that is, a too-low value of K  and a too-high value of 
and Figs. 1–5). The underlying mechanism is the central niche effect introduced f mαfh/αff). Equation (18) also suggest that an equilibrium with µ  > 1 and α /α  > 1 by Stump45 where a species has reduced average fitness if it has high niche overlap f mh mmwill be unstable.
with competitors.
Finally, the factors γff = ln(1 + bff β −1ff) (bff βff)  and γfβ = ln(1+ bfβ βff) (bfβ βff)−1 Fitness–density covariance. To place our new spatial coexistence mechanism in 
describe the influence of the variance-to-mean ratios bff and bfβ of the gamma the context of existing coexistence theory, we apply scale transition theory34 to our 
distribution of the crowding indices nkff and nkfβ, respectively. For high survival model version where spatial effects are the only potential coexistence mechanism 
rates during one time step (for example, >85%), the values of γff and γfβ are close to (that is, all species have the same parameters and all individuals compete equally; 
one; in this case the exponential function in equation (1a) can be approximated by β /β  = 1, B  = 1).
its linear expansion and γff = γfβ = 1.
fi ff f
Following equation (1a), the expected fitness of an individual k of a focal 
In equilibrium we have (Nf(t + Δt) ‒ Nf(t))/Δt = 0, which leads, with equation species f (that is, its expected contribution to the population after some defined 
(7), to: interval of time Δt) in the macroscale model (equation (7)) is
( α )( α )−1   
N K fh J 1 fh∗ ∗f = f − − (13) ∑αff αff λkf = rf + sff (Wk) = rf + sf exp−βff nkff + nkfi (19)
( ) i≠ f
( )
with K ln 1−rf α −1f = − s ff  and the total number of individuals being f ( )
N ( )∗∑J∗ = N∗. Rewriting equation (13) yields Kf f 1 α= fh α+ fh . For where Wk = nkff + ∑i≠fnkfi is the fitness factor of individual k, f(Wk) = exp(−βff  Wk) i i J∗ J −∗ αff αff is the fitness function and nkff and ∑i≠f nkfi are the number of conspecific and 
αfh/αff < 1 we therefore find K *f < J , which indicates that a multispecies forest would heterospecific neighbours, respectively, of individual k within distance R. The 
host more individuals than a monoculture. To estimate J* we sum equation (13) spatial average of the fitness factor over the entire plot is
over all species i and find ∑J∗ = K
∑
i ∗ αih/αii
1−αih/α J 1 α /α . Therefore, we obtain− ∑ii −i i ih ii W̄k = cNf (t) + c Ni(t) = cJ(t) (20)
S ( S )−1
∑
J K∗ i
∑
1 αih/αii SmK
i̸=f
= 1 / + 1 / = 1 Sm (14)i 1 − αih αii i 1 − αih αii + α= = where c = πR2 / A, and J(t) = ∑iNi(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 
with ∑ ∑m = 1 KK i 1 αih/αiiS 1−αih/α  and mα = S 1 α /α  being averages over the S species ii ih ii Supplementary Figs. 1 and 2), we find for the spatial average fitness λ̄f = 1.−of the communiity. i The average individual fitness λ̃f(t) of a focal species f is the average of 
All species have positive abundances at equilibrium if the two conditions λk,f over all individuals k of species f and can be estimated for the macroscale 
µf = αfhJ*/αffKf < 1 and αfh/αff < 1 are met (see equation (9)). We now show that the model (equations (6 and 7)) as λ̃f (t) = Nf (t + Δt) /Nf (t). A key ingredient 
chance that these conditions are satisfied for all species is larger if the values of µf of scale transition theory34 is that the fitness–density covariance is given by 
show little intraspecific variability. To understand this, we assume that the µf can be cov = λ̃f (t) − λ̄f . With equation (3) and γff ≈ γμ̄ fβ
 ≈ 1 and n̄fβ = n̄fh we find
approximated by their mean . In this case J∗/μ̄ is also approximately constant and 
we can replace K ∗ −1i in equation (14) by (J /μ̄)(αih/αii)  and obtain λ̃f (t) − λ̄f = (rf − 1) + sf exp(−βff(n̄ff + n̄fh)) (21)
( )
S J∗ αih ( S α )ih −1μ̄ α J∗∑ ii ∑ α Smα where the mean of the crowding indices is given by n̄ (NJ 1 ii ff f
) = ckff(Nf)Nf  and 
∗
= 1 α +
∗
ih 1 α = (15)ih n̄fh(Nf) = ckfh(J − Nf) (equation (4)). Therefore, if clustering kff and segregation i 1 − α i 1 − α μ̄ 1 + Sm= ii = ii α kfh are independent from abundance Nf, more abundant species have more 
neighbours, since
where S is the number of species in the community, and therefore
n̄ff + n̄fh = c(kff − kfh)N ∗f + ckfhJ (22)Sm
μ̄ α= 1 Sm < 1 (16)+ α 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.
Thus, in the case of a perfect interspecific balance in µf we always have a Extended Data Fig. 7 shows the quantities n̄ff + n̄fh, n̄ff , n̄fh and λ̃f − λ̄f  plotted 
feasible equilibrium if αfh/αfh < 1, and species can go extinct only if the intraspecific over abundance Nt for data generated by our spatially explicit simulation model for 
variability in µf becomes too large. The smaller the mean value of µf, the more the scenarios of stable and unstable dynamics (Extended Data Fig. 5b,c). Indeed, 
variability in µf is allowed. Equation (16) shows that μ̄ is smaller if the number S of the stable simulations show a positive fitness–density covariance, however, there is 
species in the community is smaller and/or if the mean value of mα is smaller. no such trend for the dynamics of the unstable community (Extended Data  
Equation (16) also suggests that communities with more species need to show Fig. 7g,h).
stronger species equivalence in µf because the term S m (1 + S m )−1α α  approaches Spatial patterns will act as positive fitness–density covariance if, when a 
a value of one for a large number of species S. This finding mirrors the results of species becomes rare, areas of high conspecific crowding have fewer competitors. 
analyses of Lotka–Volterra models with random interaction matrices11 that showed We tested this for the data generated by our simulation model and for the nine 
that the larger the number of species S, the more difficult it becomes to generate a forest plots (Extended Data Fig. 8). We could estimate for each focal species f the 
feasible community. covariance between the number of conspecific neighbours (that is, nkff) and the 
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles NAturE EcoLogy & EvoLutioN
total number of neighbours (that is, nkff + nkfh) and demand that the covariance network that can only be shared on request because most PIs have not made them 
should be mostly positive and larger for more abundant species. However, since publicly available. For data requests see https://forestgeo.si.edu/sites-all.
the quantity nkff appears in this test on both sides, a positive covariance can be 
expected. To compensate for this artefact, we instead used the covariance between Code availability
the local dominance of conspecifics in the neighbourhood of individuals k (that is, 
−1 The source code of the simulation model is provided in the Supplementary dkff = nkff (nkff + nkfh) ) and total number of neighbours (that is, nkff + nkfh) (Extended Information.
Data Fig. 8).
Spatially explicit simulation model. The model is a spatially explicit and Received: 17 March 2020; Accepted: 1 March 2021;  
stochastic implementation of the spatial multispecies model (equation (7)), similar Published: xx xx xxxx
to that of May et al.35,37 and Detto and Muller-Landau17, and simulates the dynamics 
of a community of S tree species in a given plot of a homogeneous environment References
(for example, 50 ha) in 5 yr time steps adapted to the ForestGEO census interval  1. Hutchinson, G. E. The paradox of plankton. Am. Nat. 95, 137–147 (1961).
(Extended Data Fig. 5 and Supplementary Figs. 1 and 2). Only reproductive (adult)  2. Hubbell, S. P. The Unified Neutral Theory of Biodiversity and Biogeography 
trees are considered, but size differences between them are not considered. During (Princeton Univ. Press, 2001).
a given time step the model first simulates stochastic recruitment of reproductive  3. Usinowicz, J. et al. Temporal coexistence mechanisms contribute  
trees and placement of recruits, and second, stochastic survival of adults that to the latitudinal gradient in forest diversity. Nature 550,  
depends on the neighbourhood crowding indices for conspecifics (nkff) and 105–108 (2017).
heterospecifics (nkfβ) (but excluding recruits). In the next time step, the recruits  4. Anderson-Teixeira, K. J. et al. CTFS-ForestGEO: a worldwide network 
count as reproductive adults and are subject to mortality. No immigration from a monitoring forests in an era of global change. Glob. Change Biol. 21,  
metacommunity is considered. To avoid edge effects, torus geometry is assumed. 528–549 (2015).
The survival probability of an adult k of species f is given by  5. Comita, L. S. et al. Testing predictions of the Janzen–Connell hypothesis: a 
( )
( )
sf exp −βff nkff + nkfβ  (equation (1a)). The two neighbourhood indices nkff meta-analysis of experimental evidence for distance and density-dependent 
and nkfβ describe the competitive neighbourhood of the focal individual k and sum seed and seedling survival. J. Ecol. 102, 845–856 (2014).
up all conspecific and heterospecific neighbours, respectively, within distance R,  6. Wright, J. S. Plant diversity in tropical forests: a review of mechanisms of 
but weight them with the relative individual-level interaction coefficients β /β  species coexistence. Oecologia 130, 1–14 (2002).fi ff
(refs. 19,21,26).  7. Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. 
Each individual produces on average r  recruits, and their locations are Syst. 31, 343–366 (2000).f
determined by a type of Thomas process28 to obtain clustering. To this end, the  8. Barabás, G., D’Andrea, R. & Stump, S. M. Chesson’s coexistence theory.  
spatial position of the recruits is determined by two independent mechanisms. Ecol. Monogr. 88, 277–230 (2018).
First, a proportion 1 – p  of recruits is placed stochastically around randomly  9. Levine, J. M., Bascompte, J., Adler, P. B. & Allesina, S. Beyond pairwise d
selected conspecific adults by using a two-dimensional kernel function (here a mechanisms of species coexistence in complex communities. Nature 546, 
Gaussian with variance σ2). This is the most common way to generate species 56–64 (2017).
clustering in spatially explicit models17,18,35–39. Specifically, we first randomly  10. HilleRisLambers, J., Adler, P. B., Harpole, W. S., Levine, J. M. & Mayfield, M. 
select one parent for each of these recruits among the conspecific adults and M. Rethinking community assembly through the lens of coexistence theory. 
then determine the position of the recruit by sampling from the kernel. Second, Annu. Rev. Ecol. Syst. 43, 227–248 (2012).
the remaining proportion pd of recruits is distributed in the same way around  11. Stone, L. The feasibility and stability of large complex biological networks: a 
randomly placed cluster centres that are located independently of conspecific random matrix approach. Sci. Rep. 8, 8246 (2018).
adults. This mode mimics spatial clustering of recruits independent of the parent  12. Saavedra, S. et al. A structural approach for understanding multispecies 
locations42 in a simple way, such as contagious seed dispersal by animals50 or forest coexistence. Ecol. Monogr. 87, 470–486 (2017).
gaps that may imprint clumped distributions of recruits of pioneer species40. For  13. Serván, C. A., Capitán, J. A., Grilli, J., Morrison, K. E. & Allesina, S. 
each species we assume a density λfc of randomly distributed cluster centres, which Coexistence of many species in random ecosystems. Nat. Ecol. Evol. 2, 
have, at each time step, a probability pfp of changing location. For each of these 1237–1242 (2018).
recruits, we first randomly select one cluster centre among the cluster centres of the  14. Kraft, N. J. B., Godoy, O. & Levine, J. M. Plant functional traits and the 
corresponding species and then determine the position of the recruit by sampling multidimensional nature of species coexistence. Proc. Natl Acad. Sci. USA 
from the kernel. For the simulation shown in Extended Data Fig. 5a, the recruits 112, 797–802 (2015).
were located at random positions within the plot.  15. Allesina, S. & Tang, S. Stability criteria for complex ecosystems. Nature 483, 
205–208 (2012).
Parameterization of the simulation model. Extended Data Fig. 5 shows  16. Barbier, M., Arnoldi, J.-F., Bunin, G. & Loreau, M. Generic assembly patterns 
simulations of the individual-based model conducted in a 200 ha area containing in complex ecological communities. Proc. Natl Acad. Sci. USA 115, 
approximately 83,000 trees with, initially, 80 species. There was no immigration. 2156–2161 (2018).
The model parameters were the same for all species, and all species followed exactly  17. Detto, M. & Muller-Landau, H. C. Stabilization of species coexistence in 
the same model rules. We selected βfi = βff to obtain no differences in con- and 
spatial models through the aggregation–segregation effect generated  
heterospecific interactions and s  = 1 (no background mortality), and we adjusted the by local dispersal and nonspecific local interactions. Theor. Pop. Biol. 112, f
parameters βff = 0.0075 and rf = 0.1 to yield tree densities (415 ha−1) and an overall 5 
97–108 (2016).
yr mortality rate (10%) similar to those of trees with dbh ≥ 10 cm in the BCI plot51.  18. Bolker, B. & Pacala, S. W. Spatial moment equations for plant competition: 
The Gaussian kernel used to place recruits around conspecific adults or around understanding spatial strategies and the advantage of short dispersal. Am. 
random cluster centres had a parameter σ = 10 m. There were 40 random cluster Nat. 153, 575–602 (1999).
centres in total for each species that had a probability of pfp = 0.3 of changing 
 19. Uriarte, M. et al. Trait similarity, shared ancestry and the structure of 
location within one census interval. The only difference between the simulation neighbourhood interactions in a subtropical wet forest: implications for 
shown in Extended Data Fig. 5b and the one shown in Extended Data Fig. 5c is community assembly. Ecol. Lett. 13, 1503–1514 (2010).
that in the former, we used a proportion p  = 0.05 of recruits to be placed around  20. Wang, X. et al. Stochastic dilution effects weaken deterministic effects of d
randomly distributed cluster centres (that is, 95% of the recruits were placed close niche-based processes on the spatial distribution of large trees in species rich 
to their parents), but in the latter, we selected pd = 0.95 (that is, 95% of the recruits 
forests. Ecology 97, 347–360 (2016).
were placed around randomly distributed cluster centres). In our simulations, on  21. Wiegand, T. et al. Spatially explicit metrics of species diversity, functional 
average, one of these cluster centres received four recruits per time step, which diversity, and phylogenetic diversity: insights into plant community assembly 
were scattered within a radius of approximately 30 m, and received approximately processes. Annu. Rev. Ecol. Syst. 48, 329–351 (2017).
13 recruits during its lifetime (at each time step it had a probability of 0.3 of  22. Godoy, O. & Levine, J. M. Phenology effects on invasion success: insights 
changing location). In contrast, in Extended Data Fig. 5a recruits were placed at from coupling field experiments to coexistence theory. Ecology 75,  
random locations within the plot. 726–736 (2014).
 23. Dieckmann, U. & Law, R. in The Geometry of Ecogical Interactions: Simplifying 
Spatial Complexity (eds. Dieckmann, U. et al.) Ch. 21 (Cambridge Univ. 
Reporting Summary. Further information on research design is available in the Press, 2000).
Nature Research Reporting Summary linked to this article.  24. Klausmeier, C.A. & Tilman, D. in Competition and Coexistence (eds. Sommer, 
U. & Worm, B.) Ch. 3 (Springer, 2002).
Data availability  25. Fahse, L., Wissel, C. & Grimm, V. Reconciling classical and individual-based 
The data that support the findings in this manuscript (and the raw data for  approaches in theoretical population ecology: a protocol for extracting 
Figs. 2–4 and Extended Data Figs. 2–4 and 8) can be found in Supplementary population parameters from individual-based models. Am. Nat. 152,  
Data Table 1. To generate this data, we used the raw census data of the ForestGEO 838–852 (1998).
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
 26. Fortunel, C., Valencia, R., Wright, S. J., Garwood, N. C. & Kraft, N. J. B.  51. Wiegand, T., May, F., Kazmierczak, M. & Huth, A. What drives the spatial 
Functional trait differences influence neighbourhood interactions in a distribution and dynamics of local species richness in tropical forests. Proc. R. 
hyperdiverse Amazonian forest. Ecol. Lett. 19, 1062–1070 (2016). Soc. B 284, 20171503 (2017).
 27. Cadotte, M. W., Davies, T. J. & Peres-Neto, P. R. Why phylogenies do not 
always predict ecological differences. Ecol. Monogr. 87, 535–551 (2017). Acknowledgements
 28. Wiegand, T. & Moloney, K. A. A Handbook of Spatial Point Pattern Analysis X.W. was supported by the Strategic Priority Research Program of the Chinese Academy 
in Ecology (CRC Press, 2014). of Sciences (Grant XDB31030000), the National Natural Science Foundation of China 
 29. Illian, J., Penttinen, A., Stoyan, H. & Stoyan, D. Statistical Analysis and (Grant 31961133027), the Key Research Program of Frontier Sciences, Chinese Academy 
Modelling of Spatial Point Patterns (Wiley, 2008). of Sciences (Grant ZDBS-LY-DQC019) and the K.C. Wong Education Foundation. 
 30. Wang, S. Simplicity from complex interactions. Nat. Ecol. Evol. 2,  T.W. and A.H. were supported by the ERC advanced grant 233066. M.C. and X. M were 
1201–1202 (2018). supported by the National Natural Science Foundation (32061123003 and 31770478). 
 31. Wilson, W. G. et al. Biodiversity and species interactions: extending L.L. was supported by the Joint Fund of the National Natural Science Foundation of 
Lotka–Volterra community theory. Ecol. Lett. 6, 944–952 (2003). China-Yunnan Province (U1902203). The FS plot project was supported by the Taiwan 
 32. Hubbell, S. P. Neutral theory and the evolution of functional equivalence. Forestry Bureau, the Taiwan Forestry Research Institute and the Ministry of Science 
Ecology 87, 1387–1398 (2006). and Technology. The BCI censuses have been made possible through support of the US 
 33. Rosindell, J., Hubbell, S. P. & Etienne, R. S. The unified neutral theory of National Science Foundation (awards 8206992, 8906869, 9405933, 9909947, 0948585 to 
biodiversity and biogeography at age ten. Trends Ecol. Evol. 26,  S.P. Hubbell), the John D. and Catherine D. McArthur Foundation and the Smithsonian 
340–348 (2011). Tropical Research Institute. We also thank the hundreds of people who contributed to the 
 34. Chesson, P. Scale transition theory: its aims, motivations and predictions. collection and management of the data from the plots. This work builds on discussion 
Ecol. Complex. 10, 52–68 (2012). of the working group sNiche (Expanding neo-Chessonian coexistence theory towards a 
 35. May, F., Wiegand, T., Lehmann, S. & Huth, A. Do abundance  stochastic theory for species-rich communities) supported by sDiv, the Synthesis Centre 
distributions and species aggregation correctly predict macroecological of iDiv (DFG FZT 118). We thank S. Lehmann, D. Alonso and S. Harpole for discussion, 
biodiversity patterns in tropical forests? Glob. Ecol. Biogeogr. 25,  and especially S. Stump for valuable feedback on earlier drafts.
575–585 (2016).
 36. Murrell, D. When does local spatial structure hinder competitive coexistence 
and reverse competitive hierarchies? Ecology 91, 1605–1616 (2010). Author contributions
 37. May, F., Wiegand, T., Huth, A. & Chase, J. M. Scale-dependent effects of T.W. and X.W. conceived and designed the project. T.W. implemented the simulation 
conspecific negative density dependence and immigration on biodiversity model, conducted the simulations, analysed the results and prepared figures and tables. 
maintenance. Oikos 129, 1072–1083 (2020). T.W. and A.H. led the writing of the manuscript. X.W. assembled and analysed the 
 38. Rosindell, J. & Cornell, S. J. Species–area relationships from a spatially plot data and conducted the spatial analyses. K.J.A.-T., N.B., M.C, X.C., S.J.D, Z.H., 
explicit neutral model in an infinite landscape. Ecol. Lett. 10,  R.H., W.J.K., J. Lian, J. Li, L.L., Y.L., K.M., W.M., X.M., S.-H.S., I.-F.S., A.W., X.W. and 
586–595 (2007). W.Y. contributed to the acquisition of the data used in the paper and in revising the 
 39. Chave, J. & Leigh, G. L. A spatially explicit neutral model of betadiversity in manuscript. All authors have given final approval to publish this manuscript and agree to 
tropical forests. Theor. Pop. Biol. 62, 153–168 (2002). be accountable for the aspects of the work that they conducted.
 40. Hubbell, S. P. et al. Light-gap disturbances, recruitment limitation, and tree 
diversity in a neotropical forest. Science 283, 554–557 (1999). Competing interests
 41. Chanthorn, W., Getzin, S., Wiegand, T., Brockelman, W. Y. & Nathalang, A. The authors declare no competing interests.
Spatial patterns of local species richness reveal importance of frugivores for 
tropical forest diversity. J. Ecol. 106, 925–935 (2018).
 42. Getzin, S., Wiegand, T. & Hubbell, S. P. Stochastically driven adult-recruit Additional information
associations of tree species on Barro Colorado Island. Proc. R. Soc. B 281, Extended data is available for this paper at https://doi.org/10.1038/s41559-021-01440-0.
20140922 (2014). Supplementary information The online version contains supplementary material 
 43. Chesson, P. & Neuhauser, C. Intraspecific aggregation and species available at https://doi.org/10.1038/s41559-021-01440-0.
coexistence. Trends Ecol. Evol. 17, 210–211 (2002).
 44. Ruokolainen, L. & Hanski, I. Stable coexistence of ecologically identical Correspondence and requests for materials should be addressed to T.W. or X.W.
species: conspecific aggregation via reproductive interference. J. Anim. Ecol. Peer review information Nature Ecology & Evolution thanks Simon Stump and the 
85, 638–647 (2016). other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer 
 45. Stump, S. Multispecies coexistence without diffuse competition; or, why reviewer reports are available.
phylogenetic signal and trait clustering weaken coexistence. Am. Nat. 190, Reprints and permissions information is available at www.nature.com/reprints.
213–228 (2017).
 46. Harms, K. E., Condit, R., Hubbell, S. P. & Foster, R. B. Habitat associations  Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 
of trees and shrubs in a 50-ha neotropical forest plot. J. Ecol. 89,  published maps and institutional affiliations.
947–959 (2001). Open Access This article is licensed under a Creative Commons 
 47. Wiegand, T. et al. Testing the independent species’ arrangement assertion Attribution 4.0 International License, which permits use, sharing, adap-
made by theories of stochastic geometry of biodiversity. Proc. R. Soc. B 279, tation, distribution and reproduction in any medium or format, as long 
3312–3320 (2012). as you give appropriate credit to the original author(s) and the source, provide a link to 
 48. Wiegand, T., Gunatilleke, C. V. S., Gunatilleke, I. A. U. N. & Huth, A. How the Creative Commons license, and indicate if changes were made. The images or other 
individual species structure diversity in tropical forests. Proc. Natl Acad. Sci. third party material in this article are included in the article’s Creative Commons license, 
USA 104, 19029–19033 (2007). unless indicated otherwise in a credit line to the material. If material is not included in 
 49. Erickson, D. L. et al. Comparative evolutionary diversity and phylogenetic the article’s Creative Commons license and your intended use is not permitted by statu-
structure across multiple forest dynamics plots: a mega-phylogeny approach. tory regulation or exceeds the permitted use, you will need to obtain permission directly 
Front. Genet. 5, 358 (2014). from the copyright holder. To view a copy of this license, visit http://creativecommons.
 50. Howe, H. F. Scatter- and clump-dispersal and seedling demography: org/licenses/by/4.0/.
hypothesis and implications. Oecologia 79, 417–426 (1989). © The Author(s) 2021
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles 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 nkff over species. b. The mean n̄fh of the heterospecific neighbourhood crowding index nkfh over species. c. The 
mean n̄fβ of the interaction neighbourhood crowding index nkfβ over species. d. The variance-to-mean ratio bkf of the conspecific neighbourhood crowding 
index nkff over species. e. The variance-to-mean ratio bkh of the heterospecific neighbourhood crowding index nkfh over species. f. The variance-to-mean 
ratio bkβ of the interaction neighbourhood crowding index nkfβ 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.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
Extended Data Fig. 3 | Correlation between different crowding indices. We estimated for all individuals of a species f the correlation between their 
crowding indices nkff and nkfβ (a) and nkfh and nkfβ (b). The crowding index nkff counts the conspecific neighbours of individual k within distance R = 10 m, 
nkfh counts the corresponding number of heterospecific neighbours, and nkfβ 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 eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles NAturE EcoLogy & EvoLutioN
Extended Data Fig. 4 | The distribution of the relative population-level interaction coefficients. a, The distribution of the relative population level 
interaction coefficients αfi/αff for the focal species of the tropical forests, resulting from equation 6. b, same as a, but for subtropical forests. c, same as a, 
but for temperate forests. Other conventions as in Fig. 3.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
Extended Data Fig. 5 | Spatial patterns can stabilize community dynamics of symmetric species. Individual-based simulations of a symmetric model 
with initially 80 identical species, simulated on an area of 200 ha without immigration for 5000 time steps (25,000 years). Different colours correspond 
to different species. left, Recruits were randomly distributed, the dynamics is unstable with 2 extinctions. middle, Recruit were mostly scattered around 
conspecific adults, the dynamics is unstable with 7 extinctions. right, recruits were mostly scattered around random cluster centre, the dynamics is stable 
without extinctions. a –c, Species abundances with the bold black line indicating the expected mean abundance J*/80. d–f, Intraspecific pattern kff. g–i, 
Interspecific pattern kfh of heterospecifics with respect to the foal species f. j–l, the mean relative interaction strength Bf of a heterospecific neighbour 
of an individual of species f. m–o, Stabilization, being the population-level heterospecific interaction strength relative to the corresponding conspecific 
interaction strength (that is, αfh/αff). Note the different scale of the y-axes. Spatial patterns were measured at a 10 m neighbourhood.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles NAturE EcoLogy & EvoLutioN
Extended Data Fig. 6 | The relationship between clustering and abundance. Shown are detailed results of the simulation of Extended Data Fig. 5b for two 
species that went extinct. a, Dynamics of the relationship between abundance Nf(t) and clustering kff for species 1. b, same for species 3. The solid line 
indicates the power law kff = 10,000/Nf(t).
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
Extended Data Fig. 7 | Mechanism underlying the fitness-density covariance in our simulation model. Quantities involved in the estimation of the 
fitness–density covariance (equation 21) for model simulations of the stable scenario (left) and the unstable scenario (right). The data were taken 
from time step 5000 of the simulations shown in Extended Data Fig. 5c (stable dynamics) and of Extended Data Fig. 5b (unstable dynamics). Panels 
show the mean number n̄ff + n̄fh of all neighbours (a, b), of conspecific neighbours n̄ff  (c, d), of heterospecific neighbours n̄fh (e, f), and the resulting 
fitness-density covariance λ̃f − λ̄f  (g, h), of all species plotted over their abundance. The red lines show the expected relationship based on equations (21) 
and (22). We fitted power law relationships kff(Nf) = a N bf  to the data. For stable dynamics we find kff(Nf) = 8.2, kfh(Nf) = 0.905 and for unstable dynamics 
kff(N ) = 3723/N 0.883f ff  and kfh (Nf) = 0.890.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
Articles NAturE EcoLogy & EvoLutioN
Extended Data Fig. 8 | See next page for caption.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
NAturE EcoLogy & EvoLutioN Articles
Extended Data Fig. 8 | Rare species advantage and abundance −clustering relationships in model simulations and ForestgeO plots. The covariance 
between the local dominance dkff of species f (that is, dkff = nkff/(nkff + nkfh)) and the total number of neighbours (that is, nkff + nkfh). A positive relationship 
of the covariance with abundance indicates that, when a species becomes rare, areas of higher conspecific crowding have fewer competitors. Shown are 
results of model simulations (a, b) and the nine ForestGEO forest plots (c–k). The p-value is for the null hypothesis that the slope of the linear regression 
(red lines) is not positive. (l–v) The corresponding relationships between clustering kff and abundance Nf. The value of b is the slope of the power law 
kff(Nf) = a N bf  (red line). We used for the analysis focal species f with more than 50 individual. For plot names see Supplementary Table 1, and for raw data 
and details of the linear regressions see Supplementary Data Table 1.
NATuRe eCOLOgY & evOLuTION | www.nature.com/natecolevol
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 
Give P values as exact values whenever suitable.
For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes
Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated
Our web collection on statistics for biologists contains articles on many of the points above.
Software and code
Policy information about availability of computer code
Data collection The study did not involve field data collection. We wrote Pascal (Dephi 5.0) code for the simulation model to conduct the example simulations 
(results reported in Extended Figures 5, 6, 7, and Supplementary Figures S1-S5). The code of the simulation model  is provided as 
supplementary file.
Data analysis The code for estimating the crowding indices and other indices of spatial patterns from the ForestGeo data (results reported in Figs. 1, 2, 
Extended Data Figures 1, and 2) is included in the code of the simulation model.  The same procedures can be used for simulated data and 
census data. The simple analyses  for Figs. 2, 3, 4 and Extended Data Figures  2, 3, 4 and 8 are provided as Supplementary Data 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 . 
1
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.
Life sciences Behavioural & social sciences  Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf
Ecological, evolutionary & environmental sciences study design
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 fully-
mapped 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.
Randomization Not applicable to this study.
Blinding Blinding was not applicable to this study. 
Did the study involve field work? Yes No
Reporting for specific materials, systems and methods
We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, 
system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. 
Materials & experimental systems Methods
n/a Involved in the study n/a Involved in the study
Antibodies ChIP-seq
Eukaryotic cell lines Flow cytometry
Palaeontology and archaeology MRI-based neuroimaging
Animals and other organisms
Human research participants
Clinical data
Dual use research of concern
2
nature research  |  reporting summary April 2020