Introduction

While the mechanisms that generate large-scale diversity patterns are still under investigation, successional dynamics and conspecific negative density dependence (CNDD) are often considered key processes for species assembly at local scales1,2,3, though other factors, such as resource specialization and immigration are not less important.

Succession is the sequential change in forest composition and structure following disturbances4. It can operate on a wide range of spatial and temporal scales, from gap-phase dynamics generated by windthrow or other locally confined mortality events to large-scale calamities such as fires, hurricanes, floods, and logging5,6,7,8. Roughly two-thirds of the global forested area is classified as secondary forest9, making the succession study particularly relevant to conservation problems. This study focuses on sequential changes that occur within a time frame of the intervals between disturbances and in the same order of magnitude as the life span of the longest-lived plants in the successional sequence4. Sequential implies that once a species is excluded from the canopy, it will not have access again unless a new disturbance occurs. The diversity of strategies in the successional sequence, referred to as successional diversity, is often linked to life-history tradeoffs, such as the tradeoff between tree high-light growth and low-light mortality or between growth and reproduction10,11. Recent modeling studies demonstrate how these tradeoffs can maintain many species similar to the competition-colonization tradeoff12,13,14.

Conspecific negative density dependence (CNDD) could affect any life stage of a tree, but most examples, including this study, refer to early-life biotically-mediated feedbacks that reduce the establishment of offspring near conspecific adults, including their mother2,15. Empirical evidence of CNDD operating in forest ecosystems is extensive, and although this mechanism was initially proposed to explain the local diversity of tropical forests, it has also been observed to operate at higher latitudes16,17,18,19,20,21,22,23,24.

There is a pressing need for mathematical frameworks integrating multiple coexistence mechanisms, as the majority consider only a single mechanism, potentially hindering progress for decades. Such frameworks are necessary to weigh the relative importance of these mechanisms and assess their interactions, which only integrative approaches can achieve. For example, CNDD is a mechanism that stabilizes and promotes coexistence by increasing intraspecific competition25,26. However, theoretical models that have investigated the importance of CNDD in maintaining diversity reached discordant conclusions27,28,29. One of the problems in doing that is that the application of this principle to predict forest diversity has been limited by the difficulty of finding a way to realistically propagate a local process that operates at early life stages across the lives of long-lived organisms that vary in size by orders of magnitude and that compete strongly for limited resources30,31.

Introducing CNDD in neutral communities provides a simplified context for understanding the potential role of density-dependent mechanisms32,33. However, this simplified scenario may not fully capture the complexity of real-world ecosystems where competitive interactions play a significant role28,34. In more realistic cases where communities undergo competitive exclusions, the importance of CNDD may vary depending on various factors such as species traits, environmental conditions, and the intensity of competition. Competitive exclusions can result in the dominance of certain species, potentially influencing the effectiveness of CNDD in promoting diversity.

A comprehensive forest diversity theory should also be framed from an evolutionary perspective to explain how plant traits are selected from a broader species pool and assembled in local communities, how species interactions drive trait distribution, and how these dynamics generate feedback on larger scales. A practical approach is to connect local forests to a metacommunity, where species traits evolve over extended evolutionary timescales.

This study proposes incorporating conspecific negative density dependence (CNDD) into a vegetation dynamic model (VDM) linked to a metacommunity. Such a process-based model can be developed using environmentally dependent physiological functions, including photosynthesis, respiration, and plant carbon economy, thereby linking life history strategies with plant traits. By doing so, we can integrate CNDD into forest demography, offering valuable insights into its interactions with other key forest processes, such as competition for light, speciation, and community structuring through successional dynamics.

The investigation of the maintenance of functional diversity in current VDMs is an important emerging field of research35,36,37,38. Previous studies using forest simulators have shown that including CNDD facilitates coexistence in these models39,40, thus, the questions are how much and how much. Detto et al.12 proposed the first analytically tractable model that predicts successional diversity and yet is consistent with individual-based forest stand simulators and the ecosystem components of some Earth System Models41,42. That study identified four possible life-history tradeoffs, each created by interspecific variation in a single plant trait and each able to maintain high successional diversity, including maximum leaf photosynthetic capacity (growing fast in high light vs. surviving in low light), fraction of carbon allocated to reproduction (growth vs. reproduction), an allometric parameter relating crown size to stem biomass (growing up vs. growing out), and age at first reproduction (shifting carbon from growth to reproduction later vs. earlier in life). In this study, we ask how all these tradeoffs interact to maintain functional diversity, i.e., the number of species with a unique combination of traits. If \(N\) species can be maintained with one tradeoff, are \(2N\) maintained with 2 tradeoffs, or even \({N}^{2}\)? Moreover, how does the addition of CNDD increase diversity? Does increased diversity, if any, entail greater or lesser successional diversity? Does CNDD alter the rate at which species are assembled on evolutionary scales? To answer these questions, we extend the model in Detto et al.12 to include CNDD and interspecific variation in multiple parameters and, thus, multiple successional tradeoffs. The goal is to derive a mathematical and graphical theory that can explain coexistence and predict richness and related macro-ecological functions for a given pool of species and as a function of the strength of CNDD.

