Beyond neutrality: adding habitat filtering to neutral models

Abstract

Background

Understanding the processes that structure species is one of the primary focuses in community ecology. Hubbell’s neutral model shows stochastic processes alone can describe the two macro-ecological patterns, species richness and species-area relationship, of the community. Although Hubbell’s neutral model can explain the macro-ecological patterns of the species at large scales, it paid less attention to construct the spatial structure of the community. Previous studies suggest that such spatial structures are mostly due to habitat filtering processes work at the intermediate spatial scales. Therefore, Hubbell’s neutral model does not explain the full picture of the community structuring due to its fully stochastic nature.

Results

In this study, we proposed a two-schema model that has the habitat filtering component and the stochastic component to construct the species assemblages seen in the community level. The proposed model uses one additional parameter (i.e. number of individuals in habitat) in addition to Hubbell’s three-parameter neutral model (i.e. fundamental bio-diversity number (θ), dispersal limitation (m) and speciation (v)). The proposed model works at two spatial scales: habitat filtering at the intermediate scales and stochastic processes at the large and very small spatial scales. The model coupled the local community dynamics with the meta-community dynamics. The local community has a fixed area with carrying capacity that is proportional to the local community size. The number of habitats in the proposed model can vary. Individuals are placed into habitats with probabilities according to the habitat suitability. Species richness and species composition in each habitat were calculated. The model is fitted for different θ values, m values, and a different number of habitats.

Conclusions

We assume that habitat filtering plays an important role together with stochastic processes to structure species in forests. Therefore, the proposed model with only four parameters can explain a large proportion of the species structuring of the communities. We found that more species can be maintained in a heterogeneous environment than a uniform environment. Therefore, habitat conservation is highly important for maintaining species diversity in forest communities.

Background

All species compete for a common set of resources, such as water, nutrient and space. The competitive exclusion principle emphasises, in long run, two-species cannot coexist if they use a common set of resources (Gauss, 1934). Hence, the competitive exclusion principle states that, in the long run, one species becomes dominant by winning the full space while all the other species go extinct (Gauss, 1934). However, the coexistence of thousands in of species in tropical forests challenges this theory. Although majority of species, in tropical forests, are rare or singletons, their survival increases by habitat filtering, negative density dependence and dispersal syndromes (Janzen 1970; Connell 1971; Harms et al. 2004; Gunatilleke et al. 2006; Seidler and Plotkin Chuyong et al. 2011; Liu et al. 2018). Habitat filtering hypothesis states that habitat suitability affects birth rates, death rates and speciation rates of species. This theory also assumes such species-habitat associations maintain species richness by delaying the species mono-dominance.

In contrast, Hubbell’s neutral theory, a fully stochastic process, proposed species are ecologically equivalent. Therefore, per capita birth rates, death rates and speciation rates are the same for all the species (Hubbell 1997; Hubbell 2001). The neutral theory suggests species richness in a local community can be maintained by the stochastic birth-death process at the local community and immigration of species from the meta-community. In the meta-community, speciation introduces new species (Hubbell 1997; Hubbell 2001). This theory successfully predicts two macro-ecological patterns (i.e. species abundance distribution and species-area relationship) using three parameters (i.e. the fundamental bio-diversity number, the immigration rate or the fundamental dispersal number and the speciation rates) (Hubbell 1997; Hubbell 2001; McGill 2003). This theory has great popularity among ecologist due to its simplicity (Rosindell et al. 2011, Rosindell et al. 2012). However, Hubbell’s neutral model that placed individuals randomly does not incorporate the spatial arrangement of the species. Spatial patterns of the species hide large information that is highly important to understand the underlying processes on species structuring at different spatial scales (Hubbell 1979; He et al. 1996; He et al. 1997; Condit et al. 2000; Wiegand and Moloney 2004; Gunatilleke et al. 2006; Wiegand et al. 2007; Getzin et al. 2008; Wang et al. 2008; Wang et al. 2010; Wiegand et al. 2012).

