## Abstract

Understanding the structure and dynamics of highly diverse tropical forests is challenging. Here we investigate the factors that drive the spatio-temporal variation of local tree numbers and species richness in a tropical forest (including 1250 plots of 20 × 20 m^{2}). To this end, we use a series of dynamic models that are built around the local spatial variation of mortality and recruitment rates, and ask which combination of processes can explain the observed spatial and temporal variation in tree and species numbers. We find that processes not included in classical neutral theory are needed to explain these fundamental patterns of the observed local forest dynamics. We identified a large spatio-temporal variability in the local number of recruits as the main missing mechanism, whereas variability of mortality rates contributed to a lesser extent. We also found that local tree numbers stabilize at typical values which can be explained by a simple analytical model. Our study emphasized the importance of spatio-temporal variability in recruitment beyond demographic stochasticity for explaining the local heterogeneity of tropical forests.

## 1. Background

The high diversity of tropical forests is one of the enduring puzzles in ecology [1,2] that challenges ecological theories. Previous modelling studies seeking to understand the basic dynamics of tropical forests focused on non-spatial patterns (e.g. [3–5]) and explored equilibrium properties such as the species abundance distribution or the species–area relationship (e.g. [6,7]). Prominent examples are neutral models that consider dispersal limitation and ecological drift, and assume that species are functionally equivalent [3,8]. They provide realistic predictions of different biodiversity patterns [4,7,9]. Extended models have also been used to explain spatio-temporal patterns of tropical forests such as species extinction and temporal shifts in species clusters, or to assess environmental variation and temporal storage effects [5,10–12].

Here we aim to investigate the mechanisms behind two spatio-temporal patterns that have rarely been studied with dynamic models: the variation in local tree numbers and in local species richness found in tropical forests [13–15]. For example, tree species richness in the 50 ha forest dynamics plot on Barro Colorado Island (BCI; Panama) ranged from approximately 20 to 80 species at the local scale (20 × 20 m^{2}; figure 1*a,b*). We argue that these patterns can provide important additional information on the spatial and temporal dynamics of species rich communities that has not been investigated until now.

We propose here a novel approach to explore the spatial and temporal dynamics of species rich forests using simple simulation models that are built around the local spatial variation of demographic rates. A specific feature of our approach is inclusion of three spatial scales: the 20 × 20 m^{2} quadrat scale (where forest state and dynamics are analysed), the 50 ha community scale (the extent of our spatially explicit simulation models), and the regional scale represented by a simple metacommunity. Our approach allows us to test different biological hypotheses. We accomplish this by combining field data analysis with a dynamic null model approach that simulates recruitment and mortality in each quadrat by maintaining certain aspects of the observed distributions (of the number of recruits per quadrat and the mortality rates per quadrat, respectively) while randomizing other aspects. In total, we test six simulation models that make alternative assumptions on local mortality and recruitment in their ability to explain the observed spatio-temporal variation in local tree and species numbers. We underpin our simulation results by an analytical model for local tree numbers and additional data analyses.

Our simplest hypothesis is that variations in local tree numbers and in local species richness emerge solely as a consequence of random birth and death processes [3,6] (i.e. the demographic stochasticity hypothesis). However, additional mechanisms not considered in classical neutral theory may contribute to the local dynamics of tropical forests. We focus here on the additional contribution of spatio-temporal variability in demographic rates. One possible mechanism is that local clusters of high tree density (and species richness) may be related to locally favourable environmental conditions for recruitment and/or survival (habitat hypothesis) [16]. Alternatively, heterogeneity in recruitment (or mortality) may be caused by spatially and temporally varying biotic habitat conditions at the local scale (e.g. temporal occurrence of tree fall gaps [15,17,18]) that might result in randomly moving local clusters of high tree density and species richness (moving cluster hypothesis) [11]. A similar effect for recruitment would be expected in the case of strong spatio-temporal variation in seed production [19]. Within the framework of our approach we can analyse such different hypothesis.