Methods

Successional forest model

Successional dynamics are usually studied with forest patch models, which have a well-understood theoretical and empirical basis43,44, but often require stochastic simulations or numerical integrations4,7. The model in Detto et al.12 is based on a set of integro-partial differential equations governing the age- and size-structured dynamics of a multi-species community of trees in an infinite successional mosaic of patches in a steady state. However, this model can be accurately defined and analyzed with simple methods that are primarily verbal and graphical.

Trees are assumed to interact only with neighbors within the same successional patch by overtopping and shading one another, and patches are dynamically coupled through seed dispersal. Both reproduction and growth of plants in the understory are assumed to be negligible. Patches are subjected to random disturbances, which kill all the trees but leave a fraction of the seeds or seedlings, resetting succession.

Recruits close the canopy early in succession, so post-disturbance regeneration is governed by individuals present in the site as seeds or seedlings when the disturbance occurs. This implies that, in each patch, all canopy individuals of a given species at any one time have identical sizes (but not individuals that fall into the understory and effectively stop growing). Moreover, all canopy individuals started growing simultaneously – when the most recent disturbance removed the canopy that was shading them.

Crucial to the model are the patch dynamics between disturbances. These dynamics are marked by sequential exclusions from the canopy as taller plants overtop smaller plants and send them into the understory45,46. As succession progresses, growth leaves room for fewer and fewer individuals in the canopy.

Reproduction and density dependence effects

Reproduction of canopy individuals is continuous, size-dependent, and species-specific, and dispersing seeds have an equal probability of landing in any patch where they accumulate. Plants allocate a fixed fraction of net primary productivity, \({{\rm{\varphi }}}\), to seeds. Seeds may germinate immediately, after which seedlings suffer understory rates of mortality and growth until the next disturbance, or they may suffer seed mortality and germinate at the start of the next disturbance. Disturbances kill all the trees but leave a fraction of the seeds or seedlings, resetting succession.

The seeds and seedlings are subjected to conspecific density-dependence (CNDD), such as those associated with specialized natural enemies. CNDD can operate at seed production or at the seedling stage. In either case, as the number of adult trees increases, fecundity or seedling survival in the understory decreases. We express the mean density of recruits of species-\(i\) in the understory \({R}_{i}\), as:

$${R}_{i}=\frac{{f}_{i}}{{m}_{i}}\frac{{p}_{i}}{1+{{\rm{\beta }}}{p}_{i}}$$
(1)

where \({p}_{i}\) is the probability of species-\(i\) to have access to the canopy space in the forest mosaic, \({f}_{i}\) and \({m}_{i}\) are the fecundity rate per unit of crown area and the seedling mortality in the understory with no CNDD, and \({{\rm{\beta }}}\) is the Berverton-Holt CNDD parameter, which ranges as the strength of CNDD increases from zero to infinity47. \({f}_{i}\) and \({m}_{i}\) can be species-specific depending on the trade-off type. For example, under the growth-fecundity trade-off, species vary in the number of seeds produced per unit of crown area. However, trade-offs like growth ‘up’ vs. growth ‘out’ or fast growth in light vs. low survival in shade remain constant across species.

Thus, the model permits density-dependent mortality caused by shading or CNDD. It also permits density-independent mortality in the understory for all stages, including seeds, seedlings and adults, and density-independent disturbance that kills all individuals except for a fraction of seeds or seedlings. However, the model ignores all other density-independent mortality of canopy trees, which is probably smaller than any of the sources of mortality that are included.

Plant light interactions

Exclusions from canopy access can occur in two ways. In the first, early successional dominants have the fastest initial height growth rate but eventually either slow their height growth or die and cede the canopy to slow-growing shade-tolerant species. In the second, growth trajectories never cross (Fig. 1b). Early successional specialists accumulate lifetime reproductive success faster than late successional dominants but have slower height growth throughout their lives. This kind of tradeoff can be created in many different ways. For example, because carbon allocated to stem tissue does not go to seeds or flowers, interspecific variation in height allometry causes a tradeoff between the accumulation rate of lifetime reproductive success and the height growth rate. Interspecific variation in reproductive allocation does the same thing. A less obvious example causes succession from shade-tolerant species to giant, long-lived shade-intolerant species. Although less familiar, this type of successional dynamics can be quite common and leads to canopy dominance by long-lived or giant pioneers48,49, including large, majestic, exceptionally old trees50. Succession from shade-tolerant species to long-lived shade intolerants can be caused simply by variation in the full sun photosynthetic rate and associated leaf respiration. In this paper, to maintain mathematical tractability, we focus only on cases in which height growth trajectories do not cross and thus do not consider the familiar examples of fast-growing short-lived pioneers that dominate the canopy early but then give way to slower-growing species.

