Abstract

There is no consensus over which model is correct, because the degree to which different processes can be discerned from SAD patterns has not yet been rigorously quantified. We present a power calculation to quantify our ability to detect deviations from neutrality using species abundance data. We study non-neutral stochastic community models, and show that the presence of non-neutral processes is detectable if sample size is large enough and/or the amplitude of the effect is strong enough. Our framework can be used for any candidate community model that can be simulated on a computer, and determines both the sampling effort required to distinguish between alternative processes, and a range for the strength of non-neutral processes in communities whose patterns are statistically consistent with neutral theory. We find that even data sets of the scale of the 50 Ha forest plot on Barro Colorado Island, Panama, are unlikely to be large enough to detect deviations from neutrality caused by competitive interactions alone, though the presence of multiple non-neutral processes with contrasting effects on abundance distributions may be detectable.

Author Summary

A classic idea in Ecology is that species coexist because they occupy different “niches”. However, random processes such as dispersal could also explain species coocurrence, without invoking niche differentiation. “Neutral” models embody this idea, omitting niche differentiation and assuming all species are identical. Such models are mostly statistically consistent with the relative abundances of tree species in tropical forests, but statistical procedures always contain an element of uncertainty and many other models could also be consistent with a particular data set. We compute how strong the non-neutral processes would need to be in order for their effect to be detectable in data sets of different sizes. We find that the largest ecological data sets currently available, such as the 50 hectare plot on Barro Colorado Island in Panama, are not large enough to distinguish between neutral and non-neutral models, unless multiple non-neutral processes are at work. This means that other types of pattern need to be studied, or larger data sets collected, in order to understand the mechanisms behind forest biodiversity.

Introduction

The appearance of common patterns in species abundance distributions (SADs) for different communities suggests that the same ecological mechanisms structure these communities [5, 6]. However, it is now thought that many patterns describing communities are rather insensitive to these processes [7—11]. For example, empirical SADs are in many cases found to be statistically consistent with Hubbell’s neutral theory [12—17], but this does not mean that communities are truly neutral because non-neutral models can predict similar [8—10, 18, 19] or even identical [4] patterns. This raises the question of whether anything can be inferred from fitting it to SAD data [3, 4, 20].

Neutral theory has many virtues [22—24] and in many ways it is more complete in scope than competing niche theories [25]. It describes community dynamics at the individual level, treating births, deaths, and dispersal as stochastic processes. It is susceptible to rigorous statistical tests, because unlike many other demographic models the likelihood of obtaining a particular community or sample can be computed exactly [26—30]. Even the controversial neutral assumption that interactions between individuals do not depend on species identity is inspired by biological reality; Hubbell observed that, in tropical forests, all species compete for light—and, therefore, space [31]. This means that neutral theory should be a good starting approximation for communities of sessile species that compete for a common resource, such as space (e.g. tropical trees or coral reefs). More realistic models will include non-neutral processes, such as interactions that depend on species identity [32, 33], but neutral theory can act as a null model for assessing the weight of evidence for such processes.

While there may be patterns or scales for which some processes are undetectable, e.g. due to central limit theorem-like effect [19, 35, 36] , strong interactions between individuals can structure communities and it is in some cases possible to detect their existence from inspection of the SADs [3, 33, 34]. If a data set is found to be consistent with neutral theory, we should therefore be able to infer that some particular non-neutral processes are not present in that community, or at least are not strong enough to produce detectable deviations from neutrality in a data set of this size.

Our purpose is to estimate an upper bound for the strength of non-neutral processes in tropical forest data sets [37—39] that have been found to be consistent with neutral theory [13]. To do this, we fit the standard neutral model (SNM) to data sets generated by a non-neutral model, and compute the probability of rejecting neutral theory. We test the neutral null hypothesis using a maximum likelihood approach (using an exact expression [26] for the likelihood of a sample from the SNM), where p-values are evaluated by a parametric bootstrap procedure.

The power therefore depends on which alternative hypothesis is true. In this paper, we focus on two classes of non-neutral processes: interspecific competition, and intrinsic (density independent) fitness differences between species. Interspecific competition is one of the classic mechanisms that promote coexistence [33, 40, 41], whereas differences in fitness represent the fact that the mean environmental conditions in a particular area of habitat will tend to favour one species over another. These represent opposite ends of a spectrum of possible non-neutral models, because symmetric interspecif1c competition tends to lead to equal abundances among species, whereas intrinsic fitness differences tend to lead to highly uneven abundances. While these are only two examples out of an infinite set of non-neutral models, our method provides a blueprint for computing the detectability of any type of non-neutral process.

