A graphical null model for scaling biodiversity–ecosystem functioning relationships

Global biodiversity is declining at rates faster than at any other point in human history. Experimental manipulations at small spatial scales have demonstrated that communities with fewer species consistently produce less biomass than higher diversity communities. Understanding the consequences of the global extinction crisis for ecosystem functioning requires understanding how local experimental results are likely to change with increasing spatial and temporal scales and from experiments to naturally assembled systems. Scaling across time and space in a changing world requires baseline predictions. Here, we provide a graphical null model for area scaling of biodiversity–ecosystem functioning relationships using observed macroecological patterns: the species–area curve and the biomass–area curve. We use species–area and biomass–area curves to predict how species richness–biomass relationships are likely to change with increasing sampling extent. We then validate these predictions with data from two naturally assembled ecosystems: a Minnesota savanna and a Panamanian tropical dry forest. Our graphical null model predicts that biodiversity–ecosystem functioning relationships are scale‐dependent. However, we note two important caveats. First, our results indicate an apparent contradiction between predictions based on measurements in biodiversity–ecosystem functioning experiments and from scaling theory. When ecosystem functioning is measured as per unit area (e.g. biomass per m2), as is common in biodiversity–ecosystem functioning experiments, the slope of the biodiversity ecosystem functioning relationship should decrease with increasing scale. Alternatively, when ecosystem functioning is not measured per unit area (e.g. summed total biomass), as is common in scaling studies, the slope of the biodiversity–ecosystem functioning relationship should increase with increasing spatial scale. Second, the underlying macroecological patterns of biodiversity experiments are predictably different from some naturally assembled systems. These differences between the underlying patterns of experiments and naturally assembled systems may enable us to better understand when patterns from biodiversity–ecosystem functioning experiments will be valid in naturally assembled systems. Synthesis. This paper provides a simple graphical null model that can be extended to any relationship between biodiversity and any ecosystem functioning across space or time. Furthermore, these predictions provide crucial insights into how and when we may be able to extend results from small‐scale biodiversity experiments to naturally assembled regional and global ecosystems where biodiversity is changing.


| INTRODUC TI ON
Worldwide, drastic environmental changes are leading to biodiversity loss at regional and global scales (Millenium Ecosystem Assessment, 2005;Newbold et al., 2015;Tittensor et al., 2014).
Many predict that the rate of species loss will accelerate in the coming decades (Pereira et al., 2010;Pimm et al., 2014). Problematically, biodiversity supports many crucial ecosystem functions and services such as biomass production, carbon sequestration, and nitrogen cycling and retention Hooper et al., 2005;Weisser et al., 2017). That is, in small-scale experiments, ecosystem functioning increases with increasing species richness (Isbell et al., 2015;Reich et al., 2001Reich et al., , 2012Roscher et al., 2004;Schnitzer et al., 2011;Tilman et al., 2001;van Ruijven & Berendse, 2003;Zhang et al., 2012). In naturally assembled ecosystems, species richness appears to have an even larger effect on some functions like biomass production Liang et al., 2016;van der Plas, 2019;Zhang et al., 2012). Thus, ongoing species loss will likely have serious consequences for how ecosystems function (Cardinale et al., 2011;Hooper et al., 2012).
Yet, the results from local-scale experiments and observational studies may not be representative of the regional and global scales where species loss is occurring for several reasons (Craven et al., 2020;McGill et al., 2015). First, species loss in biodiversity experiments is often simulated randomly, non-random loss may have larger or smaller impacts on ecosystem functioning than random species loss (Isbell et al., 2008;Komatsu et al., 2019). Second, the relevant processes determining both biodiversity maintenance and ecosystem functioning change with increasing scale. Local processes like niche differentiation and plant soil-feedback may be most relevant at small scales while processes like dispersal within a meta-community are likely more important at regional scales (Leibold et al., 2017). Finally, the patterns of biodiversity change are likely different across scales with changes in community composition and species identity changing the most at small scales and loss dominating at the global scale (Blowes et al., 2019;Dornelas et al., 2014;McGill et al., 2015). Understanding how to link localscale results to the regional and global scales where species loss is occurring is crucial for predicting the consequences of the global extinction crisis.
Furthermore, theory predicts that the relationship between biodiversity and ecosystem functioning is likely scale-dependent.
A recent review on the theory predicting scale dependency of biodiversity-ecosystem functioning relationships suggests that there are four non-mutually exclusive process-based reasons for scale dependence (Gonzalez et al., 2020). First, based on coexistence theory, as spatial scale increases the opportunity for niche partitioning increases. This, in turn, increases the importance of species richness for ecosystem functioning. Second, if species abundance responds to environmental fluctuations in space and time then environmental autocorrelation determines the saturation rate of biodiversity-ecosystem functioning relationships. If this environmental autocorrelation is strong, then biodiversity is more important for ecosystem functioning at larger spatial scales.
Third, connectivity between local communities embedded within a meta-community can increase or decrease the reliance of ecosystem functioning on biodiversity across scales. Finally, the structure of networks within and among local communities is likely to result slope of the biodiversity ecosystem functioning relationship should decrease with increasing scale. Alternatively, when ecosystem functioning is not measured per unit area (e.g. summed total biomass), as is common in scaling studies, the slope of the biodiversity-ecosystem functioning relationship should increase with increasing spatial scale. Second, the underlying macroecological patterns of biodiversity experiments are predictably different from some naturally assembled systems.
These differences between the underlying patterns of experiments and naturally assembled systems may enable us to better understand when patterns from biodiversity-ecosystem functioning experiments will be valid in naturally assembled systems.

