A cross-scale assessment of productivity-diversity relationships

Biodiversity and ecosystem productivity vary across the globe and considerable effort has been made to describe their relationships. Biodiversity-ecosystem functioning research has traditionally focused on how experimentally controlled species richness affects net primary productivity (S→NPP) at small spatial grains. In contrast, the influence of productivity on richness (NPP→S) has been explored at many grains in naturally assembled communities. Mismatches in spatial scale between approaches have fostered debate about the strength and direction of biodiversity-productivity relationships. Here we examine the direction and strength of productivity’s influence on diversity (NPP→S) and of diversity’s influence on productivity (S→NPP), and how this varies across spatial grains using data from North American forests at grains from local (672 m2) to coarse spatial units (median area = 35,677 km2). We assess relationships using structural equation and Random Forest models, while accounting for variation in climate, environmental heterogeneity, management, and forest age. We show that relationships between S and NPP strengthen with spatial grain. Within each grain, S→NPP and NPP→S have similar magnitudes, meaning that processes underlying S→NPP and NPP→S either operate simultaneously, or that one of them is real and the other is an artifact. At all spatial grains, S was one of the weakest predictors of forest productivity, which was largely driven by biomass, temperature, and forest management and age. We conclude that spatial grain mediates relationships between biodiversity and productivity in real-world ecosystems and that results supporting predictions from each approach (NPP→S and S→NPP) serve as an impetus for future studies testing underlying mechanisms. Significance statement The relationships between diversity and productivity are central to efforts explaining global variation of biodiversity and rates of carbon sequestration. However, little is known about the relative importance of biodiversity as the driving force, or as the consequence of ecosystem-level productivity. Our analysis of a comprehensive database of North American forests reveals that biodiversity and productivity can be both cause and effect and that their relationship strengthens with spatial grain. Importantly, we show that environmental context is more important in determining biodiversity and productivity than either biodiversity or productivity alone. Productivity-diversity relationships emerge at multiple spatial grains, which should widen the focus of national and global policy and research to larger spatial grains.

1 Overview of hypotheses predicting grain dependence of relationships between net primary 1 3 8 productivity (NPP) and species richness (S). No.
Direction Mechanism of grain dependence Weakens or strengthens towards coarse grain?

