Thompson et al. - 2018 - Biodiversity-Ecosystem Function Relationship PDF

The strength of the biodiversity – rspb.royalsocietypublishing.org ecosystem function relationship...

The strength of the biodiversity – rspb.royalsocietypublishing.org ecosystem function relationship depends on spatial scale Patrick L. Thompson1,2, Forest Isbell3, Michel Loreau4, Mary I. O’Connor1,2 Research and Andrew Gonzalez5 Cite this article: Thompson PL, Isbell F, 1 Department of Zoology, University of British Columbia, and 2Biodiversity Research Centre, Loreau M, O’Connor MI, Gonzalez A. 2018 University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z4 3 The strength of the biodiversity – ecosystem Department of Ecology, Evolution and Behavior, University of Minnesota, Saint Paul, MN 55108, USA 4 Centre for Biodiversity Theory and Modelling, Theoretical and Experimental Ecology Station, CNRS and function relationship depends on spatial scale. Paul Sabatier University, 09200 Moulis, France Proc. R. Soc. B 285: 20180038. 5 Department of Biology, McGill University, Montreal, Quebec, Canada H3A 1B1 http://dx.doi.org/10.1098/rspb.2018.0038 Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 PLT, 0000-0002-5278-9045 Our understanding of the relationship between biodiversity and ecosystem functioning (BEF) applies mainly to fine spatial scales. New research is Received: 5 January 2018 required if we are to extend this knowledge to broader spatial scales that are Accepted: 11 May 2018 relevant for conservation decisions. Here, we use simulations to examine con- ditions that generate scale dependence of the BEF relationship. We study scale by assessing how the BEF relationship (slope and R 2) changes when habitat patches are spatially aggregated. We find three ways for the BEF relationship to be scale-dependent: (i) variation among local patches in local (a) diversity, Subject Category: (ii) spatial variation in the local BEF relationship and (iii) incomplete compo- Ecology sitional turnover in species composition among patches. The first two cause the slope of the BEF relationship to increase moderately with spatial scale, Subject Areas: reflecting nonlinear averaging of spatial variation in diversity or the BEF ecology, theoretical biology relationship. The third mechanism results in much stronger scale dependence, with the BEF relationship increasing in the rising portion of the species area Keywords: relationship, but then decreasing as it saturates. An analysis of data from the species richness, ecosystem functioning, Cedar Creek grassland BEF experiment revealed a positive but saturating slope of the relationship with scale. Overall, our findings suggest that the b-diversity, Jensen’s inequality, BEF relationship is likely to be scale dependent. nonlinear averaging, spatial scale 1. Introduction Author for correspondence: Much of our understanding of how biodiversity affects ecosystem functioning stems Patrick L. Thompson from hundreds of experimental studies and field observations conducted at rela- e-mail: [email protected] tively small scales of space and time, often considered ‘local scales’ [1–3]. These experiments consider biodiversity effects to be local, because this is the scale at which species interact and compete for resources (e.g. less than 200 m2). At the local scale, selection and complementarity effects generally cause ecosystem functioning to increase with species richness in a positive but decelerating fashion [1,5–7]. Do these findings, and the theory that explains them, apply to larger scales where resource complementarity may occur over regional environmental gra- dients [8,9]? This gap in our knowledge limits our ability to link biodiversity and ecosystem functioning (BEF) science to the larger spatial scales where conservation decisions are generally made, and which are relevant for the provisioning of many ecosystem services. Extending local mechanisms to broader scales is possible; recent advances in theory now allow us to scale up our understanding of the bio- diversity–ecosystem stability relationship [11,12]. Here, we address the similar challenge of providing scaling theory for the BEF relationship. Electronic supplementary material is available Scaling the BEF relationship to larger spatial extents is more complicated than simply applying the local BEF relationship to the greater number of species that are online at https://dx.doi.org/10.6084/m9. present in larger regions. Spatial variation in environmental conditions, con- figshare.c.4105337. nectivity, biotic interactions and stochastic processes causes local communities & 2018 The Author(s) Published by the Royal Society. All rights reserved.  to differ in composition and in the shape of their BEF this up to larger regions, for example, landscapes comprised 2 relationships [2,14]. Recent studies have highlighted the impor- many local habitat patches, requires that we consider the rspb.royalsocietypublishing.org tance of compositional turnover in space (b-diversity) in relationship between g-diversity, the diversity present across maintaining high rates of ecosystem functioning at landscape all patches at a given spatial scale, and the total ecosystem func- scales, because greater diversity is required to maintain ecosys- tioning at that scale. For b to remain constant across spatial tem functioning across the range of environmental conditions scales, then proportional changes in a- and g-diversities must present in larger regions [15–19]. Because of this, we might result in the same proportional change in ecosystem function- expect that the BEF relationship should become stronger at ing. Our question is therefore: what causes proportional larger spatial scales [8,9]. This importance of regional biodiver- changes in a- and g-diversities to result in different proportional sity is further supported by metacommunity models changes in ecosystem functioning, and when this occurs, does it comprising two spatial scales, where productivity is more result in an increase or decrease in the BEF slope? Furthermore, strongly dependent on regional than local diversity, but the how does the predictive power of the BEF relationship change Proc. R. Soc. B 285: 20180038 slope of the BEF relationship locally and regionally is mediated with spatial scale? And do effects of scale depend on whether by dispersal, which affects the b-diversity present in the region the BEF relationship is strong, weak, positive or negative? [20–22]. However, Cardinale et al. demonstrated that the Here we develop basic theoretical expectations for how the slope of the BEF relationship can be constant across spatial slope of the biodiversity2ecosystem functioning relationship scales (from small to large extents) if the local BEF relationship might change with spatial scale. For this, we use simulation Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 as well as the community composition is equal across all local models to explore how variation in different BEF parameters patches. Thus, in this unlikely scenario of homogeneous con- (a, S and b, of Y ¼ aSb), as well as different patterns of b-diversity, ditions across space, changes in local diversity result in causes the BEF relationship to change with spatial scale. In similar proportional changes to the biodiversity of the region. taking this approach, we do not directly consider the ecological However, Cardinale et al. and Chesson et al. found mechanisms driving this variation (e.g. environmental hetero- that the processes driving the BEF relationship change with geneity, connectivity and biotic interactions), which can also spatial scale, with selection effects becoming complementarity directly affect rates of ecosystem functioning. Rather, we effects at broader spatial scales if species tend to dominate simply ask how spatial variation in community composition local habitat patches that match their niche requirements. and the shape of the local BEF relationship might be expected To develop theoretical expectations for scaling the BEF to cause the BEF relationship to change with spatial scale. We relationship, we can start with our understanding of the do this by simulating regions composed of local habitat patches relationship at local scales for which there exists substantial where functioning in each local patch depends on the number theoretical understanding and empirical evidence. Theory of species in that patch, following Y ¼ aSb. We aggregate suggests that the relationship between ecosystem functioning these local patches to estimate the BEF relationship over and local species richness (a-diversity) is driven by comple- larger spatial scales. Regional ecosystem functioning is mentarity and selection effects [5,7,25]. This theory is assumed to be the sum of functioning in all patches, while supported by empirical evidence from hundreds of exper- g-diversity is assumed to be the number of species across all iments that report a relationship that can be described as a local patches in the region. Therefore, g-diversity will be less power law, Y ¼ aSb, where Y is the level of ecosystem function, than the sum of the a-diversity in all patches when there is com- S is the number of species present, a is a constant, and b is the positional overlap among patches. We start with a control case exponent of the BEF relationship (or the slope in log–log study (case I), in which we demonstrate how the BEF relation- space), which indicates the strength of biodiversity effects (i.e. ship remains constant across spatial scales when (i) the BEF the proportional change in ecosystem functioning per change relationship is constant across all local patches and (ii) commu- in richness [1,2,6,26]). The shape of this relationship, described nity composition is equal across all local patches, following by the value of b, can vary among locations and communities, Cardinale et al.. We then explore four case studies, each and with attributes, such as trophic level and over time. deviating from these initial criteria in one way: case II, variation In reality, biodiversity generally only explains a fraction of in local a; case III, variation in local S; case IV, variation in local the variation in the ecosystem functioning at local scales b; case V, b-diversity among local patches. In cases II–V, we because other factors, such as soils or climate, are also impor- explore how the slope and explanatory power of the BEF tant for functioning. For example, in experiments relationship changes with spatial scale. Cases III–V, but not I designed to test for the BEF relationship, it is typical for less or II, change the slope of the BEF relationship with spatial than half of the variance in ecosystem functioning to be scale. Of these, case V (b-diversity) shows the largest scale explained by local diversity (e.g. ). If environmental factors dependence, whereas the changes in b with increasing spatial change with spatial scale, then we may expect the explanatory scale in cases III and IV are relatively minor in comparison. power of the BEF relationship (i.e. the coefficient of determi- We then compare these theoretical expectations with empirical nation, R 2) to also change with spatial scale. In scaling up, it data from the Cedar Creek biodiversity experiment to show is unclear whether biodiversity will still be an important pre- how our expectations hold in empirical communities where dictor of ecosystem functioning. Therefore, in considering the local BEF relationship arises through biological interactions. how the BEF relationship changes with spatial scale, we must consider both changes in the slope of the relationship as well as changes in the explanatory power of the relationship. This power-law relationship Y ¼ aSb allows us to consider 2. Simulation model what would be required for the BEF slope, b, to change with We simulate the BEF relationship at multiple spatial scales by spatial scale. In previous applications of this relationship, S is modelling regions composed of local habitat patches governed by equal to a-diversity, the diversity present at the local scale where species interact and compete for resources. Scaling Yi ¼ ai Sbi i , ð2:1Þ  (a) b = 1.25 (b) b = 0.75 3 rspb.royalsocietypublishing.org 60 10 ecosystem functioning 40 b=1 b = 0.5 5 20 Proc. R. Soc. B 285: 20180038 b = 0.25 b=0 0 0 b = –0.25 0 10 20 30 0 10 20 30 species richness species richness Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 Figure 1. Illustration of the shape of the BEF relationship with different values of b. Higher values of b are shown in panel (a), lower values of b are shown in panel (b). This separation was done for clarity to show the curvature of the BEF relationship for low values of b. where Yi is the level of ecosystem function in patch i, Si is the rounded to the nearest integer and, excluding non-positive number of species present, ai is a coefficient that determines the cases. magnitude of the BEF relationship in a patch, and bi is We find that bA and the R 2 of the BEF relationship remain the exponent of the BEF relationship in that patch. We scale constant across spatial scales (figure 2; electronic supplemen- ecosystem functioning to larger spatial areas by tary material, figure S1). This is true regardless of whether we assume that all patches contain the same set of species or X A YA ¼ Yi , ð2:2Þ completely unique species (same relationship as figure 2). i¼1 where A is the number of local patches comprising the region. SA is the number of species in the region. At each spatial scale, A, we 4. Case II: variation in ai across local habitat estimate the slope of the linearized BEF relationship, bA , given patches YA ¼ aA SbAA , ð2:3Þ In case II, we draw values of ai from a normal distribution which we estimate using linear, least-squares, regression on a with a mean of 5 and a standard deviation of 2, excluding log–log scale across 2000 simulated replicate regions. We non-positive values. We assume that bi is equal across all repeat this procedure 100 times to obtain replicate estimates of patches, and that all patches contain the same species. bA to estimate variability across simulation runs. To test the sen- We find that variation in ai does not cause the mean BEF sitivity of results to different possible BEF relationships, we slope, bA, to change with spatial scale (figure 2a). This is contrast seven different values of bi ranging from 20.25 to because a is a linear coefficient in the local BEF relationship 1.25. This incorporates the various shapes of the BEF curve and so there is no potential for nonlinear averaging to cause observed in empirical data [1,2,28]: negative decelerating the BEF relationship to change with spatial scale. However, ð1 , bi , 0Þ, no relationship ðbi ¼ 0Þ, positive decelerating variation in ai does reduce the explanatory power of the ð0 , bi , 1Þ, positive linear ðbi ¼ 1Þ and positive accelerating BEF relationship (figure 2b), but this explanatory power ðbi. 1Þ (figure 1). In our simulations, we have assumed that increases as spatial scale is increased. Aggregating patches at the BEF relationship follows a power law. However, our larger spatial scales reduces noise in the local BEF relationship conclusions should not depend on the specific functional associated with variation in ai. form used, but rather whether the BEF relationship is linear or nonlinear, and whether the curve is concave up or down. All simulations were performed in R v. 3.4.2. 5. Case III: variation in a richness across local habitat patches In case III, we draw values of Si from a normal distribution with 3. Case I: control—no variation in ai, a mean of 10 species and a standard deviation of 3, rounded to a richness, bi, and no (or maximum) the nearest integer and excluding non-positive values. We assume that bi is equal across all patches, and that all patches b-diversity across patches contain unique species, to avoid incomplete compositional turn- In case I, we assume that all patches in a region have equal ai, Si over. Note that it is not possible to have complete compositional and bi, and that b-diversity is either zero (all patches with the overlap between patches when Si varies across patches. same composition) or maximal (all patches contain different We find that the mean BEF slope, bA, changes with spatial species). We generate variation in diversity across replicate scale, and this variation across scales depends on the shape of regions by drawing different values of Si from a normal distri- the local BEF relationship (figure 2a). When the local relation- bution centered on 10 with a standard deviation of 3, ship is positive and decelerating ð0 , bi , 1Þ, as is most often  (a) case I case II case III case IV 4 rspb.royalsocietypublishing.org relative BEF slope (bA − bi) 0.02 bi −0.25 0 0 0.25 0.5 0.75 −0.02 1 1.25 Proc. R. Soc. B 285: 20180038 −0.04 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 scale (b) case I case II case III case IV 1.00 Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 R2 of BEF relationship 0.75 0.50 0.25 0 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 scale Figure 2. Illustration of how the BEF relationship changes with spatial scale in cases I– IV. Panel (a) shows how the BEF slope, bA changes relative to the mean local slope bi , at different spatial scales (see electronic supplementary material, figure S1 for raw values of bA). Panel (b) shows the R 2 of this relationship. The four cases are shown in the different panels—no variation in local species richness or local bi (case I), variation in local ai (case II), variation in local S (case III) or variation in local bi (case IV). The solid line shows the median across 100 replicate simulations each consisting of 2000 replicate regions at each scale. Inter-run variability omitted for clarity in panel (a), but shown in electronic supplementary material, figure S1. In panel (b), the interquartile range is smaller than the width of the lines and so is not shown. (Online version in colour.) the case in empirical data [2,28], bA increases as a saturating when the BEF relationship is accelerating (bi. 1). Conse- function of increasing spatial scale. When the local relationship quently, whether bA increases or decreases with spatial scale is linear, either due to no relationship ðbi ¼ 0Þ or a linear depends on local BEF relationships’ shape. In summary, relationship ðbi ¼ 1Þ between species richness and function, when the relationship is concave down (e.g. positive decelerat- bA does not change with spatial scale (electronic supplemen- ing; figure 1), as is typical of local BEF relationships, bA tary material, figure S1). When the local relationship is increases with spatial scale. When the relationship is concave negative and decelerating ð1 , bi , 0Þ or positive and accel- up (e.g. positive accelerating or negative decelerating), bA erating ðbi. 1Þ, bA decreases as a saturating functioning of decreases with spatial scale. increasing spatial scale. Across multiple local patches, the strength of biodiversity effects changes with spatial scale to a greater or lesser extent, as shown by the interquartile range 6. Case IV: variation in bi across local habitat around the mean trends in (electronic supplementary material, figure S1), due to variation in local species richness across repli- patches cate random draws of species richness within patches. Across In case IV, we hold Si equal across all patches in a region but all possible values of bi , except when the local BEF relationship now draw bi for each patch from a normal distribution is linear (bi ¼ 0 or 1Þ, there is a slight decrease in the R 2 of centred on bi , with a standard deviation of 0.15; this variance the BEF relationship with scale, with the greatest decrease corresponds roughly to what has been observed in exper- occurring with the smallest values of bi (figure 2b). imental grassland plant communities [2,14] and in local Variation in a richness across a region leads to a saturating forest communities. We generate variation in diversity change in bA with increasing spatial scale because of nonlinear across replicate regions as in case I. To avoid incomplete com- averaging (i.e. Jensen’s inequality [30,31]) of the contribution of positional turnover, we assume that all patches contain each local patch to the functioning of the region (electronic sup- unique species, but our estimated values of bA are the same plementary material, figure S2a). That is, because the local BEF if we assume that all patches share the same set of species, relationship is nonlinear, changes in diversity result in greater as this would also meet this requirement. changes in ecosystem functioning in low-diversity compared We find that the mean BEF slope, bA, increases as a saturat- with high-diversity patches, although the opposite is true ing function of increasing spatial scale, regardless of the value  (a) (b) 5 1.00 rspb.royalsocietypublishing.org 0.4 0.75 R2 of BEF relationship BEF slope (bA) 0.3 0.50 B1 −0.75 −0.5 −0.25 0.2 0.25 −0.1 Proc. R. Soc. B 285: 20180038 −0.05 −0.025 −0.01 0 0 0 10 20 30 40 50 0 10 20 30 40 50 scale scale Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 Figure 3. The strength of biodiversity effects, bA (panel (a)) and R 2 (panel (b)), at different spatial scales when there is incomplete compositional turnover across local patches. Different degrees of compositional turnover are indicated by colour (low values of B1 correspond to low turnover, B1 ¼ 0 indicates complete turnover). The solid line indicates the median across 100 replicate simulations each consisting of 2000 replicate regions at each scale. The bands show the interquartile range. (Online version in colour.) of bi (figure 2a). In this case, there is considerable variation in bi communities that have been combined. In our simulations, we across local patches (electronic supplementary material, figure set b0 ¼ 5, and we contrasted a range of b1 (20, 20.05, 20.1, S1), as dictated by our assumptions, but because this variation 20.25, 20.5 and 20.75) to explore how the BEF scaling is averaged at larger spatial scales, we find that the dispropor- relationship depends on the species accumulation rate across tionate effect of biodiversity change on local function in space. The value of b1 determines the shape of the species communities with high bi (electronic supplementary material, area relationship (SAR), with high and low values result- figure S2b) causes the average bA to increase. This variation ing in steep and shallow SARs, respectively (electronic at local scales, and the averaging out at larger spatial scales supplementary material, figure S3). causes the R 2 of the BEF relationship to increase in a saturating We find that the BEF slope, bA, increases with spatial manner with spatial scale (figure 2b). This occurs because vari- scale, peaking at intermediate scales and then decreasing ation in bi adds noise to the local BEF relationship, but this (figure 3a), while variation in ecosystem functioning explained noise is reduced by aggregating patches at larger spatial scales. by g-diversity decreases as a saturating function with spatial In this case, variation in the slope of the BEF relationship scale (figure 3b). These effects are consistent across all values across patches, bi, leads to a positive saturating increase in bA of bi , although the pattern flips when bi is negative (electronic with spatial scale (figure 2), now, because of nonlinear aver- supplementary material, figure S4). The magnitude of the aging of the contribution of each patch to regional ecosystem peak in bA and the scale over which it occurs depends on the functioning. A given change in a richness in patches with degree of compositional turnover. The peak is greatest, and higher than average bi will have a proportionally greater occurs at the lowest spatial scales, when compositional turnover effect on regional ecosystem functioning than would the is low. It decreases in magnitude and shifts successively same change in a richness in patches with lower than average to higher spatial scales as compositional turnover increases. bi (electronic supplementary material, figure S2b). Therefore, Likewise, the scale at which bA falls below b1 increases as com- as we aggregate across larger and larger regions, we are more positional turnover increases. The magnitude and speed of the likely to capture a greater range of bi, and so changes in g decline in R 2 of the BEF relationship with spatial scale also richness result in a greater proportional change in ecosystem increase as functional turnover increases. When compositional functioning than they would at small scales. turnover is high, bA remains equal or greater than b1 across the full range of scales considered. Consistent with case I above, with complete compositional turnover (b1 ¼ 0), there is no effect of scale on bA. 7. Case V: b-diversity across patches In this case, the slope of the BEF relationship changes with In case V, we explore how different levels of b-diversity affect spatial scale because changes in mean a richness do not result how the BEF relationship changes with spatial scale. We in the same proportional change in g richness. Incomplete com- again assume equal Si across all patches in a region, and an positional turnover drives two mechanisms, which together equal bi. However, across replicate simulated regions, Si is determine how b changes with spatial scale: (1) proportional drawn from the same normal distribution as in the previous changes in g richness are always less than proportional changes cases. We now determine the regional species richness SA by in a richness; and (2) the correlation between a and g richness combining the patches in the region, one at a time, and deter- becomes weaker at larger spatial scales. mining whether species are already present in the region Mechanism 1 causes bA to increase with spatial scale by p ¼ eðb0 þb1 SA^ Þ =ð1 þ eðb0 þb1 SA^ Þ Þ, where p is the binomial prob- because it means that a proportional change in g richness ability that a new species is unique. b0 and b1 are the intercept allows for a greater change in ecosystem functioning compared and slope of the logit model that determines the rate at which with the same proportional change in a richness. This mechan- this probability decreases with SA^ , the number of species in the ism is strongest when b-diversity is low. Mechanism 2 causes  bA to decrease with spatial scale because it erodes the signal of relationship, bA, as we did in our simulations. We then esti- 6 the local BEF relationship, which is ultimately responsible for mated the median and interquartile range of bA across all rspb.royalsocietypublishing.org driving bA at larger spatial scales. This is because as spatial years. To isolate the mechanism causing the BEF relationship scale increases, so does the range of possible a-richness to change with spatial scale, we performed the same simu- values that can lead to the same g-richness. For example, in a lation, but ignored the species identities and assumed that all simulated region of 10 patches, a g-richness of 10 can result patches contained unique species. This eliminates incomplete from each patch containing a single unique species, from all compositional turnover (i.e. maximizes compositional turnover) patches containing the same 10 species, and every scenario in and makes g-diversity unbounded. between these extremes. However, these two regions would Consistent with our theoretical predictions from cases III– have very different levels of biodiversity-driven ecosystem V, we find that bA increases with spatial scale but appears to functioning, and so the relationship with g-richness is weak. saturate at the largest spatial scale considered (figure 4a, dark The strength of mechanism 2 increases with spatial scale as grey). At the same time, the BEF relationship R 2 decreases Proc. R. Soc. B 285: 20180038 the number of possible values of a-richness for each level of asymptotically as scale increases, reaching close to zero by g-richness increases. However, it depends on the level of 2000 m2 (roughly 24 aggregated patches; figure 4b). Because b-diversity in the region; greater b-diversity leads to a slower the limited species pool in the experiment precludes us from increase in the strength of the mechanism with spatial scale. estimating the BEF relationship at larger spatial scales, we Because of mechanism 2, the local signal of the local BEF cannot determine whether the saturation in the slope, bA, and the R 2 of the BEF relationship is due to biological processes Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 relationship can be lost entirely at larger spatial scales and so bA decreases to 0. Together these two mechanisms cause bA or constraints to do with the size of the species pool. From to have a hump-shaped relationship, increasing at small spatial this test of the model predictions, we cannot identify which scales because the first mechanism dominates, but then of the three mechanisms (cases III–V) may or may not be decreasing at larger scales as the second mechanism increases responsible. We hypothesize that the pattern is dominated in strength. by the effect of compositional turnover as the increase in bA A variation on the pattern described above occurs when disappears when we assume that local patches have unique compositional turnover is extremely low (p ¼ 0.75). In this species (figure 4a, light grey). By assuming that patches have case, the decline in bA at large spatial scales is slower than it unique species, we remove any compositional overlap and so is when compositional turnover is higher (figure 3a) because artificially elevate g-diversity. only low-diversity plots have differences in composition and so g richness increases in regions with low, but not high a-diversity. As in the other cases, this results in a steepening of the BEF relationship. However, this increase in slope persists 9. Discussion at larger spatial scales as most of the changes in g richness Our work shows that the relationship between BEF is expected occur when the first few low-richness plots are aggregated. to change with the spatial scale at which it is observed. We have identified three mechanisms which can drive this scale depen- dence: (i) variation in species richness across local habitat 8. Experimental comparison using the Cedar patches in a region (case III), (ii) variation in the strength of the BEF relationship across local habitat patches in a region Creek biodiversity experiment (case IV) and (iii) incomplete compositional turnover across We analysed the Cedar Creek biodiversity experiment data local habitat patches in a region (case V). The first two mech- for above-ground biomass to ask whether the BEF relation- anisms result from the effects of noninear averaging on ship across spatial extents observed in this empirical dataset the effects of diversity change in low- versus high-diversity is consistent with our theoretical predictions. The experiment patches or in patches with weak versus strong BEF relation- included 18 species, drawn at random in combinations of 1, 2, ships, respectively. The third mechanism, which results in the 4, 8 and 16 species. To estimate the BEF relationship at spatial greatest changes in the BEF relationship with spatial scale, scales larger than a single plot (81 m2) we simulated regions results from the fact that, with incomplete compositional turn- by randomly drawing individual plots. We excluded all plots over, proportional changes in mean g-diversity are always less with 16 species because g-diversity in the experiment was 18 than proportional changes in mean a-diversity, and because so these plots contained almost all species. Excluding these the explanatory power of g-diversity on changes in ecosystem plots allowed for a greater range of species richness resulting functioning decreases with spatial scale. When the local BEF from random draws of plots at each spatial scale. We also relationship is positive and decelerating, as has generally excluded all plots that were not sampled in all years of the been found in experiments and often in field observations experiment. We then summed the above-ground biomass [1,2,28], the first two mechanisms cause the slope of the BEF and calculated the number of unique species that were planted relationship to increase with spatial scale, while the effect of in the plots in each simulated region. We did this at increasing the third mechanism is more variable and depends on spatial extents spanning the range from one plot to 30 plots. We the degree of compositional turnover present in a region. How- were unable to extend our approach to larger spatial scales ever, in general, the slope of the BEF relationship is expected to because the limited species pool in the experiment caused increase with spatial scale when species accumulate with space, g-diversity to converge on one or two values at all larger but to decrease when this accumulation has saturated. Overall, spatial scales, preventing us from accurately estimating bA. our findings provide an expectation for a scale-dependent We repeated this process 5000 times at each spatial scale for BEF relationship. each of the 14 years in the dataset. For each random draw of While all three mechanisms probably contribute to scale plots, we estimated biomass and species richness in each year dependence of the BEF relationship, we expect that effects of of the experiment. We then estimated the slope of the BEF incomplete compositional turnover should dominate the  (a) (b) 7 rspb.royalsocietypublishing.org 0.8 actual unique species R2 of BEF relationship 0.7 0.2 BEF slope (bA) 0.6 0.5 0.1 0.4 Proc. R. Soc. B 285: 20180038 0.3 0 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 scale (m2) scale (m2) Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 Figure 4. The strength of biodiversity effects, bA ( panel (a)) and R 2 ( panel (b)), estimated for each year at each spatial scale in the Cedar Creek simulations. Values of bA were estimated at each spatial scale by drawing 5000 replicate combinations of local plots and estimating the slope of the relationship between their aggre- gated biomass and their unique species in log – log space. Local scale plots were 81 m2 and the max spatial scale corresponds to 24 plots aggregated together. The dark grey shows the results when the species identities from the experimental plots were used to estimate g-diversity, the light grey shows the results when all plots were assumed to consist of unique species. The solid line indicates the mean across all 14 years in the dataset, the bands show the interquartile range. other two scaling mechanisms in most landscapes, and in con- peak at small spatial areas, and then decline with increasing trolled experiments, such as in our analysis of the Cedar Creek area sampled. In such scenarios, the slope of the BEF relation- experiment. This expectation is based on the fact that the over- ship at large spatial extents may be lower than it is at local all magnitude of change in the BEF slope was up to 10 times scales. By contrast, when the SAR is slow to saturate with greater in case V (incomplete compositional turnover) than in increasing area sampled, we expect the BEF slope to rise more the other two cases. These findings are consistent with other slowly with increasing area, and remain high over a larger studies that have demonstrated the importance of b-diversity range of scales. This pattern is consistent with what has been in maintaining ecosystem functioning at landscape scales shown for the invariability–area relationship , where [16–19]. Our approach with the Cedar Creek experimental invariability in ecosystem functioning increases at local scales, data, where we created simulated landscapes by randomly but then saturates at regional scales. Of course, our analysis combining local scale plots, is similar to the approach taken has not considered very large scales (e.g. continental), where by Pasari et al. , van der Plas et al. and Hautier et al. species area relationships have been found to steepen again. However, we have done this over a range of spatial [33,34]. In this case, we may expect to see a further steepening scales, rather than at a single larger spatial scale, which of the BEF relationship, as has been found for invariability. allows us to more generally consider how scale affects the Our results provide a theoretical expectation for how and BEF relationship. Further to the direct effect of incomplete com- why the BEF relationship should depend on the spatial scale positional turnover, these empirical studies also highlight how at which it is observed. They suggest that the way that this spatial differences in dominance and productivity of individ- scale dependence is realized in real ecosystems will depend ual species can lead to positive BEF effects. Such spatial on the abiotic and biotic processes that determine patterns insurance effects may be important drivers of how the BEF of a-, b- and g-diversities, and rates of ecosystem functioning. relationship changes with spatial scale , and are not Because our goal was simply to develop basic theoretical included in the mechanisms that we have identified here. In expectations for how and why the BEF relationship should addition, variation in local species richness or the local slope change with spatial scale, we intentionally omitted these of the BEF relationship would increase the relative contri- processes from our simulations. Further development of butions of these scaling mechanisms (cases III and IV) to BEF this theory should now consider the ecological processes scaling across the landscape. However, this variation would that would determine how this scale dependence will play have to be much higher, probably unrealistically high, to out in real landscapes. An obvious next step would be to reveal effects on scale dependence of similar magnitude as develop expectations for how environmental heterogeneity those that can result from incomplete compositional turnover. across time and space , should cause the BEF relationship Therefore, we expect that increases in the BEF slope are most to change with scale, via its influence on b-diversity. Linking likely to be driven in nature by the loss or gain of species that these expectations to observed environmental conditions and are shared across multiple local sites in a region, because it is patterns of b-diversity could give us predictions about how these shared species that cause changes in a-diversity to the BEF relationship should change with scale in a given proportionally exceed changes in g-diversity. landscape. Furthermore, metacommunity processes such as The range of spatial areas over which we should expect the dispersal and spatial connectivity are known to alter patterns slope of the BEF relationship to increase remains unresolved of biodiversity, community composition and the strength of and depends in part on the shape of the SAR, which captures the BEF relationship (i.e. spatial insurance) [13,20,21]. how species accumulate across space. When the SAR saturates For example, dispersal can both increase and decrease quickly, reflecting low b-diversity, we expect the BEF slope to b-diversity, and can be an important determinant of the  number of species present at a site, and whether they are well where the BEF relationship is not scale dependent. However, 8 suited to the prevailing conditions. Future progress on this few (if any) real landscapes conform to these strict criteria. rspb.royalsocietypublishing.org topic should focus on understanding how these ecological For example, environmental heterogeneity, spatial distance processes lead to scaling via the mathematical mechanisms and stochastic factors lead to compositional turnover in space that we have identified in this paper.. Therefore, we suggest that our expectation should be for To accomplish these goals, we see three obvious avenues a scale-dependent BEF relationship. forward: (i) simulation models that allow composition to Our findings provide a theoretical understanding of depend on the local abiotic and biotic conditions in each why we should expect the BEF relationship to vary with the local habitat patch and the spatial gradients of abiotic spatial scale at which it is observed. They suggest that we conditions and connectivity of patches; (ii) experiments that cannot simply apply our current understanding, which is vary a, b and g-diversities across a range of spatial scales; almost exclusively based on small-scale experiments and obser- and (iii) field observations (e.g. [17,28]) where community vations [1–3], to understand the consequences of biodiversity Proc. R. Soc. B 285: 20180038 properties can be nested at multiple spatial scales to estimate change at larger spatial scales without a theoretical framework the effects of biodiversity on ecosystem functioning. Simulation for scale-dependent change. We have identified three mechan- models offer the opportunity to further develop this theory and isms, which drive the scale dependence of the BEF relationship. incorporate additional ecological processes that we expect We must now apply this understanding to real landscapes , should be critical in determining the scale dependence of the where patterns of biodiversity, composition and ecosystem Downloaded from https://royalsocietypublishing.org/ on 07 August 2024 BEF relationship. Experiments and field observations offer functioning are determined by abiotic and biotic gradients, the opportunity to test and isolate the contribution of the rather than by statistical probabilities, as they are in our different scaling mechanisms that we have identified here. Fur- simulations. To this end, our findings provide an important thermore, field experiments can also provide information on step towards linking our understanding of BEF science to the how much natural ecosystems vary in the parameters that we spatial scales that are relevant to conservation decisions and have found to be drivers of BEF scaling. In our analysis of the provisioning of ecosystem services. the Cedar Creek data, we were able to identify that the BEF relationship depends on spatial scale in empirical data, but unfortunately this experimental design did not allow us to Data accessibility. All code required to reproduce the findings in tease apart the specific contributions of the different mechan- this paper are available at https://github.com/plthompson/BEF- spatial-scaling. The Cedar Creek data are available at https:// isms, although it appears that compositional turnover is www.cedarcreek.umn.edu/research/data. important. However, experiments could be designed that Authors’ contributions. All authors jointly conceived of the research ques- would differentiate these mechanisms, for example, by manip- tion. P.L.T. designed and performed the simulations and data ulating levels of spatial compositional turnover in species analyses and wrote the initial draft of the manuscript. All authors richness. Variation in the slope of the BEF relationship could contributed to revisions. be removed by calculating the predicted levels of ecosystem Competing interests. We declare that we have no competing interests. functioning in each local patch based only on the number of Funding. P.L.T. is supported by Killam and NSERC Postdoctoral Fellowships. A.G. is supported by the Liber Ero Chair in Biodiversity species in the patch, and the average BEF relationship in Conservation, NSERC Discovery grant. This research was partially the experiment. supported by funding from the Quebec Centre for Biodiversity Our conclusion, that a scale-dependent BEF relationship is Science. M.L. was supported by the TULIP Laboratory of Excellence expected, extends previous theory on this topic , which (ANR-10-LABX-41) and by the BIOSTASES Advanced Grant, funded suggested that the BEF relationship should be consistent by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement across spatial scales. In their case, Cardinale et al. assumed no. 666971). M.I.O. is supported by an NSERC DG. F.I. was sup- that all species in a region were present in every patch, and that ported by the US National Science Foundation’s Long-Term the strength of the BEF relationship was constant across Ecological Research (LTER) program (DEB-1234162) and the LTER patches. These assumptions meet those of case I in our study, Network Communications Office (DEB-1545288). References