Our models are similar in structure to Hubbell’s standard neutral model (SNM), in that we consider stochastic population dynamics in a local community where strong density dependence regulates the total community size to I individuals, coupled by immigration to a much larger metacommunity. We consider two models of interspecific competition: one, which we shall denote HL, is a multi-species stochastic Lotka-Volterra model similar to that studied by Haegemann and Loreau [33]; the other, denoted by PC, has stochastic Ricker-like dynamics as studied by Pigolotti and Cencini [42]. Our model of intrinsic differences in fitness, denoted by IF, assumes that the fecundity of each species is a randomly generated variable. Each of our models has a single parameter that determines how strong the non-neutral processes are. In model HL, parameter 7/ represents the relative difference in strength between interspecif1c and intraspecif1c interactions, so that when 7/ = 0 the dynamics are neutral whereas when 7/ > 0 coexistence is promoted. In model PC, parameter c determines the difference between inter and intra-specif1c density dependence, so c = 0 corresponds to neutral interactions and non-neutrality becomes stronger as c increases. In model IF, the fitness of each species is generated from a Gamma distribution with shape factor l/k, so that when k = 0 all f1tnesses of the species are the same and the local dynamics are neutral.

The proportion of immigrants of different species follow their relative abundance in the metacommunity. We consider two cases: in case LOGS the metacommunity is described by a logseries with fundamental diversity constant 6, and in case EVEN the metacommunity has ST species which all have equal abundance. A logseries distribution can arise from many processes, including but not restricted to neutral dynamics [43]. We considered the even metacommunity because it represents a meta-community limit of our local community dynamics, and as a result represents a contrasting, extremely non-neutral, limit to the logseries. When coupled to the LOGS metacommunity, each of our models should be equivalent to the SNM when the local dynamics are neutral (when 7/, c, or k equals zero). As the dynamics are made more non-neutral, the deviations from the SNM should become stronger, and we expect the power of the test of the neutral null model to increase. However, when coupled to the EVEN metacommunity, the models are not equivalent to the SNM even when the local dynamics are neutral. In this case the power of the test could be high even if the local dynamics were neutral, though if] is very small the statistical power could still be low.

First, we explore how the parameters of our models affect the probability of detecting non-neutral processes. For a real community we do not know a priori the appropriate parameters to use, so we need to choose the parameters so that the alternative model gives comparable patterns to the empirical data. In the second part of our study, we estimate the power of tests of neutrality for empirical data from three New and Old World tropical forests, including Barro Colorado Island (BCI) in Panama. Our power calculation provides an estimate of the smallest sample size that is needed to detect non-neutrality of known intensity, and of the range strengths of non-neutrality needed to reject neutrality for a given species abundance data set.

Results

Power calculation for fixed non-neutral model parameters

When interspeciflc interactions are non-neutral (models HL and PC, top and middle row of Fig. 1), we see a simple pattern: as the strength of non-neutral processes is increased (7/ or c increases from zero), the power of the test increases. In addition, for these models the power of the test increases as the local community size I is increased. This shows that any strength of non-neutrality, however weak, can in principle be detected provided the data set is large enough. However, the system sizes needed to give a significant power may be too large to be empirically accessible when local dynamics are nearly neutral (7/ or c close to zero). For the LOGS metacommunity, and when the local dynamics are strictly neutral (7/ = 0 for model HL or c = 0 for model PC), the models are equivalent to the SNM, and the power is equal to the threshold p-value for statistical significance (0.05 in our study). However, for the EVEN metacommunity the power can be higher than this threshold even when the local dynamics are neutral, because the immigration process makes the model no longer equivalent to the SNM (though this is not visible for the relatively small values of m used with models HL and PC in Fig. 1)

The power is again low when the local dynamics are neutral (k —> 0) and the metacommunity follows a logseries (bottom left panel). However, the power does not increase monotonically as the non-neutrality parameter k is increased. This is because strong selection rapidly leads to dominance by a single species [44], especially in small communities, which is a pattern that can also arise from the SNM if the immigration parameter m —> 0. Moreover, the power no longer increases monotonically when the local community size I increases; for example the power for I = 2000 in Fig. 1, bottom left, is higher than for I = 200, 5000, or 20000. This appears counterintuitive because statistical power should increase monotonically with sample size. However, I represents more than just the amount of data available: it is a parameter which interacts non-linearly with the model dynamics. In the IF model, for instance, it determines whether the dynamics are in the strong or weak selection limit, and it also plays a nonlinear role in the SNM.