I NPP→S and S→NPP
Spatially asynchronous demographic stochasticity impacts small populations (or small grains) and averages out over large grains.
Both NPP→S and S→NPP strengthen towards coarse grains II NPP→S At larger grains, higher NPP is associated with increased heterogeneity and/or dissimilarity of local patches, allowing for greater regional coexistence.
NPP→S strengthens towards coarse grains (28,35,36) III NPP→S A statistical interaction between NPP and grain in their effect on S emerges as a consequence of increasing occupancy with NPP.
NPP→S weakens towards coarse grains (37) IV NPP→S At very large grains (thousands of km 2 and larger), high productivity increases occupancy and population size, thus increasing the probability of reproductive isolation and speciation NPP→S strengthens towards coarse grains V S→NPP Stochastic sampling effects dominate at small grains, resource partitioning at larger grains ('spatial insurance'), and their relative magnitude determines the grain dependency.
Both strengthening or weakening possible (39,40) VI S→NPP Functionally redundant species at the regional grain can compensate for low richness at local grains.
S→NPP strengthens towards coarse grains VII S→NPP With incomplete compositional turn-over, proportional changes in larger-grain richness are always less than proportional changes in smaller-grain richness such that the explanatory power of richness on changes in functioning decreases with spatial scale. S → NPP strengthens strengthens towards coarse grains until species richness saturates (42) Running title: Cross-scale diversity-productivity Craven et al. 7  Spatial patterns in productivity (NPP) and richness (S) emerged at coarser spatial grains, with higher S 1 4 5 and NPP usually observed in the eastern USA than in the western USA (Fig. 1). Biomass, a time-1 4 6 integrated measure of NPP that also influences diversity, also exhibited similar patterns (Fig. 1). Bivariate  show NPP as a response to S, panels C and D show NPP as a predictor. Solid lines are least-squares linear 1 5 3 regressions fitted at each grain, shaded areas are standard errors. Analyses were performed using stratified  Structural Equation Models (SEM). We examined relationships between species richness and net 1 5 8 primary productivity (NPP) across spatial grains using two SEMs for each spatial grain: the first 1 5 9 (S→NPP) testing the direct effect of S on NPP and the indirect effect of NPP on S (via biomass), and the 1 6 0 second (NPP→S) testing both the direct and indirect effects of NPP on S (Fig. 3). In both SEMs, 1 6 1 environmental variables (e.g., mean annual precipitation (MAP), mean annual temperature (MAT), Running title: Cross-scale diversity-productivity Craven et al. 9 temperature seasonality, and elevation range), size of the species pool, forest age, and management were 1 6 3 used to explain variation in S, biomass, and NPP. At the intermediate and coarse grains, we also included 1 6 4 area (of each spatial unit) to account for variation in species richness due to the effects of area (see 1 6 5 Methods). variables (e.g., mean annual precipitation, mean annual temperature, temperature seasonality, and 1 7 1 elevation range), size of the species pool, forest age, and management, in forests across the contiguous 1 7 2 USA at three spatial grains. Both models fit the data well at all spatial grains (P-value of the Chi-square 1 7 3 test > 0.1; Table S1). Boxes represent measured variables and arrows represent relationships among 1 7 4 variables. Solid blue and red arrows represent significant (P< 0.05) positive and negative standardized 1 7 5 path coefficients, respectively, and their width is scaled by the corresponding standardized path 1 7 6 coefficient. Solid and dashed gray arrows represent non-significant (P>0.05) positive and negative 1 7 7 standardized path coefficients, respectively. AGE is forest age, MANAGED is forest management, 1 7 8 ANN.PREC is mean annual precipitation, ANN.TEMP is mean annual temperature, TEMP.SEAS is 1 7 9 temperature seasonality, ELEV.RANGE is elevation range, S.POOL is the regional species pool, and 1 8 0 AREA is area. S, BIOMASS, NPP, and AREA were natural log transformed prior to analysis. grains. Points are standardized path coefficients and solid lines are 95% confidence intervals. Both models fit the data well for all spatial grains (P-value of the Chi-square test > 0.1; Table S1). At increasing spatial grain ( Fig. 3 & 4). At the fine spatial grain, we found a weak direct effect of S   Overall, the SEMs show that the productivity-diversity relationship increases in strength with spatial 2 0 2 grain, and both relationships (S→NPP and NPP→S) explain similar amounts of variation, albeit with 2 0 3 some differences in the direct and indirect effects. At fine spatial grains, our SEMs show greater support 2 0 4 for a strong indirect effect of NPP on S via biomass, but do not support the inverse effect of S on NPP. Towards coarser spatial grains, our SEMs do not conclusively show stronger support for one direction of  Table S1) and between biomass and NPP ( Fig. S3D, E, and F; Table S1). and NPP, and to provide an assumption-free alternative to the SEMs, we fitted two random forest models 2 1 0 for each of the three spatial grains: one with NPP and the other with S as response variables. We found 2 1 1 that NPP was an important predictor of S at the fine and intermediate spatial grains (Fig. 5A), with 2 1 2 unimodal and linear effects respectively (Fig. 5), but was less important relative to other predictors at the 2 1 3 coarse spatial grain. For S, we found that species pool, MAT, MAP, and forest age were the best coarse grains. For NPP, we found that species richness was one of the weakest predictors relative to other  and NPP at three spatial grains, which is the mean decrease in squared error caused by each of the  The first important result is the similar magnitude of the S→NPP (18) and NPP→S (10, 26, 43) 2 3 0 relationships at all grains. This reflects, in part, that both productivity and species richness have many 2 3 1 environmental and geographical drivers in common (44), which complicates distinguishing correlation 2 3 2 from causation, even when using SEMs (45, 46). There are two possible interpretations of this result: (i) it 2 3 3 may indicate that diversity's causal effects on productivity and productivity's causal effects on diversity 2 3 4 operate simultaneously, which was suggested by (18), but never demonstrated on observational data from experiments that manipulate diversity in ways that mimic biodiversity change (i.e. species gains and 2 3 9 losses) in real-world ecosystems (11, 48-50), we see little hope for resolving this with contemporary data 2 4 0 and approaches.