Fig. 1: A perfect tradeoff emerges from competition dynamics.
figure 1

a The relationship between break-even time and height growth in a metacommunity of species that vary randomly in the fraction of net plant productivity allocated to reproduction (\({{\rm{\varphi }}}\)) and the allometric parameter regulating crown size \((a)\) to biomass (\(b={{\rm{\phi }}}{a}^{\frac{c}{c+1}}\)), shown in the inset. The species that will potentially coexist in the forest patch mosaic constitute the left-upper envelope of the scatterplot and represent the fastest height growth for a given break-even time (red circles). The inset shows the relationship among the physiological traits \(f\) and \({{\rm{\phi }}}\). b Height as a function of age for the potential coexisting species in (a). Note that the trajectories never cross.

Full size image

The “perfect” tradeoff

Species differ in their height growth rates (\(g\)) and break-even time (\(r\)), which is the age at which an individual growing in full sun achieves a lifetime reproductive success (\({LRS}\)) of one, assuming no CNDD51. An individual reaches \({LRS}=1\) when at least one of its offsprings survive a disturbance event and begins growing in full sun. Typically, allocation tradeoffs constrain tall species that grow fast in height by taking longer to reach \({LRS}=1\)52,53,54,55. A variation in a single trait or a combination of traits can achieve this tradeoff. To illustrate this concept, we consider a pool of species that varies in two independent traits: the fraction of net plant productivity allocated to reproduction, \({{\rm{\varphi }}}\), and an allometric parameter \({{\rm{\phi }}}\) that regulates the crown size (\(a\)) relative to structural biomass (\(b\)), i.e., \(b={{\rm{\phi }}}{a}^{\frac{c}{c+1}}\), where \(c\) is an allometric exponent. Random variation in these two traits generates a tradeoff between \(g\) and \(r\) but with some variability (Fig. 1a). For a given break-even time, there is an optimal trait combination that confers fast growth (red cycles). In a competitive environment, these species will overtop the suboptimal species (points below the envelope) before they can reach \({LRS}=1\), making it impossible for them to compete. To see that, imagine a two-species system with \({g}_{1} > {g}_{2}\). At equilibrium, the highest species-1 closes the canopy before reaching \({LRS}=1\) (i.e., at time less than \({r}_{1}\)). A necessary condition for the coexistence of species-2 is to reach \({LRS}=1\) before species-1 closes the canopy, i.e., \({r}_{2} < {r}_{1}\). In practice, this implies that a species that grows slower than its competitor must reach \({LRS}=1\) faster: it cannot be defective in both aspects. Consequently, for the successional model to maintain coexistence, species must satisfy the following necessary (but not sufficient) conditions:

$${r}_{1} \, > \, {r}_{2}\, > \, \ldots \, > \, {r}_{n}\,,\,{g}_{1} \, > \, {g}_{2}\, > \, \ldots \, > \, {g}_{n}$$
(2)

Inequalities (1) represent a “perfect” tradeoff and reduce a multidimensional life-history trait space to a single axis, even if there is no physiological or structural tradeoff between \(g\) and \(r\), with early successional species having small \(r\) and late successional big \(r\). The condition for coexistence is more stringent and must be determined by the exact sequential times the canopy closes (see below). This is why the condition is necessary but not sufficient.

Population equilibrium and coexistence condition

At equilibrium, the population of species \(i\) is governed by the following equations (Supplementary Appendix A):

$${p}_{i}={t}_{i-1}-{t}_{i}$$
(3)

where \({t}_{i}\) is the time when species-\(i\) individuals start to be overtopped and \({t}_{i-1}\) is the time when the last individual enters the understory (Fig. 2). At equilibrium, \({t}_{i} < {r}_{i} < {t}_{i-1}\), because otherwise, all individuals of species-\(i\) would have \({LRS} > 1\) or \({LRS} < 1\), so the system would not be at equilibrium (Fig. 2). \({t}_{i}\) are computed from the condition that the canopy is closed, i.e. the sum of the crown area of all individuals in the canopy equals patch size. To derive (3), time is transformed so the patch-age distribution becomes uniform between zero and one.

Fig. 2: Lifetime reproductive success, LRS, is the number of offspring that survive disturbances and are produced by an individual with full access to light until age t.
figure 2

At \({r}_{i}\), the individual of species-\(i\) has produced enough seeds during his life that, on average, one will turn into a mature tree (i.e., \({LRS}=1\)). In our model, \({LRS}\) is equal to \({\left(t/r\right)}^{c+1}\), where and \(c\) is an allometric constant shared by all species, which depends on the scaling exponent between tree mass and crown area (fixed at 1.5). The niche shadow is the interval that cannot be invaded by shorter competitors and ends at the time (\({t}_{i}\)) species-\(i\) starts being overtopped by taller competitors. The safety distance represents how far break-even time is from the niche shadow of the taller competitor (\({t}_{i-1}\)).

Full size image