In this special case, there is a single dominant species, relative to which all other species have zero fitness. All local recruits will therefore be of the dominant species, though other species will also be present due to immigration. This case is particularly simple because the species identity of each individual in the local community is the dominant species with probability 1—m(1—1/ST), and each of the other species with probability m/ST. We find that the power of the test of neutrality is low at small I, increases to a maximum at an intermediate value of I, and then decreases as I increases again (Fig. 2). Because, for the non-neutral model in this limit, I is nothing more than a sample size—a community of size 2 I can be constructed by adding two communities of size I—this non-monotonic relationship between power and I must be due to the nonlinear role played by I in the community dynamics in the neutral null model. The two other model parameters (immigration rate and diversity of the metacommunity) also strongly affect the chance of rejecting the neutral hypothesis. The parameter m describes Local community size (J) the probability that a newly born individual in the local community is an immigrant from the metacommunity, so the local community resembles the metacommunity more closely as m is increased; when m = 1, the local community is effectively a random sample from the metacom-munity and local dynamics are irrelevant. Increasing 6 in model LOGS, and increasing ST in model EVEN, lead to more diverse metacommunities, and as a result tend to increase diversity in the local community.

It is clear why increasing m should reduce the power for models HL/LOGS and PC/LOGS, because the local community is then more like a logseries, and hence more like the SNM. The power changes very little between m = 10—4 and m = 10—3, reflecting the much greater importance of local dynamics on the patterns when m is small. It is less clear why increasing 6 should reduce the power of the test, though it is worth noting that both increasing 6 and increasing m have the effect of increasing the local richness. Increasing m or ST in models HL/EVEN and PC/EVEN also increases the local richness, and while it is not obvious why this should make the model resemble the SNM, we find that it also reduces the power of the test.

3 and 4). For model IF/LOGS when k x 10—3, for Table 1. Fitted parameters for CTFS data sets and probability of rejecting the neutral null hypothesis, for model PCILOGS.

On the whole, however, both m and the community diversity tend to increase the power of the test, which is the opposite effect from what is seen in models HL and PC. This is because the local dynamics tend to lead to monodominant states, which as explained before are indistinguishable from the SNM even if their origin is highly non-neutral. Processes which increase the local diversity allow the non-neutral features of this model to be better detectable.

This would be a very useful relationship if it were found to hold in general, because power calculations for large I are very expensive compu-tationally and this would enable us to estimate power by varying c as a proxy for I. We find, however, that this behaviour is not preserved for other parameter values. We did not find any simple way to summarise the dependence on model parameters evident in Figs. 1, 3, and 4 that would enable us to estimate the power outside of the parameter range we tested explicitly.

Power calculation for large forest surveys

SNM has been found to be statistically consistent with several large tropical forest data sets [12—17], but this does not mean that SNM is an exact description so a power calculation gives us, in principle, an upper bound for the degree of non-neutral processes in these systems. We do not know a priori the appropriate non-neutral parameter values for these forests, but we can choose model parameters so that the model data match a number of features of the empirical data. Specifically, we chose model parameters so that the community size, mean species richness, Table 3. Fitted parameters for CTFS data sets and probability of rejecting the neutral null hypothesis, for model IFILOGS.

No solutions were found for k 2 0.01. Pasoh is not included because no parameter sets consistent with the empirical data were found with k 2 10—4. and mean Shannon diversity of sample data sets from the non-neutral model are the same as in the empirical data set. More details are given in the Methods section.

The determination of appropriate model parameters, and the power calculation itself, is very computationally expensive, so we performed this for three candidate strengths of non-neutrali-ty, and for models PC and IF only. The results are summarised in tables 1—4.

This is because the non-neutral process tends to increase the evenness in the community, so when c is increased the fitted immigration rate needs to be increased in order to match the Shannon indeX in the empirical data. In other words, the parameters of this model need to be close to neutral (either c small, or m close to 1) in order to agree with the empirical data, so the power of the test is low.

This is because both the metacommunity and the local dynamics are non-neutral, but the local dynamics tends to lead to very uneven (i.e. monodomi-nant) communities while the metacommunity tends to increase the evenness of the community. These processes need to be in balance in order for the model to match the diversity and richness of real data, and as a result the fitted model is far from neutral. Table 4. Fitted parameters for CTFS data sets and probability of rejecting the neutral null hypothesis, for model IFIEVEN.