Importance of habitats to explain the species coexistence is widely documented (Harms et al. 2004; Valencia et al. 2004; Gunatilleke et al. 2006; Lai et al. 2009; Wang et al. 2009, Zhang et al. 2012, 2013). For example, an addition of different habitat preferences for members of a species with different genotypes will enable alleles that facilitate survival in those habitats to be maintained as balanced polymorphisms in the population (Levene 1953, Barton 2010). In contrast, the importance of neural or stochastic processes to explain the diversity is extensively discussed. For example, Kimura (1983) suggests that “genetic drift” is a neutral process that can change the allele frequencies over time. Hubbell (1997, 2001) used a similar term called “ecological drift” to explain the macro-ecological patterns in diverse communities (see Leigh (2007) for the evolution of the neutral theory). Some authors suggest that species-habitat association and dispersal patterns (stochastic processes) together can play important roles for species structuring at the intermediate scales (Webb and Pert 2000; Harms et al. 2004; Valencia et al. 2004; Gunatilleke et al. 2006; Seidler and Plotkin 2006; Wang et al. 2009; Qi et al. 2011; Liu 2014; Guo et al. 2017; Liu et al. 2018). However, their stochastic processes differ from Hubbell’s ecological drift or immigration of species from meta-communities. Some suggest that spatial pattern at small scales is due to pathogens and species interactions (Connell 1971; Wright 2002; Gilbert and Lechowicz 2004). Neutral theory suggests stochastic processes at small spatial scales cause species diversity to be unpredictable at small scales.

There are studies that quantify the relative importance of neutral and niche assembly processes to describe the species richness at different spatial scales (Tuomisto et al. 2003; Gilbert and Lechowicz 2004; Jones et al. 2008; Li et al. 2011; Lin et al. 2013; Shen et al. 2013; Lin et al. 2013; Baldeck et al. 2013; Wang et al. 2014; Van der Plas et al. 2015; Yu et al. 2016; Yang et al. 2016; Henkel et al. 2019). However, only few studies highlight the importance of a general mathematical framework to link the processes (e.g. Azaele 2015; Azaele et al. 2016; Matthews and Whittaker 2014; Munoz et al. 2014; Van der Plas et al. 2015).

Here, we asked how many species can be maintained by a neutral model to which one new level of complexity, i.e. space, has been added. For that, we proposed a two-stage model that used species-habitat association and Hubbell’s stochastic processes (i.e. birth-death and dispersal limitation and speciation rates) to describe the species richness. We add habitat filtering into Hubbell’s spatial explicit neutral model to capture ecological processes (i.e. stochasticity, species-habitat association and dispersal limitation) at different spatial scales. Our proposed model uses four parameters, where three parameters come from Hubbell’s neutral model, to generate the species richness at a local community. Our model can generate species richness (and also species composition) for many numbers of habitats and different types of habitats (isolated habitats or connected habitats) by adding one additional parameter (we called this “number of individuals each in habitats”).

Methods

Habitat conservative neutral model

We incorporate Hubbell et al. (2001), Etienne et al. (2005) and Munoz et al. (2007, 2014) to develop our model. We used two stage schemas to generate the regional species pool and the local community. First, we simulate the reference pool with labelling the species. Then immigrants from the regional pool were filtered according to its habitat preference under neutrality. In general, this newly proposed model can be used for the local community with any size with any number of habitats.

We have the following two steps in the model.

Meta-community dynamics

Meta-community size is very large (J_M = 10⁷). Species diversity of the meta-community is driven by two parameters. Fundamental biodiversity number θ = 40 and 50, and speciation rate (ν <  <  <  <  < 1). At the beginning, the number of species in the meta-community is assumed zero. The first individual in the meta-community has the species label s = 1 and individual number i = 1. The second individual of the meta-community can be an offspring of the first individual or new species. If it is an offspring of the previous species, it has the same species label s = 1 and a new individual number i = 2. If it is a new species, then it has a new species label s = 2 and a new individual number i = 2. The rate of new species coming into the regional pool (meta-community) depends on the relative species abundance of the meta-community. When individuals in the regional pool increase, the speciation rate is slowed down. This was incorporated into the model by generating a random number called $R1=\frac{\theta }{\theta +i-1}.$R1 is called the random number generator, where θ is the fundamental biodiversity number, a constant. When individuals in the regional pool increase, i.e. when i increases, then R1 decreases. If the generated random number is less than the R1_, then individual receives a new species label (s = s + 1). Otherwise, a species label is randomly assigned to the individual from the previous species (therefore, species label number depends on the relative species abundance of species in the regional pool at that time). Hence, each individual is tested for whether it is a new species or offspring from previous species.

Local community dynamics