The interval \({{{\rm{\xi }}}}_{i}={t}_{i-1}-{r}_{i}\) defines the safety distance of species-\(i\) from the canopy exclusion time set by taller competitors (\({t}_{i-1}\), Fig. 2). The further a species’s \(r\) from this line, the more patches will be occupied by individuals with \({LRS} > 1\), leading to a higher population density. Conversely, the closer a species’s \(r\) is to \({t}_{i-1}\), the lower the population density is, the greater the risk of extinction (in a stochastic system). The interval \({{{{\rm{\varepsilon }}}}_{i}=r}_{i}-{t}_{i}\) represents the non-invasibility interval, or niche shadow14, of species-\(i\). No species with a break-even time \(r\) within this interval can successfully invade a community of taller species at equilibrium. Because competition is asymmetric, the condition for stable coexistence given by the mutual invasion criterion56 of species-\(i\) in a community of taller competitors reduces to:

$${r}_{i}\, < \, {t}_{i-1}$$
(4)

The metacommunity and evolutionary adaptation

The species that can coexist in the forest mosaic are selected from a pool of \({J}_{M}\) species that satisfy the inequalities given in Eq. (1). The number and composition of species in the metacommunity vary on evolutionary time scales, so they can be considered ‘frozen’ regarding the dynamics in the forest mosaic57. Most of our results are obtained assuming a uniform distribution of break-even times among the species in the metacommunity. This distribution is bounded between \({r}_{\min }\) and \({t}_{0}=1\), representing the minimum and maximum age that is biologically feasible given the disturbance regime (i.e. an individual growing exclusively in full sun will have expected \({LRS} > 1\)). These assumptions allow us to derive species richness, species abundance distribution, and the distribution of break-even time analytically (Supplementary Appendix B).

In the final analyses of the paper, we introduce an evolutionary model with genetic drift and trait selection that can generate an evolutionarily stable distribution of break-even times to study how fast richness is generated, how species composition changes over evolutionary time, and how CNDD affects these dynamics. We will also explore the differences between mutations in a single trait and two traits to generate diversity using allocation to reproduction and the allometric parameter regulating crown size to biomass as physiological traits. In this model, species evolve from a common ancestor. A speciation event marks the origin of a species with a new \(r\) and \(g\). Physiological traits evolve as a Gaussian jump process on a continuous trait axis58. Selection is based on competition outcomes in the local community as new species may successfully establish and cause resident species to be excluded. Species that fail to invade successfully or are excluded from the local community risk extinction in the metacommunity, creating a feedback loop between local and metacommunity dynamics.

Statistics and Reproducibility

As a case study, we used the forest inventory data of the 50-ha plot of Barro Colorado Island to fit the species abundance distribution model. The inventories contain all trees with a diameter at breast height (DBH) greater or equal 1 cm, mapped with coordinates relative to the plot, with species identified and DBH measured59. Since our model does not track individuals in the understory, we selected trees based on their canopy access by calculating the fraction of each tree’s crown area with direct light exposure. We assumed the crowns to be perfect circles centered at the stem location, with crown radius and height estimated from site-specific allometric equations42,60. This approach allowed us to determine the extent to which each tree’s crown area is not overlapped by a taller neighboring tree. We selected only trees with at least 25% exposure and DBH\(\ge \)10 cm, resulting in an average of 7163 individuals and 203 species per census across eight censuses conducted in 1982, 1985, 1990, 1995, 2000, 2005, 2010, and 2015.

The frequency histogram was computed by the modified Preston’s binning method given by the following formula61:

$${F}_{0}=\frac{1}{2}{s}_{1}$$
$${F}_{1}=\frac{1}{2}{s}_{1}+\frac{1}{2}{s}_{2}$$
$${F}_{2}=\frac{1}{2}{s}_{2}+{s}_{3}+\frac{1}{2}{s}_{4}$$
$$\ldots $$
$${F}_{i}=\frac{1}{2}{s}_{{2}^{i-1}}+{s}_{{2}^{i-1}+1}+\ldots +{s}_{{2}^{i}-1}+\frac{1}{2}{s}_{{2}^{i}}$$

where \({s}_{i}\) is the number of species with \(i\) individuals. The frequency histogram was calculated for each of the eight available censuses and then averaged. We fit the simple model without CNDD. We compute the number of species with \(i\) individuals by discretizing Eq. (15) in Supplementary Appendix B as:

$${s}_{i}=\frac{{{\rm{\lambda }}}}{2\sqrt{1-2i\Delta z}}{e}^{-{{\rm{\lambda }}}{r}_{-}\left(i\Delta z\right)}\Delta z+\frac{{{{\rm{\lambda }}}}^{2}}{4}\left[{Ei}\left({{\rm{\lambda }}}{r}_{-}\left(i\Delta z\right)\right)-{Ei}\left({{\rm{\lambda }}}{r}_{+}\left(i\Delta z\right)\right)\right]\Delta z$$