The largest value of c for which the model is consistent with the empirical forest data was found to be when ST —> oo, in which limit the immigration process behaves effectively as a speciation process. It was found that when the model was fitted to the Pasoh data with c = 0.7523, the probability of rejecting the neutral null hypothesis was 0.20. The power of the test was much lower when c = 0.1, and was always low when the model was fitted to the Lambir data. The model could not reproduce the richness and Shannon diversity of the BCI data set unless c < 0.1.

Here, both the non-neutral local dynamics and the metacommunity tended to lead to highly uneven abundance distributions, and as a result the Shannon index produced by the model was very low unless k was small. It was not possible to fit the model to Pasoh when k 2 10—4, or to BCI or Lambir when k 2 0.01. While there were two discrete parameter sets each that matched the richness and evenness of Pasoh and Lambir (one where 6 was low and m high, and one where 6 was higher and m lower), the power of the test with these parameters was always low.

Discussion

This contradicts the suggestion that the SAD for large samples will approach the same canonical form, and that larger sampling efforts would consequently be futile [35]. Indeed, the SNM has been rejected using SAD data from very large phytoplankton communities [45], and we found that the SNM could also be rejected for the tropical tree species abundances of Yasuni National Park (see Methods). Our results also show that independent niches can be distinguishable from neutrality, contrary to suggestions by Chisholm and Pacala [36], because the test is most powerful when the species in our model HL undergo independent stochastic logistic dynamics (i.e. when 7/ approaches 1, see Fig. 1). Our HL model with y = 1 differs from the independent-niche models of Chisholm and Pacala [36] and Haegeman and Etienne [19] by having strong intra-specif1c density dependence, so the marginal distributions are very different from the SNM.

As shown in Tables 1 and 2, statistical power remains very low as c (the parameter measuring the intensity of inter versus intra-specif1c competition) is varied. This is because the model parameters required to give the same richness and evenness in the empirical data are themselves close to neutral, either because c is small or because the metacommunity follows a logseries distribution and m is close to 1 (in which case the local community strongly resembles a neutral-like metacommunity). This result is in agreement with the good fits of some species-indepen-dent neutral models to a large number of SADs for very diverse communities [46]. It is interesting to note that Volkov et al. [47] estimated that interspecif1c species interactions are many times smaller than intraspecif1c interactions in tropical forests, but we cannot apply their results to our models because they did not include immigration from a metacommnunity.

In model IF/EVEN, intrinsic local fitness differences tend to decrease richness and evenness, whereas the even local community increases richness and evenness. This means that both processes can be strong while still producing levels of richness and abundance consistent with the empirical data. Du et al. [10] noted that non-neutral processes which have opposing effects on relative abundance distributions can lead to abundance distributions that resemble neutral theory, but our investigation shows that they can still be distinguished from SNM in some cases. We can therefore conclude that such a combination of strongly non-neutral processes is not present in data sets for which the SNM is not rejected, such as the three CTFS forests we studied in this paper.

In our analysis, effect size is encoded in the parameter values of the non-neutral model under consideration (respectively, 7/, c, and k for models HL, PC, and IF). When density dependence is non-neutral (models HL and PC), power increases as interactions become more non-neutral (see Fig. 1). However, for non-neutral intrinsic fitness (model IF) the power of the test depends non-monotonically on k (Fig. 1, bottom row). These patterns can be understood from the effect that these parameters have on diversity patterns—strong non-neutrality (k large) in model IF leads to monodominance, which is indistinguishable from the neutral model with strong dispersal limitation (m very small).

2). In most standard statistical tests, a “sample” consists of a number of statistically independent measurements. In an ecological community (or a model thereof), the individuals are not statistically independent because of their interactions (whether within or between species). This is true even in the SNM: an equilibrium community of size I can be generated as a hypergeometric subsample of a community with larger I [48], but the individuals are not independent because this represents sampling without replacement. This means that the community size I plays a nonlinear role and is not a straight analogue of the sample size in standard statistical tests, so statistical power does not necessarily increase monotonically with I.