Local community has J number of individuals. The number of habitats (h = 1,  2,  …) in the local community is pre-determined. One of the assumptions in the Hubbell (2001) neutral theory is that the local community is saturated. The number of individuals in each habitat is fixed (J_h), and it is proportional to the area occupied by each habitat $\left({J}_h=J\times \frac{\mathrm{area}\kern0.5em \mathrm{of}\kern0.5em {h}_i}{\sum \limits_{h=1}^k\mathrm{area}\kern0.5em \mathrm{of}\kern0.5em {h}_i}\right)$. It is assumed that in a community with a saturated number of individuals, the number of individuals (J) in an area is proportional to the area sample (A) (i. e. J = ρA). Fundamental immigration number (I₀) is defined by ${I}_0=\frac{m\times \left({J}_i-1\right)}{\left(1-m\right)}$_, where m is the dispersal limitation (0 ≤ m ≤ 1). Dispersal limitation depends on the strength of local community couple with the meta-community. If the local community is fully isolated, then m = 0. We assumed that m depends on the distance between habitat to the meta-community. Here we assumed that distance from habitat to meta-community is the same for all the habitats. Therefore, m is the same for all the habitats. Hence, the fundamental immigration number (I₀) is proportional to the number of individuals in that habitat.

Simulation of local community

The first individual of the local community is randomly drawn from the regional pool (the probability of species s in the local community is proportional to the relative species abundance of the regional pool), and it is placed on one of the randomly selected quadrate (quadrates were selected randomly without replacement). Species label is the same as the meta-community species label.

In the next step, another habitat is randomly chosen without replacement. If the habitat is empty, an individual is randomly drawn from the meta-community and received the same species label. If it is not an empty habitat, there are two possibilities. The second individual in the local community may come from the offspring of the already established species in the local community (if so, then it has the same species label) or species from the regional pool proportional to relative species abundance of the regional pool. However, species from the regional pool may be a new species or the same species of the regional pool.

Whether species is an offspring of the species in that habitat or species from regional species pool depends on the fundamental immigration number. We select a random number from a unit length, and if this random number is less than the $R2=\frac{I_0\left[k\right]}{I_0\left[k\right]+j-1}$ replacement is from the regional species pool and otherwise, it is an offspring of the randomly drawn individual from the local community. When the local community grows, the rate of new species that comes into the local community decreases.

The above three steps (i.e. 3.5–3.7) were performed until the local community is saturated (J = 125,000, i.e. all the habitats are completely occupied by the individuals).

Results

The new model describes the effect of habitat filtering and neutral components on species richness. First, we fixed two parameters and change the other parameter to understand the effects of each parameter on species richness. We used the same meta-community for all three scenarios.

Effects of number of individuals in habitats on species richness of local communities

We fixed the parameters θ and I (see Table 1). Then we defined two habitats with different sizes. We assumed that the habitats are fully saturated and habitat preference is proportional to the habitat size. Therefore, in two habitats with equal sizes, their habitat preferences are equal to 0.5. When individuals in the second habitat are two times greater than the individuals in the first habitat, habitat preferences for the first and second habitats are 1/3 and 2/3, respectively. We observed species richness increases with the number of individuals in the habitat.

Table 1 Local species composition for fundamental immigration number (I = 100) and the fundamental biodiversity number (θ = 50)

Effects of fundamental biodiversity number (θ) on species richness of local communities

We fixed the number of individuals in two habitats (i.e. h₁ = 50, h₂ = 100) and the fundamental immigration number (i.e. I = 100) (see Table 2). Results show species richness increases in each habitat and the local community when fundamental biodiversity number (θ) increases.

Table 2 Simulated species richness for fixed meta-community size (θ = 50, J_M = 1000) and local community size (J_m= 150), I = 100, h₁ = 50, h₂ = 100

Effects of immigration parameter (I) on the species richness of the local community

Fundamental biodiversity number (θ = 50) and the number of individuals in each habitat is fixed (h₁ = 50, h₂ = 100) (see Table 3). Results show species richness increases in each habitat and the local community when fundamental immigration number (I) increases.

Table 3 Simulated species richness values for fixed meta-community size (θ = 50, J_M = 1000) and local community size (J_m= 150), h₁ = 50, h₂ = 100

Effects of number of habitats on the species richness of the local community

Figure 1 shows species rank abundance in the local community when the neutral model is used. Neutral models were fitted for five different sizes of local communities.

Fig. 1

Log species abundance distribution for five local communities with sizes 20, 200, 2000, 20,000 and 200,000 according to neutral model with parameters θ = 40, I = 100

Figures 2, 3 and 4 show species rank abundance in the local community when the habitat-based neutral model is used. Habitat-based neutral models were fitted for five different sizes of local communities (Figs. 5 and 6). Figure 2 shows local community species rank abundance when two equally sized habitats are present. Figure 3 shows community species rank abundance when four equally sized habitats are present in the local community. Figure 4 shows local community species rank abundance when five equally sized habitats are present. All these scenarios θ and I were fixed.