where \({E}_{i}\left(\cdot \right)\) is the integral exponential function and \({r}_{\pm }(z)\) are the solution of \({2r}^{2}-2r+z\).

This model has only two free parameters: \({{\rm{\lambda }}}\) and the scaling factor \(\Delta z\). We fit these parameters by minimizing the sum of square errors between observations and model:

$${\min }_{{{\rm{\lambda }}},\Delta z}{\Sigma }_{i}{\left({F}_{i}^{{obs}}-{F}_{i}^{{\mathrm{mod}}}\right)}^{2}$$

using the function fmincon in Matlab (R2019a).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Results

A graphical method

The model belongs to a family of tractable competition models of ecological assembly whose stable equilibrium is independent of changes in invasions’ order, size, and timing62. Figure 3a shows a dynamic simulation of six species quickly reaching an equilibrium where two species go extinct (gray lines). The determination of species that coexist in the forest mosaic from the metacommunity can be illustrated with a graphical method. The graph in Fig. 3b represents the \({LRS}\) of plants with full access to light as a function of the patch age (time since last disturbance), where each line corresponds to each of the six species in the metacommunity. The break-even times (dots) are the age when \({LRS}=1\) (horizontal dashed line). As mentioned above, the age axis is transformed so the maximum patch age is \({t}_{0}=1\). Because of the tradeoff between break-even time and height growth, species on the graph are ranked from shortest to tallest—the method proceeds from right to left.

Fig. 3: A graphical method to explain coexistence of multiple species.
figure 3

a Numerical simulation of the dynamic system shows that the four coexisting species (black lines) persist while the two go extinct (gray lines). b Graphical method to determine which species from the metacommunity stably coexist in the forest mosaic. Each black curve represents a species-specific \({LRS}\). A competitively excluded species (light gray) cannot persist either because: (1) it has a break-even time (light gray dot) that is greater than \({t}_{0}\), which means that its expected \({LRS}\) is less than one, even without any competitors present, or (2) it cannot reach its break-even time (gray dot) before a taller species has already closed the canopy over it (at the vertical red line immediately to the left of the gray dot).

Full size image

Only species with \(r\) smaller than \({t}_{0}=1\) can be established, so all species with \(r\) greater than \({t}_{0}\) are rejected a priori (e.g., the last species on the right of Fig. 3b). The first species to establish, ranked 1, which is also the tallest species in the mosaic, is the species with the greatest \(r\) less \({t}_{0}\) (\({r}_{1}\) in Fig. 3b). At equilibrium, species-1 will close the canopy at time \({t}_{1}\)(vertical red line), setting a new boundary for the next species in the successional sequence. The next species that will persist is the species with the greatest \({r < t}_{1}\) (\({r}_{2}\) in Fig. 3b). At equilibrium, the first two species will close the canopy at time \({t}_{2}\), which is the new boundary for the third species in the sequence and so on. This simple recursive algorithm can be repeated until no more species in the metacommunity are left. Note that, in the graphical example \({t}_{2}\) precludes one species from establishing (the species depicted by a gray curve between the curves for species-2 and species-3).

Species richness and the role of CNDD

The safety distance \({{{{\rm{\xi }}}}_{i}=t}_{i-1}-{r}_{i}\) is a random variable that depends solely on the distribution of \(r\) in the metacommunity. For example, if the distribution is uniform, \({{\rm{\xi }}}\) is exponentially distributed with rate \(\frac{{J}_{M}}{1-{r}_{\min }}\). The niche shadow \({{{{\rm{\varepsilon }}}}_{i}=r}_{i}-{t}_{i}\) is also a random variable and can be derived from the population Eq. (3). Interestingly, the niche shadow and the safety distance are approximatively linearly proportional (b and Supplementary Appendix A):

$${{{\rm{\varepsilon }}}}_{i}\cong \frac{c+1-2{{\rm{\beta }}}{r}_{i}}{c+1+2{{\rm{\beta }}}{r}_{i}}{{{\rm{\xi }}}}_{i}$$
(5)

This relation allows us to compute several macroecological properties using the renewal process analogy (Supplementary Appendix B). For example, in the case without CNDD (\({{\rm{\beta }}}=0\)), the intervals between break-even times in the mosaic are twice those in the metacommunity because \({{{\rm{\varepsilon }}}}_{i}\cong {{{\rm{\xi }}}}_{i}\) (Fig. 5b). Consequently, the number of species that coexist in the mosaic is roughly half the number in the metacommunity \({J}_{M}\). For sufficiently large \({J}_{M}\), the general formula for successional diversity subjected to CNDD can be approximated by a linear function (Supplementary Appendix B):

$$s\cong 1+\left(1+{{\rm{\beta }}}\frac{{r}_{\min }+1}{c+1}\right)\frac{{J}_{M}}{2}$$
(6)

Equation (6) and Fig. 4a show that \({{\rm{\beta }}}\) interacts with \({J}_{M}\) in promoting diversity: the larger the species pool, the larger the increase in diversity with CNDD. Another consequence of Eq. (5) is that CNDD reduces recruits more for late-successional species (Fig. 4b), thus influencing species distribution along succession, as we will see later.