To test the performance of the alternative model versions we assess how well they reproduced several spatial and temporal patterns observed at the 50 ha plot of neo-tropical forest BCI in Panama [20]. Our primary focus is on the observed quadrat-scale distributions of tree abundance *N* and species richness *S* per quadrat. As additional spatio-temporal patterns we test, for example, the changes in the number of trees and species per quadrat and time (Δ*N* versus *N*_{0} and Δ*S* versus *S*_{0}), and to explain the empirical pattern of Δ*N* versus *N*_{0} we use an analytical model. We show that the dynamics of local tree density were stabilized mainly by local recruitment being largely independent of the local tree density. The observed large spatio-temporal variability in local recruitment was the key factor for emergence of the observed variation in local tree density, and limited dispersal was important for reproduction of the observed pattern of local species richness.

## 2. Material and methods

### (a) Study site

We used data from the 50 ha (1000 × 500 m^{2}) forest dynamics plot on BCI [15,20,21]. The plot is located on an island in the Panama Canal (9°15′ N, 79°86′ W) and the vegetation on BCI represents seasonally moist tropical forest with a pronounced dry season and mean annual rainfall of 2600 mm. The plot comprises old growth lowland moist forest and contains up to 305 tree and shrub species. Elevation in the plot ranges from 120 to 155 m.a.s.l. All trees with a stem diameter at breast height (dbh) of 1 cm or more have been mapped and had their diameters measured in 1982 and 1985, and from then on every five years until 2010, resulting in seven censuses.

### (b) Data

We used the six censuses from 1985 to 2010. As recruits we classified all trees which appeared for the first time in a given census (i.e. grown during the last 5 years to a dbh ≥ 1 cm). As dead trees we classified those trees that were alive in a given census and recorded as dead or disappeared during the next one.

For data analysis we divided the forest area into 1250 non-overlapping contagious quadrats with a size of 20 × 20 m^{2}. However, to test the generality of our approach we also conducted supplementary analyses using quadrat sizes of 10 × 10 m^{2} and 50 × 50 m^{2}.

For every quadrat and census we counted the number of living trees and the number of species, and determined their distribution (figure 1). We also determined the distributions of the numbers of recruits and of dying trees per quadrat and census interval (figure 2*c,d*). Finally, we determined the distribution of the *per capita* mortality rates within quadrats by dividing in each quadrat the number of dying trees between census *t* and *t* + 1 by the number of trees at census *t*. Electronic supplementary material, figure C1 shows all distributions and pairwise correlations among these variables (using rank-correlation analysis).

To assess the temporal variability in several quadrat-scale variables (e.g. number of recruits or *per capita* mortality rates) we estimated the coefficients of variation (CV) of these variables for each quadrat across the available censuses and then determined the distribution of the CV values of all 1250 quadrats. We also determined the analogous distribution for spatial variability.

## 3. Simulation models

### (a) General structure of the simulation model at the 50 ha community scale

We developed several simple, but spatially explicit simulation models describing the community dynamics at the 50 ha scale based on individual trees. All models include stochastic representations of mortality, dispersal, recruitment and immigration of trees with dbh ≥ 1 cm in a 50 ha (1000 × 500 m^{2}) forest plot. Each tree is characterized by its position in space (continuous *x*, *y* coordinates) and its species identity. For simplicity, tree size is not considered. In each time step mortality, recruitment and immigration events are calculated.

As simple reference model we used a spatially explicit simulation model, which corresponds to a spatially explicit implementation of neutral theory [3]. The number of trees that die per time step (i.e. the 5-year census interval) is calculated based on a *per capita* mortality rate, which was directly calculated from the data (*m* = 0.15 per 5 years). When a tree dies, a new tree (recruit) is added to the community (i.e. the zero-sum assumption [3]). For determining the position of the recruit, we first randomly select a parent tree among all trees at the 50 ha community scale and then apply a dispersal kernel. In the more complex model versions we make alternative assumptions on mortality and recruitment (see ‘Submodels for hypothesis testing’). To investigate the consequence of dispersal, we use an isotropic Gaussian dispersal kernel [22,23] with specified mean dispersal distance (10, 20, 30, 40 and 50 m). We use periodic boundary conditions to avoid edge effects. In addition, we simulate a long-distance dispersal (LDD) scenario, where we select a random position for the recruit and a random parent anywhere in the plot.