4.
Synthesis. This paper provides a simple graphical null model that can be extended to any relationship between biodiversity and any ecosystem functioning across space or time. Furthermore, these predictions provide crucial insights into how and when we may be able to extend results from small-scale biodiversity experiments to naturally assembled regional and global ecosystems where biodiversity is changing.

K E Y W O R D S
grasslands, productivity, species richness-area relationship, statistical scaling, upscaling in scale dependence of relationships between diversity and ecosystem functioning.
Importantly, understanding the potential impact of biodiversity loss on ecosystem functioning at the global scale requires two different types of scaling. First, biodiversity experiments often occur at small spatial scales. To apply this research to the global problem of biodiversity loss, we must understand how these relationships are likely to change with increasing sampling extent (i.e. from 1 m 2 to 1 km 2 , area scaling). Second, biodiversity experiments are abstractions of naturally assembled systems where diversity is manipulated, and other factors are held constant. Understanding the consequences of biodiversity change on ecosystem functioning also requires a conceptual scaling from experiments to naturally assembled systems. Recent work has attempted to bridge these gaps in our understanding with conceptual advances (Gonzalez et al., 2020;Isbell et al., 2017), empirical assessments using observational biodiversity gradients at multiple spatial scales (Chisholm et al., 2013;Craven et al., 2020;Luo et al., 2019) and theory Thompson et al., 2018). The major take-home message of most of these studies has been that the importance of biodiversity should increase with increasing area sampled (reviewed by Gonzalez et al., 2020).
Here, we provide a simple graphical null model that uses speciesarea curves and biomass-area curves ( Figure 1) to predict how biodiversity-ecosystem functioning relationships should change with increasing sampling extent. The goal of this graphical null model is to provide a simple scaling framework that enables us to use small-scale experimental and observational datasets and extrapolate these patterns across spatial scales that are unmeasured.
This extrapolation represents the expected real-world scaling relationship, given that the patterns we observe at these small scales hold across others. This graphical null model focuses primarily on scaling up with increasing sampling extent but may have important implications for how to scale from experimental to naturally assembled systems.
F I G U R E 1 Graphical depiction of the underlying math of the plant species richness-biomass relationship across different spatial scales. (a) Locally, species richness has a nonlinear relationship with sampling extent while biomass has a linear relationship with sampling extent. For a specific area A i , the species-area curve represents the average number of species S i found on plots of area A i . Similarly, for a specific area A i the biomass area curve yields the average total biomass B i of species found on that area. The corresponding coordinate pairs (S i ,B i ) form the graph of the species richness biomass production relationship (b). These coordinate pairs reveal an underlying accelerating nonlinear relationship between species richness and biomass production across scales. At any given scale, the species richness biomass production curve contains the origin (because at zero species richness, biomass production is zero and vice versa), as well as the point (S i ,B i ) of measured averages. If the species richness biomass relationship is linear (see Figure S1 for alternative interpretations of this shape) and goes through these two points, then the predicted slope of the species richness-biomass relationship is the average biomass for that area divided by the average species richness for that area (c) and the slope of the species richness biomass production relationship increases with increasing spatial scale of sampling. (d) When biomass is measured in terms of per unit area, it is not likely to change with increasing sampling extent. (e) If this is the case, then the species-area curve determines the slope of the species richness-biomass relationship with increasing spatial scale (f). Thus, we expect the per unit area version of the species richness-biomass relationship to have a decreasing slope with increasing area and that the slope will decrease at a decelerating rate. See Supporting Information S1 for mathematical derivation of these predictions 2 | MATERIAL S AND ME THODS