Fig. 4: Effect of CNDD on successional diversity and recruitment.
figure 4

a Analytical results (lines, Eq. (6)) of successional diversity, as a function of CNDD parameter \({{\rm{\beta }}}\) for different metacommunity size \({J}_{M}\). b The expected density of recruits as a function of break-even time declines with CNDD. This decline is stronger for late-successional species. Near \(r=1\) recruits drop to zero because the break-even time is too close to the disturbance interval. Other parameters \({J}_{M}=100\), \(c=1.5\), and \({{\rm{f}}}/{{\rm{m}}}=1\).

Full size image

The graphical method offers an intuitive explanation of how CNDD operates in this system. In Fig. 5a, species-i (black line) closes the canopy at time \({t}_{i}\), which, in the absence of CNDD, can be found by taking \({{{\rm{\varepsilon }}}}_{i}\) roughly equal to \({{{\rm{\xi }}}}_{i}\). This boundary prevents the earlier successional species (gray line) from establishing because its break-even time (gray dot) is greater than \({t}_{i}\). However, in the presence of CNDD, \({{{\rm{\varepsilon }}}}_{i}^{{\prime} }\) is smaller than \({{{\rm{\xi }}}}_{i}\), Eq. (5), shifting \({t}_{i}\) to \({t}_{i}^{{\prime} }\) and allowing the earlier successional species to coexist.

Fig. 5: A graphical method explains the role of CNDD in reducing competitive exclusion and generating diversity.
figure 5

a Without CNDD the species represented by the gray line is excluded because its break-even time (gray dot) is greater than \({t}_{i}\) With CNDD, \({{{\rm{\varepsilon }}}}_{i} < {{{\rm{\xi }}}}_{i}\) and \({t}_{i}\) shits to \({t^{\prime} }_{i}\) (blue line), allowing coexistence. b Simulations show that without CNDD \({{\rm{\varepsilon }}}\approx {{\rm{\xi }}}\) (red); with CNDD \({{\rm{\varepsilon }}} \, < \, {{\rm{\xi }}}\). Dashed lines represent Eq. (5) with \({r}_{i}=1/2\) and \(c=1.5\).

Full size image

Species break-even time and abundance distributions

Given the properties of \({{\rm{\xi }}}\) and \({{\rm{\varepsilon }}}\), we can now compute the probability distribution of break-even times and species abundance using the renewal process analogy (Supplementary Appendix B):

The distributions for different values of \({{\rm{\beta }}}\) are shown in Fig. 6. Without CNDD, the break-even time distribution is flat and increases sharply near \(r=1\). With CNDD, the peak at \(r=1\) is reduced, and the probability increases linearly with \(r\) (Fig. 6a).

Fig. 6: Species distributions along the successional and abundance axes.
figure 6

a Break-even time and (b) species abundance probability distributions for different values of CNDD parameter \({{\rm{\beta }}}\). as indicated by the color bar (\({J}_{M}=25,{c}=1.5,\,{r}_{\min }=0.2\)). c Fit SAD model to tree species abundance data from the Barro Colorado Island plot. Only trees with access to the canopy with DBH\(\ge \)10 cm were selected for the analysis (see Statistics and Reproducibility). The frequency distributions are plotted using Preston’s binning method61. Fitted parameters: \(\frac{{J}_{M}}{1-{r}_{\min }}=360\) and scale parameter \(\frac{m}{f}=7.2\times {10}^{-3}\) ha.

Full size image

The species abundance distributions (Fig. 6b) show a clear peak that shifts on the left with increasing CNDD (the x-axis is logged so that a small shift might represent a relatively large difference in abundance). This is consistent with the increase in richness; their abundance decreases as more species are fitted into the community. Note also that the shape of the distribution is minimally affected by CNDD, for example, the skewness increases by 7% passing from \({{\rm{\beta }}}=0\) to \({{\rm{\beta }}}=0.75\).

As an example, we fitted the SAD model without CNDD to tree species abundance data from the Barro Colorado Island plot (Fig. 6c). Because our model tracks only individuals with access to the canopy, we include only trees with at least 25% canopy access using a crown overlapping algorithm based on relative position and site-specific allometry (see Supplementary Appendix B for details).

Adaptive evolution

So far, we have connected the local community with a static metacommunity. To explore how traits evolve over time and how they interact with the local community dynamics, we run several simulations using an evolutionary model.