Fig. 2

Log species abundance distribution for five local communities each divided into two habitats with equal sizes. Local community sizes are 20, 200, 2000, 20,000 and 200,000. The two parameters are θ = 40, I = 100

Fig. 3

Log species abundance distribution for five local communities each divided into four habitats with equal sizes. Local community sizes are 20, 200, 2000, 20,000 and 200,000. Two parameters are θ = 40, I = 100.

Fig. 4

Log species abundance distribution for five local communities each divided into five habitats with equal sizes. Local community sizes are 20, 200, 2000, 20,000 and 200,000. The two parameters are θ = 40, I = 100

Fig. 5

Number of species vs. log species abundance for local communities with single habitat size. Local community size is 200,000 with parameters θ = 40, I = 100. This model is equivalent to neutral model

Fig. 6

Number of species vs. log species abundance for local communities with two equal-sized habitats. Local community size is 200,000 and the two parameters are θ = 40, I = 100

We noticed that when local community size increases, both species richness and number of rare species in the community are increased. This observation is very clear for large communities. For example, adding one additional habitat can increase the species richness from 120 to 140 (see Fig. 3 with Fig. 4). When four habitats are present, almost 35 species are singletons (Fig. 7). When five habitats are present, almost 45 species are singletons (Fig. 8).

Fig. 7

Number of species vs. log species abundance for the local community divided into four equal-sized habitats. Local community size is 200,000 and the two parameters are θ = 40, I = 100

Fig. 8

Number of species vs. log species abundance for the local community with five equal-sized habitats. Local community size is 200,000 and the two parameters are θ = 40, I = 100

Habitats decrease the abundance of common species in the local community (Figs.1, 2, 3 and 4). When local community size is fixed (e.g. for 200,000 individuals) dominant species in the neutral model, the two habitat-based neutral models, four habitat-based neutral models and five habitat-based neutral models show 199,238, 99,716, 50,254 and 40,162 individuals, respectively (Figs.1, 2, 3 and 4).

Discussion

Ecology has a rich history on describing species richness and species abundance distributions using simple models (Fisher 1943, Preston 1948, 1962 MacArthur 1957, Sugihara 1980, Tokeshi 1990; Hubbell 2001). Such parsimonious models are called “Occam’s razors” in Ecology (Sober 1990; Baker 2007). These simple models are broadly categorised into purely statistical models, neutral models, and niche-partitioning models (Ferreira et al. 2008; Marquet et al. 2014). However, almost all of these models fail to fully describe the species abundance distribution especially its long tail rare species distribution (Hubbell 2001). Some models (e.g. pure statistical models) lack biological meaning.

Hubbell’s developed three parameters fully stochastic model that can describe the long tail rare species abundance distribution in the tropical forests (Hubbell 2001). One of the core assumptions of the neutral theory is that all the species are equivalent with respect to their per capita birth rates, death rates and speciation rates (Hubbell 2001). Although many argued that the ecological equivalence of species is not realistic, this reality simplification assumption helps to form a null model for ecologist. Hence, the falsifiability of this assumption left more lights for ecologist and the assumption of “ecological equivalence” cannot be regarded as a weakness of the neutral theory (Alonso et al. 2006).

Some authors relaxed this “ecological equivalence” assumption so that they could study the species coexistence in more general (Zhou and Zhang 2008; He et al. 2012, Zhang et al. 2012, Zechen et al. 2012). For example, Zhou and Zhang (2008) show that subtle changes to species fecundity drastically reduced the species richness and changes to species abundance distribution. They argued that even strong dispersal limitations cannot compensate for the subtle changes to species fecundity to maintain the species richness. In a similar direction, Zechen et al. (2012) show that, to observe neutral-like species abundance, distribution patterns under slightly asymmetric competition need severe dispersal and recruitment limitations. Therefore, dispersal and recruitment limitations are not sufficient to compensate for the difference in species fitness (Zechen et al. 2012). Hence, more ingredients need to build a unified theory. However, He et al. (2012) and Zhang et al. (2012) show the assumption of “ecological equivalence” when species birth rates and death rates are tradeoff. They show that “generalised neutral models” (i.e. species birth-death rates are tradeoff) yield the same macro-ecological patterns as the neutral model. According to our view, He et al. (2012) and Zhang et al. (2012) define ecological equivalence of species at fitness level.