Increased metacommunity diversity decreases statistical power for models HL and PC (Fig. 4), which echoes the observation that higher local diversity leads to SADs that look more like those created by the SNM even in the presence of niche structuring [19, 36]. This suggests that it might be easier to quantify non-neutral interactions in less diverse forests [19]. However, this is not true for all types of non-neutral processes: for model IF the power increases when the metacommunity diversity is increased.

The power calculations in this paper are very computationally expensive, and it would be unfeasible for us to repeat them for I much larger than ~ 30000 individuals. Moreover, to do this we would need to know how parameters y (for model HL) and c (for model PC) are affected by I. Our notation tacitly assumes that each species is sensitive to mean population densities over the whole community, but in real systems, where individuals of a species are clumped together, an individual will only interact with nearby individuals so the values of y and c might depend on I. We are therefore unable to estimate the factor by which CTFS data sets would have to be enlarged for us to distinguish model PC from the SNM. In this paper, we have analysed a range of non-neutral scenarios: non-neutral density dependence affecting mortality or recruitment; non-neutral differences in intrinsic fitness; neu-tral-like or extremely non-neutral metacommunity. These processes have contrasting effects on the SAD, so arguably represent the extremes of a spectrum of possibilities. Nevertheless, there are many types of community processes which are not encompassed by our models. For example, trophic or mutualistic interactions are not present in our models and should lead to very different patterns of abundance, though these are more likely to be relevant in other systems than the tropical forests on which we focus. Similarly, dynamics that lead to multimodal SADs should be relatively easy to distinguish from neutrality [34, 49]. We have made a number of simplifying assumptions to keep the number of parameters manageable, but our framework could still be used to perform a power calculation for any type of non-neutrality that can be incorporated in a simulation model. It would also be preferable to perform power calculations for spatially explicit models, which represent a more realistic dispersal process and can readily be simulated [50]; however, for our test we would need a likelihood for the spatially explicit neutral model, which is not currently available.

[53]) that SAD data do not have sufficient resolving power to assess the importance of non-neutral processes in structuring forest communities. Our study shows that even large-scale tropical forest data sets are not large enough, or are too diverse, to detect non-neutral species interactions using the SAD alone. However, we would not expect to see a good fit to the SNM if these forests contained multiple processes with opposing effects on richness and evenness. Patterns that include more information, such as multiple samples [27], spatio-temporal changes [54, 55] or phylogenetic data [56, 57] , are likely to be much more revealing about the processes that generated them. Provided it is possible to compute the likelihood for obtaining such patterns in a neutral model, our approach can be adopted to calculate the sampling effort needed to detect and quantify non-neutral processes, and understand the forces that structure communities.

Methods

This section describes (i) our non-neutral alternative models, and methods for generating samples from them; (ii) our method for testing whether to reject the neutral null hypothesis for a particular data set; (iii) the method for combining (i) and (ii) to give a power calculation; (iv) the method for estimating model parameter values in order to estimate the statistical power of particular experiments.

As with the SNM, strong local density dependence is assumed to keep the local community size fixed at I individuals. A fraction m of all recruits immigrate from a “meta-community”, which is assumed to be large enough for the relative abundances (P,- for the i’th species) to be effectively static in time. This immigration prevents drift to local monodomi-nance. The models differ in the relationship between the local abundances and the birth and death rates, in the relative abundances in the metacommunity, and in whether the dynamics are syncronous or sequential in time. Local model HL

It can be thought of as a multispecies stochastic Lotka-Volterra competition model with immigration, where a single parameter controls the relative strength of inter-and intra-specific interactions. Our model differs from that of Haegeman and Loreau [33] in that each death event is immediately followed by a single birth event so that the local community size remains constant. We consider a local community consisting of] individuals, each of which has a species identity which is an integer between 1 and ST. Mortality is affected by inter and intra-specific density dependent mortality, so that the probability that the next individual that dies has species iden-tiy i is proportional to Where 11,- is the number of individuals of species 1' in the community, so that I = 2:1 11].. In the neutral case, 7/ = 0, the mortality rate is the same for all individuals irrespective of species, but When 7/ > 0 per capita mortality is greater for more abundant species. The dynamics proceed by choosing the individual that dies next, so that the probability that the dead individual has species identity 1' is proportional to M ,-. A recruit is then chosen to be of species 1' with probability where P1- is the relative frequency of species 1' in the metacommunity. One time step consists of I of these elemental update steps. If y = 0, the mortality rate is independent of species identity, so species interactions are neutral; when 0 < y g 1, niche differentiation tends to promote species coexistence [33]. The model is ill-defined when 7/ > 1, since that would lead to M1- < 0. At first sight, it might appear that a more general model could be obtained by using the functional forms in Haegeman and Loreau [33], which allow for density independent as well as density dependent mortality. The rates in Equation (13) of that paper correspond in our notation to where both M ,- and F,- are now rates rather than probabilities (no longer normalised so that they sum to unity) and we have introduced the factor ST P1- to allow the immigration rate to differ between species. Here, r+ and r_ denote respectively the rates of density-independent birth and mortality, Kl plays the role of a carrying capacity, and a tunes how neutral the interactions are are (neutral when 1; non-neutral competition for 0 < a < 1; mutualistic for a < 0). A little algebra shows that Equations (3) and (4) are equivalent to Equations (1) and (2) (up to overall prefactors that do not affect the sequence of processes in the simulation) with the choice of parameters The values of these parameters is Within the range for which our model is defined (0 < m g 1, defined by Equations (1) and (2) captures the apparently more general density dependence de-order for there to be a nontrivial equilibrium) and — fined in Equations (3) and (4), except for the case of strongly mutualistic interactions.