Each simulation starts with a single species at \(r={r}_{\min }\). No suboptimal strategies exist for the single trait mutation (allocation to reproduction), so each new species can potentially compete for a successional niche. However, in the first stage, a newly speciated species can establish itself only if \(r\) is larger than the resident, in which case, the new species will replace the resident. This continues, causing \(r\) to “walk” to the right until it is sufficiently close to one (i.e., \(r > {r}_{\min }+{{\rm{\varepsilon }}}\)), making coexistence between the resident and a higher-\(r\) species possible. From this point, new species begin to invade to the right, left, or between the resident species with smaller \(r\) than either resident. Species packing continues indefinitely until the simulation is stopped, with denser packing on the late-successional. Two factors generate these dynamics. First, speciation to the right of the tallest species is always successful. Second, as the species get closer to one, their limiting similarity decreases, approaching asymptotically zero as the number of species becomes infinite12. Equation (3) shows this, because \({{{\rm{\varepsilon }}}}_{i}\) is proportional to \({{{\rm{\xi }}}}_{i}\). As more and more random species are drawn from the metacommunity, the \({{\rm{\xi }}}\)’s approach zero and so the \({{\rm{\varepsilon }}}\)’s also approach zero, which means the spacing between coexist species also goes to zero. Figure 7a shows the increasing number of species as function of evolutionary time for different CNDD parameters \({{\rm{\beta }}}\). Because CNDD reduces limiting similarities and increases the probability of successful speciation, species packing is faster in the presence of CNDD. For a given evolutionary time, species richness increases linearly with the CNDD parameter (Fig. 7b).

Fig. 7: Evolutionary dynamics.
figure 7

a Adaptive evolutionary dynamics show species packing increases faster with CNDD. Simulations start from a single ancestor with \(r={r}_{\min }\) and speciation occurs as a mutation in a single trait (allocation to reproduction) evolving as a jump process. The success of establishment in the local community and extinction in the metacommunity are based on the outcome of the competition model, which is evaluated after each speciation event. b Species richness as function of the CNDD parameter for an evolutionary time equal to 150 speciation events. Each point represents an independent simulations. Linear regression, (\(n=600\)) shown for reference (solid line). c Comparison of species richness as function of metacommunity size with a single-trait evolution (allocation to reproduction) and a multi-trait evolution (allocation to reproduction and allometric parameter relating crown size to biomass). Each point represents an independent simulations. Linear regressions (\(n=600\)) shown for reference (solid lines).

Full size image

In the final analysis, we compare simulations for single-trait mutations with a multi-trait. Surprisingly, the multi-trait case does not sustain more diversity than the single trait for a given metacommunity size but is slightly lower (Fig. 7c). The reason is that the multi-trait case generates suboptimal strategies that cannot compete for a successional niche (Fig. 1a).

Discussion

Our study integrates competition for limited resources, specialized natural enemies, and evolutionary dynamics, broadly appreciated as among the most important processes governing plant community composition, structure, and function, into a unified theoretical framework. It elucidates the role of conspecific negative feedback operating at early life stages in maintaining successional diversity in forest ecosystems with a strong competitive hierarchy.

Despite the diversity of life-history strategies determined by many plant traits and parameters, coexistence ultimately relies on species adhering to a single fundamental axis. By simplifying ecological complexity into a tradeoff between height growth and break-even time, rather than increasing the number of niches, which is usually thought to create opportunities for coexistence, the model suggests that optimizing these two demographics alone is sufficient to maintain high diversity. Thus, increasing the number of life-history traits does not necessarily facilitate coexistence because it generates suboptimal strategies. This result differs from other successional models that found that two traits offer distinct mechanisms of successional differentiation13 because, in our model, the effect of different traits can be recast into the same fundamental coexisting mechanism.

Given the disparity between the number of plant traits that can influence \(g\) and \(r\) and the dimension of successional diversity, it is likely that each successional type comprises several ecologically equivalent species, which could explain why coexisting species often exhibit a wide range of functional traits that appear weakly correlated with their position along the successional axis63,64, as shown in the inset of Fig. 1a. Successional strategies not considered in the current model, such as shade-tolerant shrubs that reproduce in the understory or fast-growing species with a small maximum stature, increase the dimensionality of the problem53 but do not preclude equalization.

The addition of CNDD in the successional model does not substantially alter the fundamental reasons for successional coexistence. While it is true that CNDD makes a community more invasible by reducing limiting similarities, it is also true that evolutionary dynamics or repeated random invasions will also reduce limiting similarities to zero in the limit, even in the absence of CNDD. The effect of CNDD in the successional model can be explained as a contraction of the niche shadow in relation to the safety distance (Fig. 5). All else equal, increasing CNDD causes species to become less abundant. This delays competitive exclusion times, as fewer individuals take more time to close the canopy. Longer exclusion times mean potentially more species packing.

In the presence of ecologically equivalent species within each successional type, the effect of CNDD on successional diversity is diluted. When there are two or more species per successional strategy, CNDD promotes diversity within each strategy by reducing ecological drift33. However, in doing so, CNDD becomes less effective in maintaining successional diversity. By definition, CNDD is a function of species abundance and not the abundance of the successional type. Thus, the more ecologically equivalent species within each successional type, the weaker the average strength of CNDD is and, thus, the weaker the effect of CNDD on successional diversity (Supplementary Appendix A). This paradoxical impact on diversity generates the counter-intuitive prediction that hyper-diverse forests are less successionally diverse, which appears to be consistent with the empirical observation that the most diverse tropical forests have the least variation in demographic rates that imply successional diversity, i.e., variation between species-specific growth and mortality rates consistent with successional diversity65,66,67.