| Deriving the graphical null model
If the species-area relationship is nonlinear (Braun-Blanquet, 1932;Cain, 1938;Connor & McCoy, 1979;Gleason, 1922;Lomolino, 2000Lomolino, , 2001Preston, 1960;Rice & Kelting, 1955;Rosenzweig, 1995;Figure 1a), then when species richness and ecosystem functioning are positively related the relationship between the two must change with increasing sampling extent. This change occurs because species richness and ecosystem functioning independently depend on area but not at the same rate (Figure 1a,b; Supporting Information S1).
Whether the slope of a biodiversity-ecosystem functioning relationship increases or decreases depends on how that ecosystem function scales. If an ecosystem function increases linearly with increasing sampling extent (Figure 1a; Figure S1), then the apparent contribution of each additional species to ecosystem functioning will increase with increasing sampling extent. Regardless of the processes involved, the slope of this biodiversity-ecosystem functioning relationship will increase, as predicted by current scaling theory ( Figure 1c, as predicted by Gonzalez et al., 2020;Isbell et al., 2017;Thompson et al., 2018). However, ecosystem functioning in biodiversity-ecosystem functioning studies is commonly standardized per unit area (e.g. biomass per m 2 , e.g., Isbell et al., 2011;Meyer et al., 2016;Reich et al., 2012;Roscher et al., 2004;Spehn et al., 2005;van Ruijven & Berendse, 2003;Weisser et al., 2017).
These data are standardized to ease comparison between plots and across scales when relevant (Roscher et al., 2004(Roscher et al., , 2005. In per unit area terms, functioning should be invariant to area (Figure 1d). If ecosystem functioning is invariant to area, then we predict that the slope of the biodiversity-ecosystem functioning relationship will decrease with increasing sampling extent (Figure 1e,f).

| Empirically validating the graphical null model
To validate this model, we used empirical data from two sites: a naturally assembled savanna at Cedar Creek Ecosystem Science Reserve

| Site description: Cedar Creek Ecosystem Science Reserve
We collected data to empirically validate this null model in a naturally assembled savanna at CCESR in July of 2012. Our sites were located on Typic Udipsamment soils (Dickie et al., 2007) with a mean annual precipitation of 77.60 ± 4.57 cm (95% confidence interval for 1962-2012). CCESR consists of a matrix of prairies and savannas. We used two oak savannas (predominantly Quercus macrocarpa) of approximately equal type, age and fire frequency. We chose oak savannas with similar burning and establishment histories that had little woody biomass and coarse woody debris to maximize similarity in nutrient availability and community composition between sites.
The two oak savannas we chose were both burned in 2011 after the summer growing season; therefore, all herbaceous biomass was considered to be the current year's growth. Both savannas were burned frequently, 9 out of every 10 years and 4 out of every 5 years, respectively, and were established 48 years ago from abandoned farmland (Peterson & Reich, 2001). The dominant tree species at both sites were bur oak Quercus macrocarpa and northern pin oak

Quercus ellipsoidalis.
We established five 20 m by 20 m plots across these two oak savanna sites (three in one and two in the other). We avoided portions of the savannas with high amounts of woody biomass because many of the common woody savanna species are resistant to fire and we could not be certain that any woody biomass was indicative of this season's growth. We established all plots >3 m from any adult trees.
Each plot was divided into 25 4 m × 4 m subplots.

| Biomass harvest: CCESR
To measure biomass, we harvested the above-ground vegetation in Vegetation was refrigerated when not being sorted. Once sorted, we placed each species in a paper envelope and dried all samples for at least 10 days in a 40°C drying oven. After 10-15 days, we weighed all samples and recorded the dry biomass.

| Site description: Barro Colorado Island 50-Hectare Plot
In addition to data collected at CCESR, we used data collected in the Barro Colorado Island 50-hectare plot from the 2015 census to validate this method outside of a grassland ecosystem (see Condit, 1998 for data collection methods).

| Barro Colorado Island 50-hectare plot: Subsampling method
We assigned each 5 × 5 m subplot within the BCI 50-hectare plot a random number. We then randomly subsampled (using the command 'sample') these subplots to assemble non-spatially contiguous subsamples of the following areas: 25, 50, 75, 100, 150, 200,225, 300 and 400 m 2 . We randomly selected 70 of each subsample size for inclusion in our analysis. We used 70 as our sample size because it allowed for the inclusion of 400 m 2 sampling areas which are visually comparable to the largest sampling area included in our data from CCESR. We subsampled 70 of smaller sampling extents to prevent differences in sample size which may be relevant for our dataset from CCESR where sample sizes varied from 125 at the smallest sampling extent to 5 at the largest sampling extent due to the nested nature of our sampling design. Furthermore, we used spatially noncontiguous subsamples because the shape of the species-area curve changes between randomly sampled landscapes and landscapes where nested samples have been taken (Rosenzweig, 1995). Thus, to ensure that our results were not dependent upon the nested sampling schema used at CCESR, we subsampled randomly within the BCI 50-ha plot.

| Data analysis
To validate whether the observed relationship between species richness and biomass production was dependent on sampling extent in the way we predict, we also examined the relationship of both species-richness and biomass production individually with increasing sampling extent. We fit linear and 'Michaelis-Menten' (nonlinear) models to our data, two commonly used functions for describing biodiversity-ecosystem functioning relationships (Hooper et al., 2005;Roscher et al., 2004;Tilman et al., 2001), and selected the model with the lower Akaike Information Criterion (AIC) to determine whether the relationships between biomass production and area, species richness and area, and species richness and biomass production were linear or nonlinear. When we were unable to calculate parameters for a Michaelis-Menten curve using the 'SSmicmen' command in R and an 'nls' command would not converge with uniform starting values (i.e. 1 for both parameters), we assigned an NA to the nonlinear model and considered the linear model to be a better fit.
We then used a mixed effects ANCOVA to determine whether the relationship between species richness and biomass production changed with increasing sampling extent using the command 'lme' when our sampling design necessitated the inclusion of random effects and 'gls' when random effects were not necessary. Our sampling design at CCESR was nested with subplots within plots across two sites. We therefore included a nested random effect of plot within site at Cedar Creek. Our randomized sampling design at BCI did not necessitate the inclusion of random effects. We tested all models for heteroscedasticity using a Breusch-Pagan test (Breusch & Pagan, 1979). When models had significant heteroscedasticity, we added a weighted variance structure to the model which was specified as an exponential function of the fitted value.
This variance structure reduced heteroscedasticity significantly resulting in a non-significant Breusch-Pagan test for all models.
All analyses were conducted in R Statistical Software and plotted using 'ggplot2'.

| Species richness-biomass relationships match macroecological expectations
We produced species-area relationships ( Figure 2) and biomassarea relationships (Figure 3)  Note that with regards to our results we use the term biomass to refer to our measures and biomass production to refer to the general ecosystem function. We distinguish between these two terms because while biomass is often a good proxy for biomass production, the latter requires a measure of turnover through time which we do not present here for BCI. We found that, at the relatively small spatial scales measured here, species richness was best fit with a nonlinear curve with increasing sampling extent at both sites (Table S1; Figure 2a,b). Alternatively, biomass increased linearly with increasing sampling extent at both sites (Table S1; Figure 3a,c). By contrast, biomass per m 2 was invariant to sampling extent at both sites (Table S1, Figure 2b,d).  Figure 4). That is, we assumed that when an ecosystem has no species it also cannot produce biomass. At each sampling extent, we examined the species richness-biomass relationship in two different ways: (a) the relationship between total species richness per plot and total biomass per plot and (b) the relationship between total species richness per plot and biomass per m 2 (e.g. Isbell et al., 2011).
In accordance with our predictions, the slope of the species richness-total biomass relationship increased with increasing sampling extent at both sites (Figure 5a,c). In fact, no predicted slope for the species richness-biomass relationship was outside of the confidence limits of the observed slope (Table S2). Mathematically, F I G U R E 3 Biomass-area curves in a Minnesota savanna and the Barro Colorado Island 50-hectare plot. (a) When biomass was summed, it scaled linearly with increasing area. (b) When biomass was examined per m 2 (as is commonly done in biodiversity-ecosystem function experiments), there was no significant relationship between biomass and plot area. (c) When biomass was summed in randomized subsamples of the BCI 50-ha plot, it increased in a linear fashion with increasing plot area. (d) Biomass per m 2 in the BCI 50-ha plot did not significantly change with increasing plot area

F I G U R E 4 Predicted slopes of the species richness-biomass relationships at each scale for both absolute biomass and biomass per m 2 .
When biomass is summed, we predict that the species richness-biomass production relationship will have an increasing slope with increasing area, for both recently burned savannas (a) and the randomized subplots of the BCI 50-ha plot (c). For biomass per m 2 , we predict that the slope of the species richness-biomass production relationship will decrease with increasing area for both recently burned savannas (b) and the BCI 50-ha plot randomized subplots (d) additional species contributed more to biomass production at the largest sampling extent, whereas the smallest sampling extent had the smallest relative gains in total biomass with each additional species. In terms of biomass per m 2 , we found that the slope of the species richness-biomass relationship decreased with increasing sampling extent (Figure 5b,d). That is, mathematically, each additional species contributed less to biomass per m 2 at the largest sampling extent. The smallest sampling extent had the largest relative gains in biomass per m 2 with each additional species. At BCI, the slope decreased significantly with increasing sampling extent, while at Cedar Creek this change was not statistically significant (Table S3).

| D ISCUSS I ON
Our primary goal in this study is to test whether a simple graphical null model-which leverages the widely observed and broadly studied relationships between biodiversity and area, and between biomass and area-can be used to characterize how relationships between biodiversity and ecosystem functioning change across spatial scales. We predict that the null expectation for the observed relationship between biodiversity and ecosystem functioning should be both scale-dependent and predictable, based on the difference in shapes between the largely nonlinear species-area relationship, versus the largely linear biomass-area relationship.
Importantly, the resulting graphical null model relating biodiversity and functioning does not imply a causal link between the two, nor does it suggest that extrapolations made using such a null model will be accurate. For example, if the underlying environmental conditions are different between the scales, these predictions may not be ecologically relevant. Rather, this graphical null model demonstrates the necessary outcome of the shared underlying relationships between spatial scale and both biodiversity and ecosystem functioning. Thus, if a scaling relationship that is consistent with our graphical null model is observed in empirical data, this result should not be taken as an indication of specific underlying causal links or biologically relevant processes, but rather, as an expected outcome of well-established mathematical scaling relationships.

| Standardizing ecosystem functioning measures across sampling extents
Our null model predicts that the slope of the species richness-total biomass relationship should increase with increasing sampling extent and alternatively that the slope of the species richness-biomass F I G U R E 5 Actual species richness-biomass relationships and predicted species richness-biomass relationships in recently burned savannas and randomly selected subplots of the Barro Colorado Island 50-hectare plot. Throughout, dashed lines represent predictions and solid lines represent linear models from data. (a) When biomass is summed, the slope of the species richness-biomass production relationship increases with increasing spatial scale as predicted in five recently burned savannas. (b) When biomass was measured as biomass per m 2 , as is commonly calculated from biodiversity-ecosystem function experiments, the slope of the species richness-biomass production relationship decreases with increasing spatial scale. (c) Furthermore, from randomly subsampled subplots of the BCI 50-hectare plot, we found similarly that the species richness-biomass production (in terms of total above-ground biomass) relationship had an increasing slope with increasing spatial scale. (d) When biomass production was calculated as biomass per m 2 from randomly subsampled subplots of the BCI 50-hectare plot, we found that the slope of the species richness-biomass production relationship decreased with increasing spatial scale per m 2 relationship decreases with increasing sampling extent. The biomass per unit area (i.e. m 2 or ha) approach is the typical method of comparison in biodiversity-ecosystem functioning studies (e.g. Duffy et al., 2017;Isbell et al., 2011;Reich et al., 2001Reich et al., , 2012Roscher et al., 2004;Spehn et al., 2005). Scaling theory in biodiversity-ecosystem functioning research, however, relies almost entirely on the assumption that ecosystem functioning at larger spatial scales is the sum of ecosystem functioning at smaller scales (see Gonzalez

| Applying macroecological predictions from local to global scales
A critical conclusion from our findings is that we can use macroecological patterns to accurately predict the scaling of the species richness-biomass production relationship at local spatial scales.
However, the local to global scale species-area curve (both logtransformed) is likely triphasic following a power law. As the area considered expands, the rate of increase in species richness slows at local scales (as seen here). At regional scales, species richness increases approximately linearly with increasing area. At landscape to global scales, species richness increases at an accelerating rate with increasing area (Lomolino, 2001;Preston, 1960;Rosenzweig, 1995). Given this underlying macroecological pattern, we expect that the species richness-biomass relationship will also change as we move from local to global scales (when both species richness and biomass are log-transformed, Figure 6, see also (Thompson et al., 2018). In terms of biomass production ( Figure 6a), we expect that, as area increases, the slope of the relationship between species richness and biomass will increase locally, stay the same regionally and decrease globally. Importantly, this global slope will not decrease below the slope of the highest local slope (Figure 6b). That is, if species richness contributes to biomass production, each additional species will contribute proportionally the most to biomass production at regional scales and the least at the smallest local scale. Alternatively, in terms of biomass production per unit area, we expect that the species-area curve will change nonlinearly from local to global scales and the biomass production per unit area will be invariant to scale (Figure 6c). In this case, the slope of the relationship between species richness and biomass production will decline as area increases (Figure 6d).
The general scaling relationship found here can also be extended for any other analysis that combines a nonlinear relationship with F I G U R E 6 Graphical predictions for how the log(species richness)-log(biomass production) relationship will scale with increasing log(area) from local to regional and global scales depending on the way in which biomass production is calculated. (a) In log space, the species-area curve is thought to be triphasic as log(area) increases from local to regional and then global scales (black curve). The log-transformed total biomass production likely monotonically increases with increasing spatial scale (grey line). (b) The slope of the log(species richness)-log(biomass production) relationship increases with increasing spatial scale at local to the regional scales, does not change at regional scales and decreases slightly at global scales. This would indicate that the exponents of the non-log-transformed relationships increase with increasing spatial scale. (c) In log space, the species-area curve is triphasic (black line). Log-transformed biomass per unit area (e.g. per m 2 or per ha) is invariant to area. (d) The slope of the log(species richness)-log(biomass production) relationship decreases with increasing area at local, regional and global scales. This would indicate that the exponents of the non-log-transformed relationships decrease with increasing spatial scale a linear relationship. For example, the species richness-time relationship, while less well studied than the species-area relationship is also nonlinear (Adler & Lauenroth, 2003;Preston, 1960;White et al., 2006). If we assume that biomass production per unit area remains relatively constant over time in an equilibrium community, then we predict that the slope of the species richness-biomass relationship will decrease with increasing time and each additional species will be less important for biomass production per unit area with increasing time (Figure S2a,b).

| Experiments versus naturally assembled systems
We report here a mathematical outcome of combining the speciesarea curve and biomass-area curve. While mathematically simple, this outcome is useful for scaling biodiversity-ecosystem functioning relationships. Importantly, this outcome does not imply causality.
That is biodiversity need not cause enhanced ecosystem functioning in naturally assembled systems for this scaling to be valid. This mathematical scaling addresses one aspect of the type of scaling necessary to apply findings from biodiversity-ecosystem functioning experiments to naturally assembled systems-going from small area plots to larger areas. However, examining the underlying macroecological patterns may also highlight how experiments may operate under different assumptions from naturally assembled systems and therefore inform the second type of scaling required-going from experiment to naturally assembled system at the same spatial extent.
These differing underlying macroecological patterns may enable us to understand when results from biodiversity-ecosystem functioning experiments are likely to be most applicable to naturally assembled systems at the same spatial extent.
While the underlying macroecological patterns of biodiversity experiments may differ from naturally assembled systems in many ways, one important way is how biodiversity and biomass are likely to change through time . These differences in the underlying changes in biomass and biodiversity through time likely result in contrasting patterns in biodiversity-ecosystem functioning relationships through time. If we view biodiversity experiments through time, biomass per m 2 often increases on average (e.g. Reich et al., 2012;Weisser et al., 2017), Figure 2c). Furthermore, initial planted species richness is often reported in publications from biodiversity experiments rather than realized species richness (e.g. Reich et al., 2012;Weisser et al., 2017). This initial species richness is constant over time. If species richness remains constant (or decreases) through time and biomass per unit area increases, then the slope of the relationship between the two will increase with increasing time, as found by Reich et al. (2012, Figure S2d). Alternatively, in naturally assembled systems at the same spatial extent that are not undergoing succession, the general expectation is that biomass remains constant while species richness increases nonlinearly (the speciestime relationship, Adler & Lauenroth, 2003;Preston, 1960;White et al., 2006). If biomass remains constant while species richness increases nonlinearly, then we predict that over time each additional species will contribute less to biomass production. In biodiversity experiments, we expect the opposite ( Figure S2a,b).
Understanding these underlying changes may help us better apply biodiversity ecosystem functioning research to naturally assembled communities where biodiversity-ecosystem functioning relationships are more varied and often negative (van der Plas, 2019).
However, biodiversity change in naturally assembled systems differs from the biodiversity change mimicked by experiments in two crucial ways. First, biodiversity loss when it occurs in naturally assembled systems is non-random (Chen et al., 2020;Isbell et al., 2008). Rather, specific functional groups may be more or less prone to local extinctions than others (Harpole et al., 2016;Isbell et al., 2008;Komatsu et al., 2019;McDowell et al., 2020;Suding et al., 2005). The consequences of this non-random loss may be entirely different than the consequences of random loss. Second, species additions often offset local extinctions (Bannar-Martin et al., 2018;Blowes et al., 2019) and thus consistent net loss may be an unrealistic expectation in naturally assembled systems at local and regional spatial scales. More likely, communities will change in composition with some species increasing or decreasing in abundance and no concomitant changes to the identities of the species present (Leibold et al., 2017). These changes in abundance likely also have different effects on ecosystem functioning that are not predicted here. In fact, our goal here is to offer a framework for making null predictions about changes in biodiversity-ecosystem functioning expectations with scale.
When actual patterns deviate from these expectations, this offers an opportunity to explore the ecological mechanisms at play (e.g. biotic homogenization, species dispersing into habitats for which they are well suited, changes in species composition). For example, if biodiversity loss continues at the global scale while biotic homogenization continues to play out at regional spatial scales, then we expect that the species-area curve will begin to saturate rather than increase at global scales. If this happens then we predict that the importance of species richness for ecosystem functioning will decrease with increasing spatial scale.

| CON CLUS IONS
How biodiversity-ecosystem functioning relationships vary with increasing spatial and temporal scale has critical implications for how biodiversity is managed to provide ecosystem services. Current conservation, restoration and valuation efforts depend on our ability to scale from local experimental evidence to larger regional and even global scales (Isbell et al., 2017;Naidoo et al., 2008). The current capacity for upscaling is limited (Bockstael et al., 2000). Upscaling the local scale knowledge of diversity-productivity relationships to larger scales relevant to management is challenging because human-induced diversity changes represent a large variety of types of changes in addition to changes in species richness.
Our results reveal the mathematical inevitabilities of how biodiversity-ecosystem functioning relationships vary with sampling extent. For any ecosystem function that scales linearly with increasing area, the contribution of each additional species must increase with increasing sampling extent due to the nonlinear species-area relationship if there is a positive relationship between species richness and that ecosystem function. Furthermore, all measures of biodiversity and all measures of ecosystem functioning have individual relationships with space and time. The shape of these underlying macroecological relationships will determine how they scale.
Understanding these underlying scaling relationships may allow us to make generalizable predictions for the consequences of species loss on any well-supported biodiversity-ecosystem functioning relationship at global scales. Additionally, understanding how the units of the measures that we use alter relationships with scale enables us to make better predictions. Many management-relevant functions and services like net primary production and harvestable bole volume are used by managers in per unit area/volume or per unit time terms. This simple standardization which allows comparison between management units also reverses our expectations for how these functions change across spatial scale. Finally, applying this simple mathematical scaling framework may enable us to focus more precisely on the mechanisms that may cause these relationships to diverge from our theoretical expectations.

ACK N OWLED G EM ENTS
Funding for the collection of field data at Cedar Creek Ecosystem

AUTH O R S ' CO NTR I B UTI O N S
K.E.B. wrote all initial drafts, provided draft revisions, conducted all the analysis and coordinated the collaboration for this manuscript with in depth feedback from A.J.W.; The original concept for this study was devised by K.E.B., A.J.W., G.A.P. and I.G.L.; S.A.S. and P.B.R. provided feedback on the initial concept and designed and enabled fieldwork at CCESR; J.W.S. and K.Y. assisted in initial design, conducted fieldwork at CCESR and provided first analyses of data.
A.T.C., J.C., A.S.M. and L.W. provided feedback on initial concepts and significantly contributed to revision of the manuscript. All authors provided feedback on all drafts of the manuscript.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns. com/publo n/10.1111/1365-2745.13578.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data collected for the purpose of this analysis can be found in the Dryad Digital Repository https://doi.org/10.5061/dryad.3xsj3 txfk (Barry et al., 2020). Data from the BCI 50-Ha plot are publicly available http://ctfs.si.edu/webat las/datas ets/bci/. The Dryad Digital Repository also includes all codes used to produce the analyses and figures presented in this paper.