Local model PC

As in model HL, we assume that interspecific interactions are weaker among het-erospecifics than among conspecifics, but in model PC we assume a Ricker-like functional form that acts on fecundity so that the number of local propagules of species 1' E {1, . . ., ST} is or 11,- exp(—an1- — 192,- 75 J- nj), Where nk is the local abundance of species k and a and b(< a) are constants. The fraction of local propagules that are of species 1' is then where c = ](a — b). If the interactions are sensitive to the average local density of the different species, i.e. all species are spread throughout the community, then for the same pool of species we expect c to be independent of]. Spatial effects could lead to a focal species only being sensitive to the dynamics of nearby species, in which case the effective value of c would depend on I, although this can only be modelled correctly using a spatially explicit model. Because a fraction m of recruits are immigrants from the metacommunity, the probability that a new recruit is of species 1' is In this model, we assume that the generations are discrete and non-overlapping: at each time-step, we compute the Rl- from the current configuration, and generate a new configuration using

Sequential updating (as in model HL above) is a more faithful biological description of triopical forest dynamics, but it is known that in the neutral limit the sequential model (Moran process) and the synchronous model (Wright-Fisher process) give indistinguishable equilibrium statistics for the large community sizes we are interested in. The syncronous model is therefore well suited to our goal of exploring the detectability of departures from neutrality. The only circumstances where the synchronous model behaves qualitatively differently from the sequential one is where the Ricker-like dynamics tend to lead to limit cycles or chaos, but that does not affect the results in this paper because we always choose c < ZST (see Metacommunity model EVEN below). Local model IF

By contrast, model IF could be interpreted as considering environmental variability between communities. In any given local community, each (out of a total number ST) species has a different intrinsic fitnessfi, different from the other species and different from its intrinsic fitness in other commuities. Our model introduces intrinsic fitness differences in a similar way to Chesson and Warner [58] , though our model is otherwise different because our fitness differences do not fluctuate in time and our dynamics are stochastic rather than deterministic. For each realisation, we generate the from a Gamma distribution with mean shape factor l/k, so that k“2 is the coefficient of variation among the f1tnesses, and all the f1tnesses are equal when k —> 0. The fraction of local propagules that are of species 1' is then and the probability that a recruit is of species 1' is As was the case for model PC, for computational efficiency we assume discrete generations, and simulate the model using a multinomial pseudorandom number generator with probability vector {Ri}.

Metacommunity model LOGS

Following the SNM, in model LOGS we assume that this abundance distribution follows a Fisher logseries with diversity parameter 6. This is often a good description of empirical data, and can arise from several models of community dynamics including Hubbell’s neutral model [31, 59]. Note that the metacommunity represents the pool from which immigrants can be drawn, which could comprise many disparate communities. Therefore, the metacommunity does not necessarily correspond to the large-I limit of a single local community model. This means that one reasonable scenario is that the metacommunity follows a canonical form form due to averaging over very large scales [35], even when the local dynamics are non-neutral. In metacommunity model LOGS we use the distribution introduced by Ewens [60] to give the number fM of species in the metacommunity with relative abundance x within the interval (x; x + dx).

In practice, we use this distribution to simulate from a very large metacommunity containing ST species; our full sampling algorithm is described in 81 Text. The results in this paper are for ST 2 2000, which is large enough to be effectively infinite (i.e. a choosing a larger ST did not have a perceptible effect on community statistics or the power of the test of neutrality, but did increase the duration of the simulation).