Importantly, CNDD modifies the frequency of successional types. The reason for the interaction between CNDD and succession dynamics lies in the power-law function describing \({LRS}\) increase with age. This function is an important regulator of population density. As we demonstrated above, CNDD reduces limiting similarities. However, because the slope of the \({LRS}\) around 1 is steeper for species with shorter \(r\) (Fig. 3b), these species can adjust their density with a smaller contraction of the niche shadow. Thus, early successional species can still maintain relatively high recruit density and large limiting similarity even in the presence of CNDD (Fig. 4b). This mechanism generates a diversity gradient along the successional axis, with more species packed in the late-successional end. Furthermore, CNDD can lead to diminishing benefit returns for species that grow larger in size, especially if its effects were proportional to biomass (in the current model is proportional to area occupied by adults), because biomass grows faster than crown and larger trees have more biomass per unit of reproductive area.

In the successional model, CNDD can substantially increase diversity, but only by an incremental fraction of what is already generated by the resource competition component of the model, which is roughly half of the available species pool. This result is in agreement with27,28, who show that while the introduction of CNDD in a neutral model increases diversity by several folds, it has only marginal effects in tradeoff models. Consequently, in this model, large diversity gradients must be generated primarily by differences in metacommunity size. However, the larger the metacommunity, the larger the increase in diversity promoted by CNDD (Fig. 4a), which might reconcile why CNDD is more important for tropical forests despite evidence that higher latitude forests also experience similar CNDD effects16,19,68.

Besides the fact that our model represents interactions among species and individuals and competitive differences—the lack of such differences was the main limitation of Hubbell’s theory—the two frameworks present some affinities: both aim to predict species diversity and species abundance distribution in a local forest community connected to a larger metacommunity and can be generalized to incorporate CNDD33. Hubbell argued that trait variation in plants may not translate into ecological uniqueness. The collapse of plant diversity into a much smaller space of ecological variation is an interesting parallel between this study and Hubbell’s neutral theory.

Both theories envision a forest mosaic subjected to disturbance69. However, in our model, patch dynamics are explicitly represented, whereas, in Hubbell’s model, the process is simplified to a single step: the patch is immediately filled by an individual randomly selected from the local community or, with some probability, from the metacommunity without accounting for a dynamic structural component or size differences among individuals.

Hubbell’s model operates on a finite domain and maintains the discreteness of individuals, which is appealing for two key reasons: it accounts for drift and allows for direct comparisons with observational data, typically collected on a delimited portion of the forest. This feature also allows computing the species-area curve. In contrast, our model represents a large number of patches and, in its current form, does not simulate drift. However, apart from drift, our theory incorporates three of the four processes influencing species composition and diversity patterns70: selection, dispersal, and speciation, whereas Hubbell’s theory includes drift, dispersal, and speciation.

Finally, we note that, besides patch size and disturbance frequency, which in the current model are scaling factors (all set to one), the model has only four governing parameters: the number of species in the metacommunity, the minimum break-even time, the allometric constant, and the CNDD parameter. However, the model without CNDD is insensitive to the crown-biomass scaling (Eq. (3) and Supplementary Fig. B6), and, in principle, the scaling can be independently estimated from allometric relationships. So, for the case without CNDD, there are only two free parameters:\(\frac{{J}_{M}}{1-{r}_{\min }}\) and a scaling factor to convert the number of individuals into density per patch area. Interestingly, it has the same number of free parameters as the classic neutral model: metacommunity size (or speciation rate) and immigration probability.

The model presented serves as an initial step towards integrating multiple mechanistic processes into a forest diversity theory that is analytically tractable and enables the prediction of macroecological properties. However, before the theory can be tested on data, some of the limitations of the current approach need to be addressed. The theory needs to include other forms of light competition and life-history strategies, such as gap specialists and plants that thrive in the understory, which can generate successional axes that are orthogonal to the axis identified in this model. Furthermore, our model does not track suppressed plants once they are folded in the understory, though this type of model allows us to represent them analytically, at least in a single understory layer71. It could also be possible to include drift and immigration by treating each successional type as a sub-neutral community in a finite area, an avenue promising to reconcile neutral and niche theories.

As mentioned above, the model has only four free parameters. Although the example in Fig. 6c is encouraging, fitting these parameters to data will require integrating demographic and species compositional datasets, including species and size abundance distributions72. For example, though not shown here, this type of model can predict forest structure, which could help constrain some of the model parameters against observed size abundance distributions73.

A comprehensive theory of forest diversity would enhance our ability to predict how ecosystems might respond to environmental shifts, such as climate change or other human alterations of natural systems. This predictive capability is essential for informing conservation efforts and sustainable management practices.