Research ArticleFOREST ECOLOGY

Positive biodiversity-productivity relationship predominant in global forests

See allHide authors and affiliations

Science  14 Oct 2016:
Vol. 354, Issue 6309, aaf8957
DOI: 10.1126/science.aaf8957

Global biodiversity and productivity

The relationship between biodiversity and ecosystem productivity has been explored in detail in herbaceous vegetation, but patterns in forests are far less well understood. Liang et al. have amassed a global forest data set from >770,000 sample plots in 44 countries. A positive and consistent relationship can be discerned between tree diversity and ecosystem productivity at landscape, country, and ecoregion scales. On average, a 10% loss in biodiversity leads to a 3% loss in productivity. This means that the economic value of maintaining biodiversity for the sake of global forest productivity is more than fivefold greater than global conservation costs.

Science, this issue p. 196

Structured Abstract

INTRODUCTION

The biodiversity-productivity relationship (BPR; the effect of biodiversity on ecosystem productivity) is foundational to our understanding of the global extinction crisis and its impacts on the functioning of natural ecosystems. The BPR has been a prominent research topic within ecology in recent decades, but it is only recently that we have begun to develop a global perspective.

RATIONALE

Forests are the most important global repositories of terrestrial biodiversity, but deforestation, forest degradation, climate change, and other factors are threatening approximately one half of tree species worldwide. Although there have been substantial efforts to strengthen the preservation and sustainable use of forest biodiversity throughout the globe, the consequences of this diversity loss pose a major uncertainty for ongoing international forest management and conservation efforts. The forest BPR represents a critical missing link for accurate valuation of global biodiversity and successful integration of biological conservation and socioeconomic development. Until now, there have been limited tree-based diversity experiments, and the forest BPR has only been explored within regional-scale observational studies. Thus, the strength and spatial variability of this relationship remains unexplored at a global scale.

RESULTS

We explored the effect of tree species richness on tree volume productivity at the global scale using repeated forest inventories from 777,126 permanent sample plots in 44 countries containing more than 30 million trees from 8737 species spanning most of the global terrestrial biomes. Our findings reveal a consistent positive concave-down effect of biodiversity on forest productivity across the world, showing that a continued biodiversity loss would result in an accelerating decline in forest productivity worldwide.

The BPR shows considerable geospatial variation across the world. The same percentage of biodiversity loss would lead to a greater relative (that is, percentage) productivity decline in the boreal forests of North America, Northeastern Europe, Central Siberia, East Asia, and scattered regions of South-central Africa and South-central Asia. In the Amazon, West and Southeastern Africa, Southern China, Myanmar, Nepal, and the Malay Archipelago, however, the same percentage of biodiversity loss would lead to greater absolute productivity decline.

CONCLUSION

Our findings highlight the negative effect of biodiversity loss on forest productivity and the potential benefits from the transition of monocultures to mixed-species stands in forestry practices. The BPR we discover across forest ecosystems worldwide corresponds well with recent theoretical advances, as well as with experimental and observational studies on forest and nonforest ecosystems. On the basis of this relationship, the ongoing species loss in forest ecosystems worldwide could substantially reduce forest productivity and thereby forest carbon absorption rate to compromise the global forest carbon sink. We further estimate that the economic value of biodiversity in maintaining commercial forest productivity alone is $166 billion to$490 billion per year. Although representing only a small percentage of the total value of biodiversity, this value is two to six times as much as it would cost to effectively implement conservation globally. These results highlight the necessity to reassess biodiversity valuation and the potential benefits of integrating and promoting biological conservation in forest resource management and forestry practices worldwide.

Abstract