The second assumption of the neutral theory is that the community undergoes a zero-sum game (Hubbell 2001). Volkov et al. 2003; Etinne et al. 2007 found that relaxing the zero-sum assumption does not affect the species abundance distribution pattern. However, we used the “zero-sum game” to relax the ecological equivalence from the local community to habitat. According to the zero-sum game (described by Hubbell 2001), the local community is saturated. Therefore we can assume that at some point, all the habitats are fully occupied by individuals. Hubbell (2001) describes for saturated communities, community size is proportional to the area of the local community. Therefore, we can assume that, in the long run, habitat preference is proportional to the number of individuals in the habitat (habitat size).

Therefore, it is very important to consider the number of individuals in each habitat as a parameter to reconcile the neutral theory with the habitat associations. We think that adding this parameter into the neutral theory increases our understanding of the underlying processes that shape the ecological communities without disturbing the degrees of freedom of the neutral models. The additional parameter allows us to place species into habitats with different probabilities. Hence, it allows us to relax one of the core assumptions in the Hubbell’s neutral model: “species are equivalent”. The neutral theory allows total freedom for species to move in the local community. However, our model relaxed this assumption and restricted species range so that off-springs can move only in parent habitat. However, our model still allows species to move within the habitat. Therefore, the independence of species is achieved at the habitat level.

This additional parameter also helps to construct the spatial arrangement of the local community. The neutral theory does not incorporate the exact spatial arrangement of the species in local communities (Armesto et al. 1986; Condit et al. 2000; Wiegand et al. 2012). We agree that such exact spatial arrangements are not important for the neutral theory to describe the aggregated statistics such as macro-ecological patterns (i.e. species richness and species abundance) at the large spatial scales. However, one major criticism is that a large amount of information loss in neutral models is due to aggregated statistics. Therefore, the neutral theory cannot explain the species composition. In this regard information on the spatial arrangement of species is very important to increase the understanding of the underlying processes in ecological communities at different spatial scales.

Our model shows the habitat increases the number of rare species in the community. Hubbell (2013) found that almost all of the tree diversity comes from 0.1% to 4.1% of individuals in tropical forests, which are rare species. Tropical forests provide shelter for thousands of rare species, and most of them are singletons. Therefore, conservation of habitats is uppermost important to maintain the species diversity of these tropical forests to slow the extinction of rare species and increase its numbers. Our results indicate species abundance distribution pattern becomes flat and much wider when the number of habitats increases. On the other hand, this model reduces the abundance of dominance species, which allows other species to increase their numbers and delay the competitive exclusion.

Our model has limitations. First, we defined habitat as spatial regions. Further, we assumed that the competition for space is equal for each and every individual (i.e. similar to Hubbell’s neutral model). Although space is one of the important and limited resources that all the species compete for, there are several other variables one can incorporate to define habitats more broadly. This may include an abundance of pollinators, predators, prey and diseases. Such density-dependent processes may act on species richness and diversity at different spatial scales. For example, species-specific pollinators can increase the abundance of most abundance species by decreasing the abundance of rare species. On the other hand predators (e.g. herbivores), diseases and pathogens can help to decrease the abundance of dominant species, thus allowing rare species to increase its numbers. Thus, incorporating such variables bring our model much closer to reality. However, adding more variables into a model increases the model complexity. As we mentioned before, ecology has a rich history of the parsimonious model to explain the complex natural phenomena and our idea was to add space into Hubbell’s neutral model to build a simple model.

Second, we defined habitat as a space where a set of species coexists together. This is a definition we adapted from studies done in large forest plots such as BCI, Sinharaja, etc. (Harms et al. 2004, Gunatilleke et al. 2006). We accept that defining habitat for each species to allow the opportunity to test our model to its logical extremes and estimate the range of possible effects better. Especially, such constraint allows species-species (interspecific) interaction to become zero. However, in the beginning, we assumed that the local community is saturated and the number of individuals in a local community is proportional to its area. Thus, allowing habitat for each species causes habitat size to be proportional to the number of individuals. In other words, we can see that habitat sizes in a local community follow the exact same pattern as relative species abundance distribution. On the other hand, we doubt that adding habitat for each species might cause habitats to be completely overlapped with the niches due to zero interaction among species.

Conclusions

Our model suggests that “the number of individuals in habitats” can enhance our understanding on species assemblages of the local communities. Our findings support the idea that the “range sizes” are different for species in tropical forests. Therefore, rare species may be restricted to a certain habitat and has a greater risk of extinction (Hubbell 2013).

We conclude that the new model with only four parameters, theta (θ), immigration parameter (I), speciation, and the number of individuals in the habitats, helps to understand the species richness in habitats. Changes to model parameters can construct species richness of any local community with different numbers of habitats.