Similar to previous models (e.g. [4]), we also considered a low probability that a recruit is an immigrant of a new species. The probability was set to 1.1 × 10^{−4} per recruitment event. This parameter reflects the observed rate of species immigrations of approximately four species per 5 years [24] and represents a simplified metacommunity. The model is written in C++ and source code for all model versions is published on the online repository GitHub (https://github.com/FelixMay/confettiSquares).

### (b) Submodels for hypothesis testing

We use three alternative submodels for recruitment (R1, R2 and R3) and two alternative submodels for mortality (M1 and M2), summarized in table 1. Other possible submodels could be excluded *a priori* by analysis of our data. We investigate all six combinations together with the six alternative assumptions on dispersal, resulting in a total of 3 × 2 × 6 = 36 alternative simulation models of spatial and temporal forest dynamics to explain the spatial variation in local tree density and species richness.

To test the demographic stochasticity hypothesis, we use submodels for recruitment (R1) and mortality (M1) based on constant *per capita* mortality and recruitment rates *m* and *r*, respectively. This model version corresponds to spatially explicit neutral models [7,25,26].

To assess the moving cluster hypothesis, we test if the observed spatio-temporal variability in recruitment contributed to the emergence of our focal patterns. We determine for each quadrat *q* the absolute number *R _{q}* of recruits by taking randomly a value from the distribution

*p*

_{rec}(

*R*) observed in the field for all six censuses together (submodel R2). Because we draw each time step a new value of

*R*from

_{q}*p*

_{rec}(

*R*), we assume here equivalence of spatial and temporal variability in

*R*. Because in this model version the number of recruits in a given quadrat is set, we determine the species identity of the recruit by ‘backward’ application of the dispersal kernel (see electronic supplementary material, appendix A). In this submodel we expect emergence of higher variability in the recruit numbers per quadrat (compared with submodel R1). To assess the relative contribution of the observed spatial variability in local mortality rates we draw in submodel M2 for each quadrat a local

_{q}*per capita*mortality rate

*m*from the observed distribution [

_{q}*p*

_{mort}(

*m*)] (similar approach as for R2).

To test for possible density dependence in recruitment we extend submodel R2 by additionally maintaining the observed covariation between the number *R _{q}* of recruits and the number

*N*of trees in quadrat

_{q}*q*(submodel R3). A detailed description of the simulation model is provided in the electronic supplementary material, appendix A.

We also test an extreme case of the habitat hypothesis for local tree numbers where *R _{q}* and

*m*are temporally constant but vary in space following the observed distributions

_{q}*p*

_{rec}(

*R*) and

*p*

_{mort}(

*m*), respectively. This model can be solved without simulation (electronic supplementary material, appendix B).

_{q}### (c) Model initialization, scheduling and simulation scenarios

All model simulations were initialized with the observed BCI census of 1995 (in total 228 819 trees, stem diameter > 1 cm). For each parameter set we conducted 30 replicate model runs. Each simulation was run for 5025 years (1005 time steps, one time step corresponds to 5 years), and the last six simulated time steps were analysed and compared with the field data. We checked that this simulation time was sufficient to reach a dynamic equilibrium for the patterns investigated in this study.

### (d) Comparison of field data and model predictions

We compared the output of our alternative models with a number of observed patterns. They included (i) the observed distributions of the number of species or trees per quadrat (figure 1*b,c*), (ii) the correlations between the 25-year change in the number of species (and trees) per quadrat and the initial numbers (i.e. the 1985 and 2010 census from the data and time steps 1000 and 1005 from the simulations) (figure 3), (iii) the correlations between the number of trees and dead trees, between the number of trees and the mortality rate, between the number of trees and recruits and between the number of dead trees and recruits, and (iv) the correlation between the number *N* of trees and the number *S* of species per quadrat (electronic supplementary material, figure C1).

Because our simulation models are rather simple (they can be seen as null models) we do not expect perfect fits of these patterns, but comparison of results allows us to infer on the aspects of recruitment and mortality that are critical for tropical forest dynamics. We compared the mean and the standard deviations between the simulated and observed distributions and applied Kolmogorov–Smirnov tests to quantitatively check their match.

## 4. Analytical model for the dynamics of tree numbers at the local scale

We develop an analytical model to explain the empirical pattern of the change in the number of trees per quadrat and time step relative to the initial number of trees (i.e. Δ*N* versus *N*_{0}) (figure 3*a*). In classical population ecology [27] this type of relationship is frequently used and it allows us to determine the stability properties of equilibriums *N** at the local quadrat scale (see electronic supplementary material, appendix B).

First, we assume similar to our simulation model M1 + R1 constant *per capita* mortality and reproduction rates *r* and *m*, respectively, and obtain
4.1

This model assumes that all recruits stem from parents in the same quadrat (i.e. strong dispersal limitation) and shows unstable dynamics. However, in our spatially explicit simulations only a certain proportion *c* of the offspring may remain in the home quadrat of their parents (i.e. a term + *crN _{t}* with

*c*< 1) whereas the rest is dispersed into other quadrats. As a consequence, each quadrat will receive an additional input of recruits

*R*from other quadrats that is independent on the local number

*N*of trees. We assume here for simplicity that

_{t}*R*does not vary in time. A generalization of model (4.1) is therefore given by 4.2

The equilibrium of equation (4.2) is given by . Additionally, at equilibrium the total number of recruits must equal the mean number of recruits per quadrat taken over all quadrats at the 50 ha community scale, such that 4.3This results in an estimate of that can be compared with the proportion of recruits remaining in the parent quadrat for different mean dispersal distances (electronic supplementary material, figure B1).

The analytical model of equations (4.2) and (4.3) shows a stable equilibrium if *c* < 1. This is a consequence of the recruitment input *R* that operates independently on local tree numbers and ‘rescues’ quadrats with low tree numbers from local extinction in a way known from metapopulation theory [28].

## 5. Results

### (a) Demographic stochasticity hypothesis

Our reference model (R1 + M1), based only on demographic stochasticity and dispersal limitation, substantially underestimated the variation in the local number of trees and species (figure 2*a,b*; electronic supplementary material, figure C3). A mean dispersal distance of 20 m produced the best agreement in the spatial distribution of local species richness (figure 2*b*; electronic supplementary material, figure C3). The assumption of constant *per capita* mortality and recruitment rates led to a substantial underestimation of the local variability in recruitment (figure 2*d*) and to a moderate underestimation of the local number of dead trees (figure 2*c*). Thus, we find strong spatial variability in local recruitment and a moderate spatial variability in mortality.

Our analytical model (equations (4.2) and (4.3)) approximates the relationship between Δ*N* (the 25 year change in the number of trees per quadrat) and *N*_{0} (the initial number of trees in 1985) well (figure 3*a*; *R*^{2} = 0.45). It allows us also to assess the proportion *c* of recruits that remain in the parent quadrat (; electronic supplementary material, appendix B), which corresponds to a mean dispersal distance of the Gaussian kernel of approximately 40 m (figure B1). Thus, recruitment at the local quadrat scale is largely independent of the local number of trees.

The analytical model (equations (4.2) and (4.3)) predicts an increase in the correlation between Δ*N* and *N*_{0} with increasing mean dispersal distance. The reference simulation model (R1 + M1) reproduced this pattern (electronic supplementary material, figure C4), but the field data show stronger correlation between Δ*N* and *N*_{0} than predicted by R1 + M1 (figure 3*a,c*).

### (b) Spatio-temporal variability in recruitment

Recruitment submodel R2 reproduces the observed variation in the local number of trees and species well when combined with a *per capita* mortality rate (i.e. submodel M1 for a mean dispersal distance of 20 m; figure 2*e,f*; electronic supplementary material, figure C3). Thus, spatio-temporal variability in the number of recruits beyond demographic stochasticity is a key factor in generating the observed variation in the local number of trees and species.

Further data analysis of the spatial and temporal variability in the local number of recruits provides hints on the underlying mechanism. We found very high temporal variability in local recruitment among the five censuses with a mean CV of 0.63 and 2.5%, and 97.5% percentiles of 0.23 and 1.27, respectively (electronic supplementary material, figure C5*a*). The temporal variability was basically identical to the analogous spatial variability (electronic supplementary material, figure C5*b*), thus making the moving cluster hypothesis more likely for recruitment than the habitat hypothesis. Additionally, implementation of the habitat hypothesis without temporal variability produced too wide a variation in local tree numbers (electronic supplementary material, figure B2*b*).

### (c) Spatio-temporal variability in mortality

We find a positive correlation between the number of dead trees and the number of trees in a quadrat (electronic supplementary material, figure C1; *r* = 0.62), which suggests that mortality can be described in first approximation by a constant *per capita* rate. Consequently, we find that model R2 + M1 with a constant *per capita* mortality rate *m* approximates the distribution in the number of dead trees per quadrat fairly well (figure 2*g*).

However, model R2 + M1 produced too high a correlation between the observed numbers of trees and dead trees per quadrat (electronic supplementary material, figure C6*c*), suggesting that the spatial variability in local mortality cannot be ignored. Indeed, mortality submodel M2 that determines *m _{q}* by random draws from the observed distribution

*p*

_{mort}(

*m*) improved the (already reasonable) results when combined with recruitment submodel R2 (cf. figure 2

*g*and figure 2

*k*), and also produced a good match of the correlation between the observed numbers of trees and dead trees per quadrat (figure C6

*b*). However, model R2 + M2 does not improve the prediction of the distribution of the number of trees and species in quadrats (figure 2

*e,f*versus

*i,j*) and the correlations between Δ

*N*and

*N*

_{0}and Δ

*S*and

*S*

_{0}(electronic supplementary material, figure C4) when compared with model R2 + M1.

The spatial and temporal variability in local mortality rates was substantially smaller than that of the number of recruits (electronic supplementary material, figure C5*a–d*). We also found that the spatial variability was very similar to the corresponding temporal variability in mortality (electronic supplementary material, figure C5*c,d*).

### (d) The density dependence hypothesis

The observed correlation between the local number of trees and the local number of recruits was low (electronic supplementary material, figure C1), and maintaining the observed covariation between number of trees and number of recruits in a quadrat (i.e. recruitment submodel R3) did not lead to any improvements with respect to submodel R2 (electronic supplementary material, figure C2). Similarly, the correlation between the *per capita* mortality rate and the number of trees in a quadrat was low (electronic supplementary material, figure C1; *r* = 0.16). Thus, density dependence of demographic rates at the community level (without regard to species) had no influence on the emergence of our focal patterns.

### (e) Results with different quadrat sizes

We repeated all simulation analyses with quadrat sizes of 10 × 10 m^{2} and 50 × 50 m^{2} (see electronic supplementary material, figures C7 and C8) and found that all main findings are robust to these changes. There was only a slight change of the optimal dispersal distance that is needed to predict the variation in local species richness (between 10 m and 20 m for 10 × 10 m^{2} and approx. 40 m for 50 × 50 m^{2} quadrats). Our findings with respect to the recruitment and mortality submodels were not affected by quadrat size.

## 6. Discussion

In this study, we provide insights into the mechanisms that drive the spatial and temporal variation of local tree density and species richness in a neo-tropical forest. For this purpose, we tested a number of biological hypotheses by using dynamic models that are built around the local spatial variation of mortality and recruitment (extracted from the census data by using 1250 non-overlapping 20 × 20 m^{2} quadrats). We then asked which of the hypotheses (i.e. model versions) are able to generate the observed spatial and temporal variation of local tree numbers, species richness and other spatio-temporal patterns extracted from the data.

### (a) Biological hypotheses

We find that our two main patterns, the variation in local tree numbers and species richness, do not emerge solely as a consequence of stochastic birth and dead processes and dispersal limitation as assumed by our simple reference model R1 + M1 (the demographic stochasticity hypothesis). Thus, additional processes not included in classical neutral theory are operating at the BCI forest in generating these patterns. We identified a large spatio-temporal variability in the local number of recruits (that operated independently on local tree density) as the main missing mechanism to reproduce our main patterns, whereas variability in local *per capita* mortality rates contributed to a lesser extent.

There are three general possibilities for the origin of the excessive variability in local recruitment: it can be spatial, temporal or spatio-temporal. Environmental stochasticity alone (e.g. [5,12]) (i.e. the second possibility) can be excluded because spatial variability in recruitment was large for individual censuses (electronic supplementary material, figure C9). One obvious hypothesis that corresponds to the first possibility is that differences in local habitat conditions would generate spatial variability in local recruitment (the habitat hypothesis). Indeed, several studies at BCI found associations of the distribution of species with habitat features such as topography or soil attributes (e.g. [16,29,30]). However, the habitat hypothesis did not fit the observed distribution of local tree numbers (electronic supplementary material, figure B2*b*), and we found that spatial and temporal variations in local recruit numbers show similar magnitudes (electronic supplementary material, figure C5*a,b*). This rather points to spatio-temporally varying biotic conditions as the underlying mechanism (i.e. the moving cluster hypothesis). Tree fall gaps, for example, are able to generate locally high recruit numbers by making light available that is necessary for seedling establishment and sapling growth [15,17]. Variability in seed production [31] can also contribute to the locally high recruit numbers. The lack of support for the habitat hypothesis for recruitment agrees with results of a study by Kanagaraj *et al*. [32] showed that local recruit assemblages at the BCI plot showed weak and temporally inconsistent habitat structuring with respect to topographic variables. Additionally, several studies have attempted to explain the variability in local species richness of large plots of tropical forests by environmental variables [13,33–35] and found that an overwhelming proportion of this variability remained unexplained by local environmental variables. Seri *et al*. [11] also found very little effect of environmental heterogeneity on the dynamics of tree clusters over time.

Recent modelling studies outlined the importance of environmental stochasticity in explaining observed patterns in tropical forests (e.g. [5,12]). They found considerable temporal variability in tree population abundances beyond what is expected based on neutral theory. Our study, focusing on the spatial pattern of tree and species numbers at the local scale, reaches a similar conclusion, but with respect to spatio-temporal variability in local recruitment (and, to a lesser extent, mortality). Our study calls for a better understanding of the causes of the high spatio-temporal variability in local recruitment. For example, this variability can facilitate coexistence of competing species through a storage effect by producing resource partitioning through space or time [19,36]. Usinowicz *et al*. [19] suggested that the strong temporal variability and asynchrony in seed and seedling production observed at BCI is able to produce storage effects. However, it is not clear if our results can be interpreted as evidence for a storage effect because submodels R2 and R3 represent strong spatio-temporal variability in the total number of recruits without regard to species, but we did not explicitly include mechanisms that desynchronize recruitment of different species as required for the storage effect. Nevertheless, local clustering of species in combination with occasional ‘recruitment pulses’ at the quadrat scale may result in spatio-temporally desynchronized recruitment of species required for the storage effect.

The observed strong spatio-temporal variability in recruitment (with 21% of the quadrats having fewer than 10 recruits and 55% fewer than 19) generates for many species recruitment limitation (i.e. failure to have viable juveniles at all available sites [15,37]). Recruitment limitation can be a powerful force for maintaining diversity in species-rich communities, because it allows many sites to be won by ‘forfeit’ by species that are not the absolutely best competitor for the site [37,38].

### (b) Dynamics of local tree and species numbers

In the field data, we found surprisingly strong negative correlations between the change in the local number of trees and species over 25 years and the corresponding initial numbers (figure 3*a,b*). Quadrats with a high tree or species number tend to lose trees or species, respectively, while quadrats with few trees or species tend to gain trees or species. This finding also indicates that local clusters of high tree density and high species richness may move stochastically through the plot.

Our analytical model (equations (4.2) and (4.3)) for the dynamics of local tree numbers quantitatively explains the observed relationship between Δ*N* and *N*_{0}. It emerges if the number of recruits at quadrat scale is only weakly associated with (or independent on) the corresponding number of trees. In the case of a density-independent constant recruit input *R*, this model leads to a stable equilibrium at the local quadrat scale. A recruit input *R* independent of the local tree number can emerge in simulation model R1 + M1 as a consequence of a wide enough kernel that locates most recruits outside the quadrat of their parents. However, R1 + M1 severely underestimated the variability in local recruitment numbers. Instead, the good fit of models R2 + M1 and R2 + M2 (figure 2) suggested that the local dynamics at BCI is driven by spatio-temporally varying local conditions that allow occasionally at some sites establishment of a high number of recruits. This puts more emphasis on the spatial and temporal pattern of ‘safe sites’ that provide just the right conditions for recruitment, a perspective shared by many gap forest models [39]. Additionally, independence between the number of recruits in a quadrat and the actual number of trees will be fostered by long time delays between seed dispersal and actual recruitment (e.g. due to seed, seedling or sapling banks). Furthermore, models R2 + M1 and R2 + M2 approximate the observed correlations between Δ*N* and *N*_{0} and between Δ*S* and *S*_{0} well, although the observed correlations were stronger.

Previous work let us expect strong correlations between the local number of trees and species [40]. However, our best models (R2 + M1 and R2 + M2) with mean dispersal distances of 20–30 m produced too large correlations in the number of species and number of trees in the quadrats, and stronger species clustering due to shorter mean dispersal distance produced weaker correlations (electronic supplementary material, figure C6*a*). The data show that only one-third of the variation in local quadrat-scale species richness is explained by the number of individuals compared with two-thirds in our best models. This result suggests existence of additional small-scale spatial structures, beyond those created by the dispersal kernel, which make this relationship more stochastic. A candidate mechanism for this is seed dispersal by animals that may cause conspecific clumps of recruits (e.g. [41,42]).

### (c) Link with forest gap models

Ecosystems—even when they can be assumed to be in a state of overall equilibrium like an old-growth forest—are constantly changing, and different parts of these systems represent different possible states. This idea has a long history in ecology. It was first proposed by Aubreville [43] and made popular by Richards' ‘mosaic theory’ for forests [44]. This notion also underlies forest-gap models [39,45,46], which exhibit small-scale disturbance through gap-creation as a major driving force of forest dynamics, or individual-based Markov models of forest dynamics (e.g. [47]). Common to these approaches is the idea that different parts of a forest can be seen as partially linked and that it is the small-scale heterogeneity of these parts that drives the patterns and dynamics of forests as a whole. We used this approach to gain insight into its dynamics over time as well as the spatial distribution of local species richness and density. The actual processes governing the dynamics of tropical forests are highly stochastic, but when one uses a high number of replicates (e.g. 1250 plots) over sufficiently long time periods (25 years in our case), patterns in the fundamental characteristics of the forest emerge. Our results confirm our basic premise that the patterns of local species richness can provide additional information on the spatial and temporal dynamics of species rich communities.

## 7. Conclusion

Hubbell's neutral theory [3] has been highly successful in explaining important characteristics of tropical forests—like species abundance curves. Here we showed that focusing on local dynamics can provide important additional insights on demographic mechanisms that drive the spatial and temporal variation of local species richness in a neo-tropical forest. We detected strong spatio-temporal variation in the local number of recruits, much beyond that expected by demographic stochasticity alone. This variability, when combined with a constant *per capita* mortality rate, was critical for predicting the variation in local tree and species numbers.

Our research reveals that, despite the observed variability in mortality and recruitment, local tree numbers and local species richness seems to stabilize at typical values. However, more research is needed to investigate the causes. Investigating a high numbers of replicate plots bears high potential to improve our understanding of species-rich plant community dynamics.

## Data accessibility

The census data for the BCI plot are publicly available at the website of the Center for Tropical Forest Science: http://dx.doi.org/10.5479/data.bci.20130603. Use of the data has been agreed upon with the principal investigators of the BCI plot: Stephen Hubbell, Richard Condit and Robin Foster.

## Authors' contributions

M.K., A.H., F.M. and T.W. conceived of the study. M.K., F.M. and T.W. analysed the data. M.K. and F.M. implemented and analysed the null models. T.W. and M.K. wrote the draft of the manuscript including contributions of A.H. and F.M. All authors discussed the results and revised the paper.

## Competing interests

We have no competing interests.

## Funding

All authors acknowledge support from the ERC advanced grant 233066. A.H. was also supported by the Helmholtz-Alliance Remote Sensing and Earth System Dynamics.

## Acknowledgments

The authors thank Friedrich Bohn, Rico Fischer, David Alonso, José Capitán, Dominique Gravel, R. Condit, S.P. Hubbell, Fangliang He and an anonymous reviewer for comments on earlier versions of this manuscript. The BCI forest dynamics research project was founded by S. P. Hubbell and R. B. Foster, and is now managed by R. Condit, S. Lao and R. Perez under the Center for Tropical Forest Science and the Smithsonian Tropical Research in Panama. Numerous organizations have provided funding, principally the US National Science Foundation, and hundreds of field workers have contributed.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3868213.

- Received July 5, 2017.
- Accepted August 22, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.