Metacommunity model EVEN

As explained above, the logseries is a reasonable candidate metacommuity model even when local community dynamics are non-neutral. However, it is also the metacommunity model in the SNM, so we want to consider the possibility that non-neutral processes are visible in the metacommunity as well. In our EVEN metacommunity model, there are ST species with equal L ST neutral and non-neutral community dynamics [33, 54, 61, 62] , though it has little empirical support. It is appropriate to use this distribution in our study because, as we shall show, it represents a metacommunity limit of our non-neutral local community models. The HL model is a form of stochastic multi-species Lotka-Volterra model, so its large-I limit is described by differential equations. When 0 g a < 1, this has a stable equilibrium with 2 relative abundance, P1. = . This distribution has been used in previous modeling studies of

When I is large, the multinomial distribution becomes sharply peaked around its mean value, so from Equation (5) (for vanishing m) the dynamics of r1. 2 follows A standard stability analysis shows that the equilibrium r1. = i is stable provided CST < 2; for higher values of c the community displays limit cycles or chaos.

However, the metacommunity represents an aggregate of many independent local communities, and we expect different species to dominate in different communities. The model assigns fitnesses independently at random to the different species, so we expect each species to have the same overall relative abundance i in the metacommunity.

Other metacommunity models could be obtained by taking the limit in different ways. For instance, if parameter c in model PC depends on I and approaches zero sufficiently rapidly in the limit I —> 00, then the metacommunity limit would be neutral and follow the same Ewens distribution as model LOGS. If the fitnesses in model IF were not i.i.d. random variables, but rather different species had different mean f1tnesses, then the metacom-munity would have a different, uneven distribution. While there is an infinite variety of possible metacommunity distributions, the EVEN metacommunity represents the most contrasting distribution to the logseries, in the sense that it has the maximum Shannon diversity index for a given species richness while the logseries is a very uneven distribution. It also has the advantage of being characterised by a single parameter (ST), whereas other commonly-used distributions (e.g. the lognormal) generally require two parameters.

Testing the neutral null model