4 1
Our second important result is that both S→NPP and NPP→S strengthen from the fine to the intermediate 2 4 2 grain, and in the case of the SEM both relationships continue strengthening towards the coarsest grain.

4 3
While grain-dependent shifts are often expected (Table 1), this had not been shown previously with 2 4 4 empirical data for S→NPP using spatial grains coarser than several hectares (25, 31, 32). If the S→NPP 2 4 5 direction is the real causal one, then our results from SEM and RF analyses support several theoretical  (Table 1) and give further impetus to efforts quantifying biodiversity effects in naturally 2 4 7 assembled ecosystems at broad spatial scales (51). If the NPP→S direction is the real causal one, then our third possibility is that both NPP→S and S→NPP are real and that they operate simultaneously, as suggested by our SEM results. In this case, we are unaware of any theory that considers how this 2 5 2 reciprocal relationship would be expected to change with increasing spatial grain. The one caveat 2 5 3 applicable to interpreting any direction of diversity-productivity relationships is that of demographic 2 5 4 stochasticity (mechanism I in Table1), which may weaken both NPP→S and S→NPP, or their synergistic 2 5 5 interplay, at fine spatial grains. In our study, the strong local effect of demographic stochasticity appears 2 5 6 plausible given the small area of the forest plots (0.067 ha) and small population sizes (12.24 ± 0.02 trees 2 5 7 per plot; range = 1-157 trees per plot) therein. This would suggest that temporal changes in local scale  The third key result is that other predictors, such as temperature and biomass, were particularly influential 2 6 0 in all our analyses. That is, the grain dependence of the relationship between S and NPP was coupled with 2 6 1 a clear increase in the effect of annual temperature (but not precipitation) on both S and NPP towards 2 6 2 coarse grains, which supports the notion that either temperature-dependent diversification (55, 56) or Running title: Cross-scale diversity-productivity Craven et al. 14 ecological limits (43) shape diversity at these spatial grains. The consistently weak effect of precipitation 2 6 4 is expected since we focus on forests, which only grow above certain precipitation thresholds (57).

6 5
Second, we found a positive, indirect effect of NPP on species richness via forest biomass at the fine 2 6 6 spatial grain, which supports multiple hypotheses (Table 1) such as the view that higher ecosystem 2 6 7 productivity enhances species diversity by enabling larger numbers of individuals per species to persist 2 6 8 due to lower extinction rates (35,37,58), particularly at fine grains where stochastic extinctions occur functioning at any spatial grain. Our results reveal that mechanisms associated with one direction of diversity-productivity relationships biomass at the fine spatial grain provides support for the more individuals hypothesis (37), although it is 2 7 9 typically tested at regional to continental spatial scales. Increasingly, macroecological mechanisms such 2 8 0 as speciation gradients (60) and water-energy variables are being examined in small-grain experimental 2 8 1 grasslands to explore their role in mediating niche-based processes (61) and biodiversity effects (62), crucial mechanism in determining spatial variation in ecosystem functioning at large spatial scales. Rather dependency of diversity-productivity relationships across spatial grains (Table 1). These recent developments in BEF research and macroecology suggest that conceptual integration between these two 2 8 9 disciplines is just beginning (65), yet further efforts to bridge disciplinary gaps are essential to deepen 2 9 0 current understanding of mechanisms that underpin the shifts in diversity-productivity relationships 2 9 1 across spatial scales. To conclude, we show that the relationship between diversity and productivity strengthens toward coarse grains. This result is in line with expectations from both BEF theory, and some (but not all) expectations 2 9 4 from macroecological studies on NPP→S, and highlights the potential of demographic stochasticity to 2 9 5 distort diversity-productivity relationships at fine grains. Moreover, we find similar support for both 2 9 6 directions of diversity-productivity relationships across spatial grains, revealing that biodiversity and 2 9 7 productivity can be both cause and effect. Future research on this relationship needs to move from fine- impacts of anthropogenic biodiversity change on ecosystem function.  county with a probability proportional to 1/sqrt(forest area+1), which will more likely select small rather 3 1 5 than large counties. This was because small counties can be merged to approach the grain of the large grain dataset) and then after it reached 98 merged spatial units (coarse grain dataset) (Fig. 1). Although 3 2 4 the algorithm substantially reduced variation in area within both spatial grains (Fig. S9), it did not 3 2 5 eliminate the variation entirely, and thus we still used area as a covariate in the statistical analyses at the  Species richness (S). For all spatial grains, we estimated diversity as species richness (S) because it is the 3 2 8 most commonly used and best understood metric of biodiversity, although other measures of diversity 3 2 9 may be better predictors of net primary productivity (67-69). We extracted S at the fine spatial grain from identified to species level. In total, our final dataset included 344 woody species and 93,771 plots. To Name Resolution Service (TNRS)(74), following the protocol described in (75). We included hybrid 3 4 1 forms but excluded any names that could not be resolved to the species level.

