The Equilibrium Theory of Biodiversity Dynamics: a general framework for scaling species richness and community abundance along environmental gradients

Large-scale temporal and spatial biodiversity patterns have traditionally been explained by multitudinous particular factors and a few theories. However, these theories lack sufficient generality and do not address fundamental interrelationships and coupled dynamics between resource availability, community abundance, and species richness. We propose the Equilibrium Theory of Biodiversity Dynamics (ETBD) to address these linkages. According to the theory, equilibrium levels of species richness and community abundance emerge at large spatial scales due to the population size-dependence of speciation and/or extinction rates, modulated by resource availability and the species abundance distribution. In contrast to other theories, ETBD includes the effect of biodiversity on community abundance and thus addresses phenomena such as niche complementarity, facilitation, and ecosystem engineering. It reveals how alternative stable states in both diversity and community abundance emerge from these nonlinear biodiversity effects. The theory predicts how the strength of these effects alters scaling relationships between species richness, (meta)community abundance, and resource availability along different environmental gradients. Using data on global-scale variation in tree species richness, we show how the general framework is useful for clarifying the role of speciation, extinction and resource availability in driving macroecological patterns in biodiversity and community abundance, such as the latitudinal diversity gradient.

and alternative stable states in both diversity and community abundance emerge from 23 nonlinear BEF relationships and the dependence of origination and extinction on population 24 size. ETBD predicts how the strength of the BEF affects scaling relationships between species 25 richness and (meta)community abundance along different environmental gradients. Using data 26 on global-scale variation in tree species richness, we show how the general framework is useful 27 for clarifying the role of speciation, extinction, resource availability, and temperature in driving 28 biodiversity patterns such as the latitudinal diversity gradient. 29 Introduction 33 The biodiversity, abundance and biomass of living systems vary extraordinarily across the 34 biosphere, through Earth history, and among major phyla. Understanding this large-scale 35 variation and its consequences for the functioning and resilience of ecological systems is crucial 36 for advancing ecological theory and forecasting the future of biodiversity and ecosystems in the 37 Anthropocene (1, 2). Given the elevated extinctions underway in the Anthropocene, identifying 38 the causes of global-scale patterns of species richness, community abundance and biomass is of 39 particular importance. 40 The underlying causes for even the most well-studied patterns still remain hotly 41 contested (3). The pattern of decrease in the number of species from the equator to the poles, 42 documented in nearly all major taxa of multicellular eukaryotes, is one prominent example. 43 Dozens of hypotheses have been formulated (4), but no clear consensus has emerged. This lack 44 of consensus is in part due to the multitudinous factors at play but also due to the lack of 45 general theoretical frameworks to test hypotheses and disentangle large-scale diversity drivers. 46 Recently, evidence is accumulating that large-scale spatial diversity patterns tend to converge 47 onto similar relationships with resource availability and climatic variables regardless of 48 differences in diversification histories (5), and the origination and extinction rates underlying 49 diversity dynamics exhibit a dependence on diversity (6, 7). These findings have been 50 interpreted as evidence of the role of region or biome-specific ecological limits to species 51 richness (8, 9), but the nature of these limits has not been entirely clear. We argue that such 52 limits should be understood as stable equilibria of diversity dynamics that shape large-scale 53 spatial and temporal patterns of diversity (10). To understand how biodiversity equilibria 54 emerge and depend on the properties of environments and biota requires addressing the links 55 between macroevolutionary dynamics and macroecological processes (11). 56 General theories of biodiversity dynamics, such as the theory of island biogeography 57 (12), neutral theory of biodiversity (13) and their derivatives (14), assume that biodiversity 58 equilibria exist but suffer from some key shortcomings preventing them from providing a 59 general framework to understanding diversity dynamics. First, they make overly stringent and 60 unrealistic assumptions about the involved processes. For example, the neutral theory of 61 biodiversity, besides assuming ecological equivalence among species that includes equal access 62 to all resources, implicitly assumes that per-species speciation rate increases linearly with 63 population size (15) despite theoretical evidence that speciation rate can exhibit more complex 64 dependencies on species abundance (16). Second, most crucially, while these theories address 65 (implicitly or explicitly) how community abundance (total biomass or number of individuals) 66 limits species richness by influencing rates of extinction and origination (e.g., 13, 15), they take 67 community abundance as an external parameter or static variable and do not address the 68 possibility that the size of the community may itself be a dynamical variable that depends both 69 on the environment and the realized species richness of the community. 70 Here we argue that the number of species in a community tends to positively affect the 71 total amount of resources used and converted into community abundance and biomass (17-72 20). The major reasons include (1) that species vary in their resource use and conversion 73 efficiencies, so that the addition of new species may increase the conversion of resource into 74 biomass or lead to exploitation of new resources, and (2) ecosystem engineering and 75 facilitation between species can augment community abundance (17,(21)(22)(23). A directional link 76 between species diversity and community abundance has been well-documented in 77 experimental studies of the biodiversity-ecosystem functioning (BEF) relationship (20,24,25), 78 suggested by macroecological and paleobiological analyses (e.g., 26,27), and demonstrated 79 with theoretical models (17, 28). There is thus a need for a theory that addresses the coupled 80 dynamics of species richness and community abundance. 81 Such coupled dynamics may have at least two non-trivial consequences for biodiversity. 82 First, the bidirectional interactions between biodiversity and community abundance could alter 83 equilibria along environmental gradients. Whereas the "More-Individuals Hypothesis" (29), 84 neutral theory of biodiversity (13), and metabolic theory of ecology (30) all predict linear and 85 unidirectional effects of community abundance and speciation rate on species richness , the bi- 86 directional relationships between diversity and community abundance may lead to sublinear 87 scaling (a power law with an exponent less than one) or superlinear scaling (exponent greater 88 than one) between energy availability, community abundance and species richness even if an 89 environmental gradient reflects the variation of only one abiotic or biotic driver of biodiversity 90 variation (e.g., only the variation of energy or resource availability). Moreover, if more than one 91 driver varies along a gradient (e.g. both energy availability and speciation rate), predicted 92 scaling patterns may be complex and depend on their covariation along the gradient. 93 Second, in dynamical systems, multiple stable states and tipping points can emerge 94 from positive feedback loops (31), such as the one between biodiversity and community 95 abundance resulting from the BEF relationship combined with the baseline positive effect of 96 community abundance on diversity. Multiple stable states and regime shifts could thus shape 97 temporal and geographic patterns of biodiversity and biomass, as well as their response to 98 environmental perturbations. Ecosystem and community ecologists have established that 99 feedbacks involving specific combinations of physical, chemical, and biological mechanisms can 100 cause tipping points and multiple stable states in specific biomes (32), ecosystems (33), 101 communities (34), and food webs (35), but the degree to which ecological tipping points can 102 operate up to scales of regional diversity or the biosphere is contested (36-38). Furthermore, The central assumption of ETBD is that a population's probability of extinction and/or 123 speciation exhibit some level of dependence on population size. There is a solid empirical and 124 theoretical basis for this assumption, in particular that extinction probability is generally 125 negatively dependent on population size (e.g., 42-45), especially at small to intermediate 126 population sizes, due to increased population size conferring increased resistance to 127 demographic and environmental stochasticity and the detrimental effects of genetic drift in 128 small populations (41, 46). A consequence of this assumption is the emergence of equilibria of 129 species richness and community abundance that depend on extinction, speciation, total 130 resource availability, and community resource utilization, which is modulated by the outcome 131 of competitive, neutral, and facilitative interspecific interactions (Fig. 1). 132 Before presenting ETBD's mathematical machinery and quantitative predictions, we 133 begin by illustrating the bones of ETBD with the scenario of biodiversity dynamics that we 134 expect universally applies to a variety of taxa and environments (Fig. 1). For communities that 135 aren't on their way to extinction, the curve quantifying the population-size (N) dependence of 136 extinction must necessarily intersect the curve quantifying the N-dependence of speciation. The 137 reason is that the origination rate cannot be higher than the extinction rate over all values of N, 138 since as N approaches 0, extinction probability necessarily approaches 1 whereas origination 139 probability approaches zero. Consequently, at this universal intersection, the (negative) slope of 140 the extinction curves must always be lower than the slope of the speciation curve ( Figure 1A). 141 This intersection of origination and extinction curves indicates a population size at which 142 origination and extinction rates are equal. This implies that in a plot of species richness S 143 against community abundance J on the vertical axis, there is a straight diagonal line called the 144 S-nullcline that has a slope equal to mean equilibrium population size ̂ (since the slope is 145 equal to J/S, which is the mean population size) and represents communities in which total 146 5 origination rate is equal to total extinction rate (i.e., dS/dt = 0; Fig 1B,  Materials for mathematical details). The equation for this S-nullcline is thus =̂. Because 148 the negative slope of the extinction curve is lower than the slope of the origination curve, this 149 nullcline exhibits sink dynamics in which communities perturbed out of equilibrium sink back to 150 the S-nullcline, which we thus call a sink S-nullcline. The reason is that points to the right of this 151 S-nullcline represent communities where mean population size (i.e., J/S) is lower than ̂ (and 152 thus extinction is generally higher than origination rate), and vice versa. All points on this 153 universal sink S-nullcline thus represent possible diversity equilibria. Origination (green) and extinction (red) rates depend on population size in various ways, but around their intersection point the extinction curve has a more negative slope than the origination rate. This intersection indicates the existence of a population size at which origination and extinction rates are equal. (B) In the phase plane for species richness S and community abundance J, this implies a diagonal line (S-nullcline, blue line) with a slope determined by equilibrium mean population size ̂. This S-nullcline indicates values of S and J in which total origination rate is equal to total extinction rate. Points directly to the right of this line represent communities where J/S < ̂; mean population size is thus smaller than ̂, leading to extinction rate being higher than origination rate and the system having the tendency to move back towards the S-nullcline. The opposite happens if J/S > ̂ (left from the S-nullcline). Equilibrium S and J occur at the intersection (solid red point) of the S-nullcline with the J-nullcline (thick red line), which quantifies the effects of biodiversity on community abundance via the biodiversity-ecosystem function relationship and demarcates the community abundance at which additions and losses to community abundance are in balance. The J-nullcline may take a variety of forms, but eventually plateaus, being constrained by energy availability (thin red line). Perturbed communities below the J-nullcline follow upward trajectories converging on the J-nullcline, whereas out-of-equilibrium communities above the J-nullcline follow downward converging trajectories. Once on the J-nullcline, trajectories flow along the J-nullcline in the direction determined by their position with respect to the S-nullcline. (C) The position of equilibria can change (red arrows) with changing (1) energy availability (top), (2) origination and/or extinction rates that affect mean equilibrium abundance and thus the S-nullcline (middle), or (3) the shape of the J-nullcline (bottom). If the J-nullcline has a more complex shape, there can be multiple equilibria, both stable and unstable (empty circle). Unstable equilibria on this S-nullcline emerge when this S-nullcline is crossed by a locally upward-accelerating J-nullcline.
The exact equilibrium position of a community along the S-nullcline is then determined 156 by the community's J, which may itself be influenced by any effects that S has on J. Instead of 157 assuming constant J, ETBD assumes that J is driven by the interplay of population dynamics of 158 all involved species and that communities with higher numbers of species better utilizes 159 resources. The J-nullcline, the line delineating the values of S and J for which J can be in 160 equilibrium (i.e., for which dJ/dt = 0), is thus an increasing function of S (Fig 1B thick red line), 161 but in a decelerating manner, as the absolute maximum J is determined by the maximum 162 energy available to the system divided by mean individual metabolic rate (thin red line). A 163 community thus has an equilibrium S and J (̂ and ̂) that is located at the intersection of the S-164 nullcline and J-nullcline. Consequently, three fundamental drivers of the system can alter the 165 position of the equilibrium (Fig 1C). increases in the efficiency by which a resource is converted into biomass by more species (e.g., 200 due to specialists having higher resource use efficiencies or "portfolio effects'' that buffer 201 perturbations)(47). 202 In Figure 1C

293
Box 1 Figure 1. Example of how a sigmoidal biodiversity-ecosystem function relationship (J-nullcline) in conjunction with nonlinearities in the origination or extinction curve lead to multiple stable states of biodiversity and community abundance equilibria. (A) When the speciation curve drops at higher N or the extinction curve increases, there can be two intersections of the curves. The mean population sizes (̂1 and ̂2 ) for which a community can be at equilibrium is largely determined by these intersections. (B) These equilibrium mean population sizes determine the heights of the resulting S-nullclines in the log-log S-J phase plane (black lines). The slope of the extinction curve relative to the slope of the origination curve at their intersection determines the nature of the associated S-nullcline shown and the phase space in which dS/dt < 0 (blue regions) versus dS/dt > 0 (pink region). Communities recovering from perturbations converge to the J-nullcline (purple line) and then follow trajectories indicated by the purple arrows. Here a sigmoidal J-nullcline is shown, leading to four equilibrium points. In A, the origination curve necessarily drops to zero at non-reproductively viable population sizes, whereas the extinction curve intercepts 1.
Scaling relationships for gradients in energy, speciation, and extinction 294 ETBD highlights that there are three fundamental drivers of biodiversity equilibria: energy 295 availability, the extinction-speciation balance, and the shape of the J-nullcline. We thus first 296 address how equilibrium S and J should change as a consequence of each of these three 297 fundamental drivers, and then present quantitative predictions for complex gradients in which 298 more than one of these drivers are varying simultaneously. 299 The ETDB framework reveals that the scaling behavior of the J-nullcline, as quantified by 300 , fundamentally affects the scaling relationships between S, J, and E/B. Over a given range of 301 S, the J-nullcline can be mathematized as a power-law function of S (e.g., 27), giving which only the position of the J-nullcline is changing: 309  330 where 0 and 0 are normalization coefficients and 1 and 1 sare scaling exponents quantify 331 the effects of 1 . Along a purely resource-driven gradient, equilibrium mean abundance 332 doesn't vary and so 1 = 0. Along a purely speciation/extinction-driven gradient, resource 333 availability doesn't vary and so 1 = 0. When 1 1 < 0, the two drivers (the positive effects of the 334 speciation/extinction balance and resource availability on diversity) are concordant along the 335 gradient, leading to what we call a concordance gradient. When 1 1 > 0, the effects are 336 discordant such that at higher resource availability, speciation is lower or extinction is higher, 337 leading to a discordance gradient. In a perfectly discordant gradient, the drivers of diversity 338 vary completely opposingly, leading to 1 1 = 1 and no change in diversity along the gradient. 339 The resulting scaling predictions for these gradients are:  The curves indicate the S-J-E scaling exponents for pure speciation/extinction gradients (light blue, r1 = 0), pure resource availability gradients (black, g1 = 0), and complex gradients in which speciation/extinction and resource availability covary. In a complex gradient, equilibrium mean population abundance ̂ may decrease with resource availability due to higher speciation or lower extinction in high resource environments (a concordance gradient, pink curves, g1/r1 < 0) or increase with resources (a discordance gradient; dark blues curves, g1/r1 > 0). Thicker curves indicate gradients having larger absolute values of g1/r1 (as shown). Note that ETBD predicts all communities have an upper stable diversity and abundance equilibrium resulting from competition for limited resources ( < 1); but depending on the form of the extinction and speciation curves, they may also have a lower stable diversity and abundance equilibrium in which facilitation ( > 1) dominates the BEF relationship in lower diversity communities.
superlinear relationship (with < − 1 / 1 ) when there is relatively strong competition but 361 still some niche complementarity (specifically, < 1 / 1 ); (2) a positive superlinear 362 relationship when there is relatively weak competition and substantial niche complementarity 363 (specifically, 1 / 1 < < 1); and (iii) a positive sublinear relationship (0 < <1) when 364 facilitation dominates. When resource availability instead varies more than the 365 speciation/extinction balance along the gradient, then 0 < < 1 when there is competition 366 ( < 1), whereas > 1 when facilitation dominates. In a perfectly discordant gradient, no 367 change in biodiversity along the gradient occurs but there is a change in J along the gradient  Figure 3 shows that the shape of the SAD, whether a 469 same-abundance distribution or lognormal distribution with low or high variance, has minimal 470 impact on overall changes in biodiversity along gradients in speciation/extinction or resource 471 availability, although it certainly modulates equilibrium mean abundance and so it could be 472 important to consider if the shape of the SAD is suspected to change drastically over a gradient.

474
Empirical application of the framework 475 We applied the scaling framework to variation in tree species richness across the globe in order 476 to clarify drivers of the latitudinal diversity gradient, using a multi-scale tree biodiversity 477 database (60)(Materials and Methods). Some authors argue that energy/resource availability is 478 a major driver of latitudinal diversity gradients (e.g., 61-64), whereas others have argued that 479 variation in speciation or extinction rate are a major driver (e.g., 65-68). Since ETBD makes representing large scales, high-quality abundance data for such large scales is rarely available. 494 We therefore first checked whether the average local number of species in plots of a given area 495 scales predictably with the biome-average regional-scale (100,000 km 2 ) richness, in order to biogeographic realms (Fig. 4C). The observed average S-J scaling is thus ∝ 3.89 , which is 517 remarkably close to the predicted 3.85 for a purely speciation/extinction driven gradient 518 Figure 4. Predicted and observed patterns of average tree richness of biomes from different biogeographic realms. The local richness down to 0.01 km 2 is proportional to regional richness (A), resulting in parallel latitudinal diversity gradients in (B). In C, D, and E, the blue lines show the predicted S-J scaling relationships and the relationships between S or J and temperature assuming a purely speciation/extinction driven gradient (assuming r1 = 0) and = 0.26. Predictions are indistinguishable from the regressions (black lines). Red lines in c represents the contrasting scenario of a purely resource-driven gradient (i.e. slope of 1; g1 = 0). Blue lines in D and E must also assume a value for g1 (quantifying the degree to which speciation/extinction changes along the gradient), which is estimated from the scaling of ln N with 1/kT. Interpretation of the observed patterns within the ETBD framework suggests that in trees the balance of speciation and extinction, rather than variation in energy availability, is the major driver of global-scale biodiversity variation, including the latitudinal gradient.
(where 1 = 0 and 1 ≠ 0), whereas the 95% confidence intervals of the observed scalings do 519 not overlap with one as expected for a purely E/B driven gradient (where 1 ≠ 0 and 1 = 0). In 520 other words, the patterns found are in accord with the idea that tree diversity variation across 521 biomes (and thus the latitudinal diversity gradient) is driven by the variation of speciation and 522 extinction rates, and not resource availability. 523 We can now use ETBD to predict how species richness should change with temperature 524 across these biomes. Since biological rates generally scale with temperature as   Table 1. Results of the Reduced Major Axis regression analysis of log species richness (S), log total number of individuals (J), and log mean population abundance (̂) in trees from terrestrial biomes and the comparison of observed slopes to ETDB-predicted slopes for a speciation/extinction gradient. to have greater species richness than the non-equilibrium communities with lower stable 617 equilibria, both because these communities have a greater "distance" to fall from and because 618 the further a community is from a stable equilibrium point, the greater its rate of increase in 619 diversity, leading to the high-equilibria community bouncing back to a higher richness more 620 quickly. While these points may seem trivial, they can help resolve longstanding debates 621 regarding the role of "cradles", "museums", and "graves" in shaping global-scale biodiversity 622 patterns (90). 623 Application of our framework to global-scale variation in tree biodiversity suggests that 624 tree diversity variation is driven by variation in the speciation and extinction curves whose  variable. For each biome-realm combination, we then calculated the mean of log regional 788 species richness (Data S1). 789 Plot-level data. We binned the plot data according to plot area, using bin widths of one 790 order of magnitude. In practice most of the bins had much lower variation in plot sizes, as the 791 bins were mostly made up of plots of similar standard areas (and minimum DBH). The bin for 792 plots of 0.01 to less than 0.1 km 2 was almost entirely made up of plots of 0.01 km 2 , so we only 793 included plots of 0.01 km 2 for this bin, in order reduce confounding error. The bin for plots of 794 0.1 to less than 1 km 2 had plot areas varying from 0.14 km 2 to 0.5 km 2 with the majority having 795 sizes of 0.2-0.25 km 2 . For each plot area bin, we calculated the mean of the following variables 796 for each biome-realm combination: log species richness, log number of individuals, mean 797 annual air temperature, and mean absolute latitude (Data S1).  119). To evaluate the sensitivity of our results to our method for calculating regional 804 species richness, we also repeated all regional species richness analyses (local versus regional 805 richness and regional richness versus absolute latitude) using raw country/state richness data 806 without implementing the normalization for variation in country/state area. In other words, for 807 each biome-realm combination, we calculated the mean log country/state species richness and 808 repeated the previously described local-regional scaling and gradient analyses on this measure 809 of regional richness. The results from this analysis were similar to our other findings (Table S1). , then source dynamics are exhibited. 1204 We have not uncovered any cases in which typical SADs (logseries, lognormal, exponential, 1205 power law) add any additional S-nullclines and thus it seems that the properties of the S-1206 nullclines generally follow the described rules, as long as (̇) isn't very close to (̇) .

1207
Consequently, we only consider in the main text the general expectations that there are either just 1208 one sink S-nullclines ( Figure 1) or two S-nullclines when is hump-shaped, leading to two ways, this diversity dependence could readily be incorporated into Eq. S1 (e.g., by making ( ) 1221 a power law function of S, which could affect the stability of the equilibrium).

1222
To help better intuit these implications, consider the population dynamics and resulting 1223 community structure acting within a time interval. It follows from the above assumptions that at  Table S1. Results of the Reduced Major Axis regression analysis using country-level species 1398 richness as a measure of regional species richness without standardizing to 100,000 sq. km.