In order to quantify Whether a particular data set is consistent with neutral theory, we adopt a maximum likelihood approach together with a parametric bootstrap as used by Walker and Cyr [45] and Rosindell and Etienne [63]. To calculate the p-value of our test, we compare the value of a test statistic for the test data set with values of the test statistic for data sets generated by the null model. We choose the maximized likelihood of the neutral model as our test statistic. The likelihood L(X|m, 9) that the neutral model would generate a data set X, for parameters (m, 9) is computed using the exact formula derived by Etienne [26]. We use code based on Tetame (http://chave.ups-tlse.fr/projects/tetamehtm), an efficient implementation in C++ of Etienne’s formula that was developed by Iabot et al. [64]. We have ported the Tetame code to C, and adapted it so it can be loaded as a dynamic library in R using the .C() function. The hypothesis test consists of the following steps: 1. For a test data set XT, we find the maximum likelihood parameter estimates (mMT, GMT , i.e. the set of parameters for which L(XT| m, 6) takes its largest value, L(X Tl m MT, GMT), the maximum likelihood estimate. 3. For the i’th sample neutral data set, compute the corresponding maximum likelihood estimate parameter set ijN) and maximum likelihood L(X | mm, 03%,) using the same procedure as used for XT.

The p-value for the test is the fraction of neutral data sets Whose maximum likelihood is lower than the maximum likelihood for the test data set, i.e. p is estimated by: Where f is the step function (f is equal to 1 if its argument is positive and 0 otherwise). 5. The neutral model is rejected if the p-value is less than the chosen threshold for statistical significance, which we take to be 0.05.

Statistical power calculation for fixed non-neutral model parameters

The statistical power can only be quantified by specifying an explicit model to represent the alternative hypothesis. The power can be computed by simulating many data sets from the alternative model, and performing a test of the null hypothesis on each data set, as explained above. The power is the fraction of cases for which the null hypothesis is rejected. A Type 11 error is defined as the failure to reject the null hypothesis when the alternative hypothesis is true, so the power is equal to 1 — fl, where fl is the probability of Type II errors.

When coupled to the LOGS metacommunity, our models are equivalent to the SNM in the limit where the non-neutral parameter (7/ for model HL; c for model PC; k for model IF) is zero, so 7/, c, or k is our effect size for these models. When coupled to the EVEN metacommunity, our models are never strictly equivalent to the SNM so the effect size cannot be defined. In general, the power of the test will also depend on other model parameters, so we need to perform power calculations for a wide range of potentially interesting parameter values. To calculate the power of tests of neutrality for a non-neutral model with a particular set of model parameters YT, we use the following procedure: 1. Generate a large number (in this study, at least 100 and usually more than 400) of equilibrium data sets from the non-neutral model With parameters YT, by simulating the model for a sufficient number of time steps. The number of time steps was chosen to be at least ten times the number of timesteps such that the species richness and Shannon diversity index in the local community appeared to have reached their equilibrium values; this number depends on the model parameters (e.g. fewer time steps are needed When m is close to 1); 2. For each data set, perform a test of the neutral null model using the parametric bootstrap method described above.

The power of the test is the fraction of non-neutral data sets for which the test was significant (i.e. the neutral null model was rejected). Because each non-neutral and neutral data set is statistically independent, the power is a binomial proportional random variable. Where shown, confidence intervals are 95% Ieffreys intervals [65].

Power calculations for particular experiments

The power of the test is a property of the ensemble of data sets that could be produced if the data were generated by a non-neutral alternative model. This ensemble of data sets depends on the model parameters which are chosen, but for a particular data set we do not necessarily know the appropriate parameters to use in the alternative model. Here, we choose parameter sets such that the model best describes a set of summary statistics of empirical data sets, specifically: species richness and Shannon index. Once the strength of the non-neutral process (or, c, or k as appropriate) and the local community size I are chosen, the model is char-acterised by two further parameters: the immigration rate m and the diversity (6 for model LOGS or ST for model EVEN) of the metacommunity. There will therefore be a discrete set of parameter values where the mean species richness and mean Shannon index of samples from the model match the empirical data sets; in most cases we found only one such parameter set, though in some cases there were two and in others there were none because the model produced a Shannon index that was always higher than, or always lower than, the empirical data. We performed this procedure, for a set of candidate values of the non-neutral parameter, to generate parameter sets resembling three tropical forest data sets belonging to the CTFS network to which neutral theory has successfully been fitted in the past [46]: Barro Colorado Island, Pasoh Forest Reserve, and Lambir Hills National Park [37—39]. For each forest, and separately for each survey year, we tested the null model that the data were generated by the SNM using the parametric bootstrap method described above. In each case we found p > 0.05, showing that the data were statistically consistent with SNM. The mean total community size, species richness, and Shannon index for these sites, averaged over the census years available at http://www.ctfs.si.edu, are given in Table 5. These sites were selected because they had higher Shannon index than a logseries distribution with the same size and richness, so we expected that a model with non-neutral interspeciflc interactions or an EVEN metacommunity would describe the data better than the SNM. Table 5. Summary statistics for the three tropical forest data sets to which our non-neutral models were fitted.

[46] have compared the SNM to three other CTFS forest sites, but the Shannon index in these sites is lower than (Korup National Park and Yasuni National Park) or almost equal to (Sinharaj a World Heritage Site) that of a logseries, so we expected it to be more difficult to find suitable model parameters. We also found that the Yasuni National Park data were not consistent with SNM (p = 0.001 for 1996 and p = 0.014 for 2003), though we found p > 0.05 for all Korup and Sinharaja surveys.

Supporting Information

Sampling from an infinite metacommunity. (PDF)

Acknowledgments

The BCI forest dynamics research project was made possible by National Science Foundation grants to Stephen P. Hubbell: DEB-0640386, DEB-0425651, DEB-0346488, DEB-0129874, DEB-00753102, DEB-9909347, DEB-9615226, DEB-9615226, DEB-9405933, DEB-9221033, DEB-9100058, DEB-8906869, DEB-8605042, DEB-8206992, DEB-7922197, support from the Center for Tropical Forest Science, the Smithsonian Tropical Research Institute, the John D. and Catherine T. MacArthur Foundation, the Mellon Foundation, the Small World Institute Fund, and numerous private individuals, and through the hard work of over 100 people from 10 countries over the past two decades. The plot project is part the Center for Tropical Forest Science, a global network of large-scale demographic tree plots.

Author Contributions

Performed the experiments: OAH DA SIC. Analyzed the data: OAH SIC DA. Contributed reagents/materials/analysis tools: DA RSE. Wrote the paper: OAH DA RSE SIC.