4 2
Filtering of species occurrences. We restricted our analyses to woody species occurring in forest. To this 3 4 3 end, we initially filtered the BONAP data to species classified as 'trees' in BONAP's taxonomic query shrubs and subshrubs (77), except for 37 species without such data for which we instead inferred 3 4 7 woodiness from online searches or assumed resemblance among congeneric species. We also filtered out unlikely to be forest occurrences, as inferred from independent species occurrences within forested pixels Supplementary Note). To make species richness data internally consistent across the different spatial 3 5 3 grains, we added a further 6,593 quality-vetted county-level forest occurrences of woody species from 3 5 4 FIA plot records to the 282,991 occurrences in the taxonomically harmonized BONAP dataset. productivity. This was measured as the change in tree C over time due to growth (gC m -2 y -1 ), and is the 3 5 8 sum of aboveground C increment of living trees between two measurements and conservatively excludes recruits and dead trees (67). Tree-level carbon was estimated by multiplying tree-level biomass (see 3 6 0 below) by 0.48, but recognize that gymnosperms may have higher carbon content than that of US. Due to insufficient data on species' dispersal abilities, we assumed that dispersal probability between 3 9 9 focal units and species' occurrences would decay with great-circle distance between the respective 4 0 0 regions' centroids. We explored five alternative exponential distance-decay functions, with scaling  Running title: Cross-scale diversity-productivity Craven et al. 20 All of the variables used in our analyses are listed and summarized in Table 2.  other parts are environmentally unique and small. We employed stratified random sampling (87) for the (2) prevent excessive statistical leverage of the large number of data points from homogeneous areas and 4 2 2 Running title: Cross-scale diversity-productivity Craven et al. 21 (3) reduce spatial pseudoreplication (autocorrelation) by increasing the geographic distance between data 4 2 3 points. We first identified 11 strata at the fine and intermediate grains respectively, using multivariate 4 2 4 regression trees with S, NPP and biomass as response variables and all covariates as predictors (Fig. 1).

2 5
We then took a random and proportionally sized sample of spatial units from each strata (fine grain, N = 4 2 6 1,000; intermediate grain, N = 500). We did not use stratified random sampling at the coarse spatial grain 4 2 7 because the number of spatial units was small (N = 98) and spatial autocorrelation was low. The spatial 4 2 8 locations of the stratified samples are in Fig. S1. All of the analyses presented here, as well as our main 4 2 9 conclusions, are based on these stratified sub-samples of the data. within models, we fitted Random Forest models (RFs) (33). The results from SEMs provide insight into 4 3 9 differences among models (i.e. between the two causal pathways per spatial grain, and among spatial covariates on S, NPP, and biomass (except for area at the fine spatial grain) ( Figure S2). An effect of area 4 4 6 Model fit can only be tested on unsaturated models, i.e. those that have at least one missing path.

5 9
Therefore, we removed the path with the lowest standardized path coefficient from the model. As SEMs 4 6 0 had equal number of paths, we could compare model fit across all models within each spatial grain using 4 6 1 their unadjusted R 2 values. After excluding the additional paths, path coefficients of S, NPP, and biomass 4 6 2 remained qualitatively the same, and model fit to the data were still accepted (Chi-square test; P>0.05).

6 3
This indicates that the models are identifiable and their results are robust. Therefore, we did not further 4 6 4 reduce the model, and models maintained the same number of paths within each scale.

6 5
To assess the differences among scales in the relationships between S, NPP and biomass for each model, we compared the standardized regression coefficients using their 95% confidence intervals. All SEMs 4 6 7 were fitted using the 'sem' function of the 'lavaan' package in R (88).