The biodiversity-productivity relationship (BPR) is foundational to our understanding of the global extinction crisis and its impacts on ecosystem functioning. Understanding BPR is critical for the accurate valuation and effective conservation of biodiversity. Using ground-sourced data from 777,126 permanent plots, spanning 44 countries and most terrestrial biomes, we reveal a globally consistent positive concave-down BPR, showing that continued biodiversity loss would result in an accelerating decline in forest productivity worldwide. The value of biodiversity in maintaining commercial forest productivity alone—US$166 billion to 490 billion per year according to our estimation—is more than twice what it would cost to implement effective global conservation. This highlights the need for a worldwide reassessment of biodiversity values, forest management strategies, and conservation priorities. The biodiversity-productivity relationship (BPR) has been a major ecological research focus over recent decades. The need to understand this relationship is becoming increasingly urgent in light of the global extinction crisis because species loss affects the functioning and services of natural ecosystems (1, 2). In response to an emerging body of evidence that suggests that the functioning of natural ecosystems may be substantially impaired by reductions in species richness (310), global environmental authorities, including the Intergovernmental Platform on Biodiversity and Ecosystem Services (IPBES) and United Nations Environment Programme (UNEP), have made substantial efforts to strengthen the preservation and sustainable use of biodiversity (2, 11). Successful international collaboration, however, requires a systematic assessment of the value of biodiversity (11). Quantification of the global BPR is thus urgently needed to facilitate the accurate valuation of biodiversity (12), the forecast of future changes in ecosystem services worldwide (11), and the integration of biological conservation into international socioeconomic development strategies (13). The evidence of a positive BPR stems primarily from studies of herbaceous plant communities (14). In contrast, the forest BPR has only been explored at the regional scale [(3, 4, 7, 15) and references therein] or within a limited number of tree-based experiments [(16, 17) and references therein], and it remains unclear whether these relationships hold across forest types. Forests are the most important global repositories of terrestrial biodiversity (18), but deforestation, climate change, and other factors are threatening a considerable proportion (up to 50%) of tree species worldwide (1921). The consequences of this diversity loss pose a critical uncertainty for ongoing international forest management and conservation efforts. Conversely, forest management that converts monocultures to mixed-species stands has often seen a substantial positive effect on productivity with other benefits (2224). Although forest plantations are predicted to meet 50 to 75% of the demand for lumber by 2050 (25, 26), nearly all are still planted as monocultures, highlighting the potential of forest management in strengthening the conservation and sustainable use of biodiversity worldwide. Here, we compiled in situ remeasurement data, most of which were taken at two consecutive inventories from the same localities, from 777,126 permanent sample plots [hereafter, global forest biodiversity (GFB) plots] across 44 countries and territories and 13 ecoregions to explore the forest BPR at a global scale (Fig. 1). GFB plots encompass forests of various origins (from naturally regenerated to planted) and successional stages (from stand initiation to old-growth). A total of more than 30 million trees across 8737 species were tallied and measured on two or more consecutive inventories from the GFB plots. Sampling intensity was greater in developed countries, where nationwide forest inventories have been fully or partially funded by governments. In most other countries, national forest inventories were lacking, and most ground-sourced data were collected by individuals and organizations (table S1). On the basis of ground-sourced GFB data, we quantified BPR at the global scale using a data-driven ensemble learning approach (Materials and methods, Geospatial random forest). Our quantification of BPR involved characterizing the shape and strength of the dependency function through the elasticity of substitution (θ), which represents the degree to which species can substitute for each other in contributing to forest productivity; θ measures the marginal productivity—the change in productivity resulting from one unit decline of species richness—and reflects the strength of the effect of tree diversity on forest productivity, after accounting for climatic, soil, and plot-specific covariates. A higher θ corresponds to a greater decline in productivity due to one unit loss in biodiversity. The niche-efficiency (N-E) model (3) and several preceding studies (2730) provide a framework for interpreting the elasticity of substitution and approximating BPR with a power function model:P = α · f(X) · Sθ (1)where P and S signify primary site productivity and tree species richness (observed on a 900-m2 area basis on average) (Materials and methods), respectively; f(X) is a function of a vector of control variables X (selected from stand basal area and 14 climatic, soil, and topographic covariates); and α is a constant. This model is capable of representing a variety of potential patterns of BPR. 0 < θ < 1 represents a positive and concave down pattern (a degressively increasing curve), which is consistent with the N-E model and preceding studies (3, 2730), whereas other θ values can represent alternative BPR patterns, including decreasing (θ < 0), linear (θ = 1), convex (θ > 1), or no effect (θ = 0) (Fig. 2) (14, 31). The model (Eq. 1) was estimated by using the geospatial random forest technique based on GFB data and covariates acquired from ground-measured and remote-sensing data (Materials and methods). We found that a positive BPR predominated in forests worldwide. Out of 10,000 randomly selected subsamples (each consisting of 500 GFB plots), 99.87% had a positive concave-down relationship with relative species richness (0 < θ < 1), whereas only 0.13% show negative trends, and none was equal to zero or greater than or equal to 1 (Fig. 2). Overall, the global forest productivity increased with a declining rate from 2.7 to 11.8 m3 ha−1 year−1 as relative tree species richness increased from the minimum to the maximum value, which corresponds to a θ value of 0.26 (Fig. 3A). At the global scale, we mapped the magnitude of BPR (as expressed by θ) using geospatial random forest and universal kriging. By plotting values of θ onto a global map, we revealed considerable geospatial variation across the world (Fig. 3B). The highest θ (0.29 to 0.30) occurred in the boreal forests of North America, Northeastern Europe, Central Siberia, and East Asia and the sporadic tropical and subtropical forests of South-central Africa, South-central Asia, and the Malay Archipelago. In these areas of the highest elasticity of substitution (32), the same percentage of biodiversity loss would lead to a greater percentage of reduction in forest productivity (Fig. 4A). In terms of absolute productivity, the same percentage of biodiversity loss would lead to the greatest productivity decline in the Amazon; West Africa’s Gulf of Guinea; Southeastern Africa, including Madagascar; Southern China; Myanmar; Nepal; and the Malay Archipelago (Fig. 4B). Because of a relatively narrow range of the elasticity of substitution (32) estimated from the global-level analysis (0.2 to 0.3), the regions of the greatest productivity decline under the same percentage of biodiversity loss largely matched the regions of the greatest productivity (fig. S1). Globally, a 10% decrease in tree species richness (from 100 to 90%) would cause a 2 to 3% decline in productivity, and with a decrease in tree species richness to one (Materials and methods, Economic analysis), this decline in forest productivity would be 26 to 66% even if other things, such as the total number of trees and forest stocking, remained the same (fig. S4). Discussion Our global analysis provides strong and consistent evidence that productivity of forests would decrease at an accelerating rate with the loss of biodiversity. The positive concave-down pattern we discovered across forest ecosystems worldwide corresponds well with recent theoretical advances in BPR (3, 2830), as well as with experimental (27) and observational (14) studies on forest and nonforest ecosystems. The elasticity of substitution (32) estimated in this study (ranged between 0.2 and 0.3) largely overlaps the range of values of the same exponent term (0.1 to 0.5) from previous theoretical and experimental studies [(10) and references therein]. Furthermore, our findings are consistent with the global estimates of the biodiversity-dependent ecosystem service debt under distinct assumptions (10) and with recent reports of the diminishing marginal benefits of adding a species as species richness increases, based on long-term forest experiments dating back to 1870 [(15, 33) and references therein]. Our analysis relied on stands ranging from unmanaged to extensively managed forests—managed forests with low operating and investment costs per unit area. Conditions of natural forests would not be comparable with intensively managed forests, because timber production in the latter systems often focuses on a single or limited number of highly productive tree species. Intensively managed forests, where saturated resources can weaken the effects of niche efficiency (3), are shown in some studies (34, 35) to have higher productivity than that of natural diverse forests of the same climate and site conditions (fig. S3). In contrast, other studies (6, 2224) compared diverse stands with monocultures at the same level of management intensity and found that the positive effects of species diversity on tree productivity and other ecosystem services are applicable to intensively managed forests. As such, there is still an unresolved debate on the BPR of intensively managed forests. Nevertheless, because intensively managed forests only account for a minor (<7%) portion of global forests (18), our estimated BPR would be minimally affected by such manipulations and thus should reflect the inherent processes governing the vast majority of global forest ecosystems. We focused on the effect of biodiversity on ecosystem productivity. Recent studies on the opposite causal direction [productivity-biodiversity relationship (14, 36, 37)] suggest that there may be a potential two-way causality between biodiversity and productivity. It is admittedly difficult to use correlative data to detect and attribute causal effects. Fortunately, substantial progress has been made to tease the BPR causal relationship from other potentially confounding environmental variables (14, 38, 39), and this study made considerable efforts to account for these otherwise potentially confounding environmental covariates in assessing likely causal effects of biodiversity on productivity. Because taxonomic diversity indirectly incorporates functional, phylogenetic, and genomic diversity, our results that focus on tree species richness are likely applicable to these other elements of biodiversity, all of which have been found to influence plant productivity (1). Our straightforward analysis makes clear the taxonomic contribution to forest ecosystem productivity and functioning, and the importance of preserving species diversity to biological conservation and forest management. Our findings highlight the necessity to reassess biodiversity valuation and reevaluate forest management strategies and conservation priorities in forests worldwide. In terms of global carbon cycle and climate change, the value of biodiversity may be considerable. On the basis of our global-scale analyses (Fig. 4), the ongoing species loss in forest ecosystems worldwide (1, 21) could substantially reduce forest productivity and thereby forest carbon absorption rate, which would in turn compromise the global forest carbon sink (40). We further estimate that the economic value of biodiversity in maintaining commercial forest productivity is$166 billion to $490 billion per year (1.66 × 1011 to 4.90 × 1011 year−1 in 2015 US$) (Materials and methods, Economics analysis). By itself, this estimate does not account for other values of forest biodiversity (including potential values for climate regulation, habitat, water flow regulation, and genetic resources), and represents only a small percentage of the total value of biodiversity (41, 42). However, this value is already between two to six times the total estimated cost that would be necessary if we were to effectively conserve all terrestrial ecosystems at a global scale [$76.1 billion per year (43)]. The high benefit-to-cost ratio underlines the importance of conserving biodiversity for forestry and forest resource management. Amid the struggle to combat biodiversity loss, the relationship between biological conservation and poverty is gaining increasing global attention (13, 44), especially with respect to rural areas where livelihoods depend most directly on ecosystem products. Given the substantial geographic overlaps between severe, multifaceted poverty and key areas of global biodiversity (45), the loss of species in these areas has the potential to exacerbate local poverty by diminishing forest productivity and related ecosystem services (44). For example, in tropical and subtropical regions, many areas of high elasticity of substitution (32) overlapped with biodiversity hotspots (46), including Eastern Himalaya and Nepal, Mountains of Southwest China, Eastern Afromontane, Madrean pine-oak woodlands, Tropical Andes, and Cerrado. For these areas, only a few species of commercial value are targeted by logging. As such, the risk of losing species through deforestation would far exceed the risk through harvesting (47). Deforestation and other anthropogenic drivers of biodiversity loss in these biodiversity hotspots are likely to have considerable impacts on the productivity of forest ecosystems, with the potential to exacerbate local poverty. Furthermore, the greater uncertainty in our results for the developing countries (Fig. 5) reflects the well-documented geographic bias in forest sampling, including repeated measurements, and reiterates the need for strong commitments toward improving sampling in the poorest regions of the world. Our findings reflect the combined strength of large-scale integration and synthesis of ecological data and modern machine learning methods to increase our understanding of the global forest system. Such approaches are essential for generating global insights into the consequences of biodiversity loss and the potential benefits of integrating and promoting biological conservation in forest resource management and forestry practices—a common goal already shared by intergovernmental organizations such as the Montréal and Helsinki Process Working Groups. These findings should facilitate efforts to accurately forecast future changes in ecosystem services worldwide, which is a primary goal of IPBES (11), and provide baseline information necessary to establish international conservation objectives, including the United Nations Convention on Biological Diversity Aichi targets, the United Nations Framework Convention on Climate Change REDD+ goal, and the United Nations Convention to Combat Desertification land degradation neutrality goal. The success of these goals relies on the understanding of the intrinsic link between biodiversity and forest productivity. Materials and methods Data collection and standardization Our current study used ground-sourced forest measurement data from 45 forest inventories collected from 44 countries and territories (Fig. 1 and table S1). The measurements were collected in the field from predesignated sample area units, i.e., Global Forest Biodiversity permanent sample plots (hereafter, GFB plots). For the calculation of primary site productivity, GFB plots can be categorized into two tiers. Plots designated as “Tier 1” have been measured at two or more points in time with a minimum time interval between measurements of two years or more (global mean time interval is 9 years, see Table 1). “Tier 2” plots were only measured once, and primary site productivity can be estimated from known stand age or dendrochronological records. Overall, our study was based on 777,126 GFB plots, of which 597,179 (77%) were Tier 1, and 179,798 (23%) were Tier 2. GFB plots primarily measured natural forests ranging from unmanaged to extensively managed forests, i.e., managed forests with low operating and investment costs per unit area. Intensively managed forests with harvests exceeding 50 percent of the stocking volume were excluded from this study. GFB plots represent forests of various origins (from naturally regenerated to planted) and successional stages (from stand initiation to old-growth). Table 1 Definition, unit, and summary statistics of key variables. View this table: For each GFB plot, we derived three key attributes from measurements of individual trees—tree species richness (S), stand basal area (G), and primary site productivity (P). Because for each of all the GFB plot samples, S and P were derived from the measurements of the same trees, the sampling issues commonly associated with biodiversity estimation (48) had little influence on the SP relationship (i.e., BPR) in this study. Species richness, S, represents the number of different tree species alive at the time of inventory within the perimeter of a GFB plot with an average size of approximately 900 m2. Ninety-five percent of all plots fall between 100 and 1,100 m2 in size. To minimize the species-area effect (49), we studied the BPR here using a geospatial random forest model in which observations from nearby GFB plots would be more influential than plots that are farther apart (see §Geospatial random forest). Because nearby plots are most likely from the same forest inventory data set, and there was no or little variation of plot area within each data set, the BPR derived from this model largely reflected patterns under the same plot area basis. To investigate the potential effects of plot size on our results, we plotted the estimated elasticity of substitution (θ) against plot size, and found that the scatter plot was normally distributed with no discernible pattern (fig. S2). In addition, the fact that the plot size indicator I2 had the second lowest (0.8%) importance score (50) among all the covariates (Fig. 6) further supports that the influence of plot size variation in this study was negligible. Across all the GFB plots, there were 8,737 species in 1,862 genera and 231 families, and S values ranged from 1 to 405 per plot. We verified all the species names against 60 taxonomic databases, including NCBI, GRIN Taxonomy for Plants, Tropicos–Missouri Botanical Garden, and the International Plant Names Index, using the ‘taxize’ package in R (51). Out of 8737 species recorded in the GFB database, 7425 had verified taxonomic information with a matching score (51) of 0.988 or higher, whereas 1312 species names partially matched existing taxonomic databases with a matching score between 0.50 and 0.75, indicating that these species may have not been documented in the 60 taxonomic databases. To facilitate inter-biome comparison, we further developed relative species richness (Š), a continuous percentage score converted from species richness (S) and the maximal species richness of a set of sample plots(S*) using (2)Stand basal area (G, in m2 ha−1) represents the total cross-sectional area of live trees per unit sample area. G was calculated from individual tree diameter-at-breast-height (dbh, in cm): (3)where κi denotes the conversion factor (ha−1) of the ith tree, viz. the number of trees per ha represented by that individual. G is a key biotic factor of forest productivity as it represents stand density—often used as a surrogate for resource acquisition (through leaf area) and stand competition (52). Accounting for basal area as a covariate mitigated the artifact of different minimum dbh across inventories, and the artifact of different plot sizes. Primary site productivity (P, in m3 ha−1 yr−1) was measured as tree volume productivity in terms of periodic annual increment (PAI) calculated from the sum of individual tree stem volume (V, in m3) (4)where Vi,1 and Vi,2 (in m3) represent total stem volume of the ith tree at the time of the first inventory and the second inventory, respectively. M denotes total removal of trees (including mortality, harvest, and thinning) in stem volume (in m3 ha−1). Y represents the time interval (in years) between two consecutive inventories. P accounted for mortality, ingrowth (i.e., recruitment between two inventories), and volume growth. Stem volume values were predominantly calculated using region- and species-specific allometric equations based on dbh and other tree- and plot-level attributes (Table 1). For the regions lacking an allometric equation, we approximated stem volume at the stand level from basal area, total tree height, and stand form factors (53). In case of missing tree height values from the ground measurement, we acquired alternative measures from a global 1-km forest canopy height database (54). For Tier 2 plots that lacked remeasurement, P was measured in mean annual increment (MAI) based on total stand volume and stand age (52), or tree radial growth measured from increment cores. Since the traditional MAI metric does not account for mortality, we calculated P by adding to MAI the annual mortality based on regional-specific forest turnover rates (55). The small and insignificant correlation coefficient between P and the indicator of plot tier (I1), together with the negligible variable importance of I1 (1.8%, Fig. 6), indicate that PAI and MAI were generally consistent, such that MAI could be a good proxy of PAI in our study. Although MAI and PAI have considerable uncertainty in any given stand, it is difficult to see how systematic bias across diversity gradients could occur on a scale sufficient to influence the results shown here. P, although only representing a fraction of total forest net primary production, has been an important and widely used measure of forest productivity, because it reflects the dominant aboveground biomass component and the long-lived biomass pool in most forest ecosystems (56). Additionally, although other measures of productivity (e.g., net ecosystem exchange processed to derive gross and net primary production; direct measures of aboveground net primary production including all components; and remotely sensed estimates of LAI and greenness coupled with models) all have their advantages and disadvantages, none would be feasible at a similar scale and resolution as in this study. To account for abiotic factors that may influence primary site productivity, we compiled 14 geospatial covariates based on biological relevance and spatial resolution (Fig. 6). These covariates, derived from satellite-based remote sensing and ground-based survey data, can be grouped into three categories: climatic, soil, and topographic (Table 1). We preprocessed all geospatial covariates using ArcMap 10.3 (57) and R 2.15.3 (58). All covariates were extracted to point locations of GFB plots, with a nominal resolution of 1 km2. Geospatial random forest We developed geospatial random forest—a data-driven ensemble learning approach—to characterize the biodiversity–productivity relationship (BPR), and to map BPR in terms of elasticity of substitution (31) on all sample sites across the world. This approach was developed to overcome two major challenges that arose from the size and complexity of GFB data without assuming any underlying BPR patterns or data distribution. First, we need to account for broad-scale differences in vegetation types, but global classification and mapping of homogeneous vegetation types is lacking (59); and secondly, correlations and trends that naturally occur through space (60) can be significant and influential in forest ecosystems (61). Geostatistical models (62) have been developed to address the spatial autocorrelation, but the size of the GFB data set far exceeds the computational constraints of most geostatistical software. Geospatial random forest integrated conventional random forest (50) and a geostatistical nonlinear mixed-effects model (63) to estimate BPR across the world based on GFB plot data and their spatial dependence. The underlying model had the following form (5)where logPij(u) and logŠij(u) represent natural logarithm of productivity and relative species richness (calculated from actual species richness and the maximal species richness of the training set) of plot i in the jth training set at point locations u, respectively. The model was derived from the niche–efficiency model, and θ corresponds to the elasticity of substitution (31). αi·Xij(u)=αi0 + αi1·xij1+…+ αin·xijn represents n covariates and their coefficients (Fig. 6 and Table 1). To account for potential spatial autocorrelation, which can bias tests of significance due to the violation of independence assumption and is especially problematic in large-scale forest ecosystem studies (60, 61), we incorporated a spherical variogram model (62) into the residual term eij(u). The underlying geostatistical assumption was that across the world BPR is a second-order stationary process—a common geographical phenomenon in which neighboring points are more similar to each other than they are to points that are more distant (64). In our study, we found strong evidence for this gradient (Fig. 7), indicating that observations from nearby GFB plots would be more influential than plots that are farther away. The positive spherical semivariance curves estimated from a large number of bootstrapping iterations indicated that spatial dependence increased as plots became closer together. The aforementioned geostatistical nonlinear mixed-effects model was integrated into random forest analysis (50) by means of model selection and estimation. In the model selection process, random forest was employed to assess the contribution of each of the candidate variables to the dependent variable logPij(u), in terms of the amount of increase in prediction error as one variable is permuted while all the others are kept constant. We used the randomForest package (65) in R to obtain importance measures for all the covariates to guide our selection of the final variables in the geostatistical nonlinear mixed-effects model, Xij(u). We selected stand basal area (G), temperature seasonality (T3), annual precipitation (C1), precipitation of the warmest quarter (C3), potential evapotranspiration (PET), indexed annual aridity (IAA), and plot elevation (E) as control variables since their importance measures were greater than the 9 percent threshold (Fig. 6) preset to ensure that the final variables accounted for over 60 percent of the total variable importance measures. For geospatial random forest analysis of BPR, we first selected control variables based on the variable importance measures derived from random forests (50). We then evaluated the values of elasticity of substitution (32), which are expected to be real numbers greater than 0 and less than 1, against the alternatives, i.e., negative BPR (H01: θ < 0), no effect (H02: θ = 0), linear (H03: θ = 1), and convex positive BPR (H04: θ > 1). We examined all the coefficients by their statistical significance and effect sizes, using Akaike information criterion (AIC), Bayesian information criterion (BIC), and the generalized coefficient of determination (66). Global analysis For the global-scale analysis, we calibrated the nonlinear mixed-effects model parameters (θ and α’s) using training sets of 500 plots randomly selected (with replacement) from the GFB global dataset according to the bootstrap aggregating (bagging) algorithm. We calibrated a total of 10,000 models based on the bagging samples, using our own bootstrapping program and the nonlinear package nlme (63) of R, to calculate the means and standard errors of final model estimates (Table 2). This approach overcame computational limits by partitioning the GFB sample into smaller subsamples to enable the nonlinear estimation. The size of training sets was selected based on the convergence and effect size of the geospatial random forest models. In pilot simulations with increasing sizes of training sets (Fig. 8), the value of elasticity of substitution (32) fluctuated at the start until the convergence point at 500 plots. Generalized R2 values declined as the size of training sets increased from 0 to 350 plots, and stabilized at around 0.35 as training set size increased further. Accordingly, we selected 500 as the size of the training sets for the final geospatial random forest analysis. Based on the estimated parameters of the global model (Table 2), we analyzed the effect of relative species richness on global forest productivity with a sensitivity analysis by keeping all the other variables constant at their sample means for each ecoregion. Table 2 Parameters of the global geospatial random forest model in 10,000 iterations of 500 randomly selected (with replacement) GFB plots. Mean and SE of all the parameters were estimated by using bootstrapping. Effect sizes were represented by the Akaike information criterion (AIC), Bayesian information criterion (BIC), and generalized R2 (G-R2). Const, constant. View this table: Mapping BPR across global forest ecosystems For mapping purposes, we first estimated the current extent of global forests in several steps. We aggregated the “treecover2000” and “loss” data (67) from 30 m pixels to 30 arc-second pixels (~1 km) by calculating the respective means. The result was ~1 km pixels showing the percentage forest cover for the year 2000 and the percentage of this forest cover lost between 2000 and 2013, respectively. The aggregated forest cover loss was multiplied by the aggregated forest cover to produce a single raster value for each ~1 km pixel representing a percentage forest lost between 2000 and 2013. This multiplication was necessary since the initial loss values were relative to initial forest cover. Similarly, we estimated the percentage forest cover gain by aggregating the forest “gain” data (67) from 30 m to 30 arc-seconds while taking a mean. Then, this gain layer was multiplied by 1 minus the aggregated forest cover from the first step to produce a single value for each ~1 km pixel that signifies percentage forest gain from 2000–2013. This multiplication ensured that the gain could only occur in areas that were not already forested. Finally, the percentage forest cover for 2013 was computed by taking the aggregated data from the first step (year 2000) and subtracting the computed loss and adding the computed gain. We mapped productivity P and elasticity of substitution (32) across the estimated current extent of global forests, here defined as areas with 50 percent or more forest cover. Because GFB ground plots represent approximately 40 percent of the forested areas, we used universal kriging (62) to estimate P and θ for the areas with no GFB sample coverage. The universal kriging models consisted of covariates specified in Fig. 6B and a spherical variogram model with parameters (i.e., nugget, range, and sill) specified in Fig. 7. We obtained the best linear unbiased estimators of P and θ and their standard error in relation to Š across the current global forest extent with the gstat package of R (68). By combining θ estimated from geospatial random forest and universal kriging, we produced the spatially continuous maps of the elasticity of substitution (Fig. 3B) and forest productivity (fig. S1) at a global scale. The effect sizes of the best linear unbiased estimator of θ (in terms of standard error and generalized R2) are shown in Fig. 5. We further estimated percentage and absolute decline in worldwide forest productivity under two scenarios of loss in tree species richness— low (10% loss) and high (99% loss). These levels represent the productivity decline (in both percentage and absolute terms) if local species richness across the global forest extent would decrease to 90 and 1 percent of the current values, respectively. The percentage decline was calculated based on the general BPR model (Eq. 1) and estimated worldwide spatially explicit values of the elasticity of substitution (Fig. 3B). The absolute decline was the product of the worldwide estimates of primary forest productivity (fig. S1) and the standardized percentage decline at the two levels of biodiversity loss (Fig. 4A). Economic analysis Estimates of the economic value-added from forests employ a range of methods. One prominent recent global valuation of ecosystem services (69) valued global forest production [in terms of ‘raw materials’ (including timber, fiber, biomass fuels, and fuelwood and charcoal] provided by forests (table S1) (69) in 2011 at US$ 649 billion (6.49 × 1011, in constant 2007 dollars). Using an alternative method, the UN FAO (25, 26) estimates gross value-added in the formal forestry sector, a measure of the contribution of forestry, wood industry, and pulp and paper industry to the world’s economy, at US$606 billion (6.06 × 1011, in constant 2011 dollars). Because these two reasonably comparable values are directly impacted by and proportional to forest productivity, we used them as bounds on our coarse estimate of the global economic value of commercial forest productivity, converted to constant 2015 US$ based on the US consumer price indices (70, 71). As indicated by our global-scale analyses (Fig. 4A), a 10 percent decrease of tree species richness distributed evenly across the world (from 100% to 90%) would cause a 2.1 to 3.1 percent decline in productivity, which would equate to US$13–23 billion per year (constant 2015 US$). For the assessment of the value of biodiversity in maintaining forest productivity, a drop in species richness from the current level to one species would lead to 26–66% reduction in commercial forest productivity in the biomes that contribute substantially to global commercial forestry (fig.S4), equivalent to 166–490 billion US$per year (1.66 × 1011 to 4.90 × 1011, constant 2015 US$, calculated by multiplying the foregoing economic value-added from FAO and the other study by 26 and 66%, respectively.) Therefore, we estimated that the economic value of biodiversity in maintaining commercial forest productivity worldwide would be 166 billion to 490 billion US\$ per year.

We held the total number of trees, global forest area and stocking, and other factors constant to estimate the value of productivity loss solely due to a decline in tree species richness. As such, these estimates did not include the value of land converted from forest and losses due to associated fauna or flora decline or forest habitat reduction. This estimate only reflects the value of biodiversity in maintaining commercial forest productivity that contributes directly to forestry, wood industry, and pulp and paper industry, and does not account for other values of biodiversity, including potential values for climate regulation, habitat, water flow regulation, genetic resources, etc. The total global value of biodiversity could exceed this estimate by orders of magnitudes (41, 42).

Supplementary Materials

www.sciencemag.org/cgi/content/full/354/6309/aaf8957/DC1

Figs. S1 to S4

Tables S1 to S2

References (79134)

References and Notes

1. Value-added has been adjusted for inflation and is expressed in USD at 2011 prices and exchange rates.
2. The elasticity of substitution (θ), which represents the degree to which species can substitute for each other in contributing to stand productivity, reflects the strength of the effect of biodiversity on ecosystem productivity, after accounting for climatic, soil, and other environmental and local covariates.
3. We used the online CPI calculator; http://data.bls.gov/cgi-bin/cpicalc.pl.
4. The tree drawings in Fig. 2 were based on actual species from the GFB plots. The scientific names of these species are (clockwise from the top) Abies nebrodensis, Handroanthus albus, Araucaria angustifolia, Magnolia sinica, Cupressus sempervirens, Salix babylonica, Liriodendron tulipifera, Adansonia grandidieri, Torreya taxifolia, and Quercus mongolica. Five of these 10 species (A. nebrodensis, A. angustifolia, M. sinica, A. grandidieri, and T. taxifolia) are listed as endagered or critically endangered species in the IUCN Red List. Hand drawings were made by R. K. Watson.
5. Elevation consists mostly of ground-measured data, and the missing values were replaced with the highest-resolution topographic data generated from NASA’s Shuttle Radar Topography Mission (SRTM).
6. Calculated from stand basal area, canopy height, and stand form factor.