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.

Global effect of tree species diversity on forest productivity.

Ground-sourced data from 777,126 global forest biodiversity permanent sample plots (dark blue dots, left), which cover a substantial portion of the global forest extent (white), reveal a consistent positive and concave-down biodiversity-productivity relationship across forests worldwide (red line with pink bands representing 95% confidence interval, right).

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).

Fig. 1 GFB ground-sourced data were collected from in situ remeasurement of 777,126 permanent sample plots consisting of more than 30 million trees across 8737 species.

GFB plots extend across 13 ecoregions [vertical axis, delineated by the World Wildlife Fund where extensive forests occur within all the ecoregions (72)], and 44 countries and territories. Ecoregions are named for their dominant vegetation types, but all contain some forested areas. GFB plots cover a substantial portion of the global forest extent (white), including some of the most distinct forest conditions: (a) the northernmost (73°N, Central Siberia, Russia), (b) southernmost (52°S, Patagonia, Argentina), (c) coldest (–17°C annual mean temperature, Oimyakon, Russia), (d) warmest (28°C annual mean temperature, Palau, United States), and (e) most diverse (405 tree species on the 1-ha plot, Bahia, Brazil). Plots in war-torn regions [such as (f)] were assigned fuzzed coordinates to protect the identity of the plots and collaborators. The box plots show the mean and interquartile range of tree species richness and primary site productivity (both on a common logarithmic scale) derived from ground-measured tree- and plot-level records. The complete list of species is presented in table S2.

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).

Fig. 2 Theoretical positive and concave-down biodiversity–productivity relationship supported by empirical evidence drawn from the GFB data.

(Left) The diagram demonstrates that under the theoretical positive and concave-down (monotonically and degressively increasing) BPR (3, 27, 28), loss in tree species richness may reduce forest productivity (73). (Middle) Functional curves represent different BPR under different values of elasticity of substitution (θ). θ values between 0 and 1 correspond to the positive and concave-down BPR (blue curve). (Right) The three-dimensional scatter plot shows θ values we estimated from observed productivity (P), species richness (S), and other covariates. Out of 5,000,000 estimates of θ (mean = 0.26, SD = 0.09), 4,993,500 fell between 0 and 1 (blue), whereas only 6500 were negative (red), and none was equal to zero or greater than or equal to 1; the positive and concave-down BPR was supported by 99.9% of our estimates.

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).

Fig. 3 The estimated global effect of biodiversity on forest productivity was positive and concave-down, and revealed considerable geospatial variation across forest ecosystems worldwide.

(A) Global effect of biodiversity on forest productivity (red line with pink bands representing 95% confidence interval) corresponds to a global average elasticity of substitution (θ) value of 0.26, with climatic, soil, and other plot covariates being accounted for and kept constant at sample mean. Relative species richness (Š) is in the horizontal axis, and productivity (P, m3 ha−1 year−1) is in the vertical axis (histograms of the two variables on top and right in the logarithm scale). (B) θ represents the strength of the effect of tree diversity on forest productivity. Spatially explicit values of θ were estimated by using universal kriging (Materials and methods) across the current global forest extent (effect sizes of the estimates are shown in Fig. 5), whereas blank terrestrial areas were nonforested.

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).

Fig. 4 Estimated percentage and absolute decline in forest productivity under 10 and 99% decline in current tree species richness (values in parentheses correspond to 99%), holding all the other terms constant.

(A) Percent decline in productivity was calculated according to the general BPR model (Eq. 1) and estimated worldwide spatially explicit values of the elasticity of substitution (Fig. 3B). (B) Absolute decline in productivity was derived from the estimated elasticity of substitution (Fig. 3B) and estimates of global forest productivity (fig. S1). The first 10% reduction in tree species richness would lead to a 0.001 to 0.597 m3 ha−1 year−1 decline in periodic annual increment, which accounts for 2 to 3% of current forest productivity. The raster data are displayed in 50-km resolution with a 3 SD stretch.

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.

Fig. 5 Standard error and generalized R2 of the spatially explicit estimates of elasticity of substitution (θ) across the current global forest extent in relation to Š.

(A) Standard error increased as a location was farther from those sampled. (B) The generalized R2 values were derived with a geostatistical nonlinear mixed-effects model for GFB sample locations, and thus (B) only covers a subset of the current global forest extent.

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.

Fig. 6 Correlation matrix and importance values of potential variables for the geospatial random forest analysis.

(A)There were a total of 15 candidate variables from three categories, namely plot attributes, climatic variables, and soil factors (a detailed description is provided in Table 1). Correlation coefficients between these variables were represented by sizes and colors of circles, and “×” marks coefficients not significant at α = 0.05 level. (B) Variable importance (%) values were determined by the geospatial random forest (Materials and methods). Variables with importance values exceeding the 9% threshold line (blue) were selected as control variables in the final geospatial random forest models. Elasticity of substitution (coefficient), productivity (dependent variable), and species richness (key explanatory variable) were not ranked in the variable importance chart because they were not potential covariates.

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*) usingEmbedded Image (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):Embedded Image (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)Embedded Image (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 formEmbedded Image (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.

Fig. 7 Semivariance and estimated spherical variogram models (blue curves) obtained from geospatial random forest in relation to Š.

Gray circles, semivariance; blue curves, estimated spherical variogram models. There was a general trend that semivariance increased with distance; spatial dependence of θ weakened as the distance between any two GFB plots increased. The final spherical models had nugget = 0.8, range = 50 degrees, and sill = 1.3. To avoid identical distances, all plot coordinates were jittered by adding normally distributed random noises.

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:
Fig. 8 Effect of the size of training sets used in the geospatial random forest on estimated elasticity of substitution (θ) and generalized R2 in relation to Š.

Mean (solid line) and SE band (green area) were estimated with 100 randomly selected (with replacement) training sets for each of the 20 size values (between 50 and 1000 GFB plots, with an increment of 50).

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.
Acknowledgments: We are grateful to all the people and agencies that helped in collection, compilation, and coordination of the field data, including but not limited to T. Malone, J. Crowe, M. Sutton, J. Lovett, P. Munishi, M. Rautiainen, staff members from the Seoul National University Forest, and all persons who made the two Spanish Forest Inventories possible, especially the main coordinators, R. Villaescusa (IFN2) and J. A. Villanueva (IFN3). This work was supported in part by West Virginia University under the United States Department of Agriculture (USDA) McIntire-Stennis Funds WVA00104 and WVA00105; U.S. National Science Foundation (NSF) Long-Term Ecological Research Program at Cedar Creek (DEB-1234162); the University of Minnesota Department of Forest Resources and Institute on the Environment; the Architecture and Environment Department of Italcementi Group, Bergamo (Italy); a Marie Skłodowska Curie fellowship; Polish National Science Center grant 2011/02/A/NZ9/00108; the French L’Agence Nationale de la Recherche (ANR) (Centre d’Étude de la Biodiversité Amazonienne: ANR-10-LABX-0025); the General Directory of State Forest National Holding DB; General Directorate of State Forests, Warsaw, Poland (Research Projects 1/07 and OR/2717/3/11); the 12th Five-Year Science and Technology Support Project (grant 2012BAD22B02) of China; the U.S. Geological Survey and the Bonanza Creek Long Term Ecological Research Program funded by NSF and the U.S. Forest Service (any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government); National Research Foundation of Korea (grant NRF-2015R1C1A1A02037721), Korea Forest Service (grants S111215L020110, S211315L020120 and S111415L080120) and Promising-Pioneering Researcher Program through Seoul National University (SNU) in 2015; Core funding for Crown Research Institutes from the New Zealand Ministry of Business, Innovation and Employment’s Science and Innovation Group; the Deutsche Forschungsgemeinschaft (DFG) Priority Program 1374 Biodiversity Exploratories; Chilean research grants Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) 1151495 and 11110270; Natural Sciences and Engineering Research Council of Canada (grant RGPIN-2014-04181); Brazilian Research grants CNPq 312075/2013 and FAPESC 2013/TR441 supporting Santa Catarina State Forest Inventory (IFFSC); the General Directorate of State Forests, Warsaw, Poland; the Bavarian State Ministry for Nutrition, Agriculture, and Forestry project W07; the Bavarian State Forest Enterprise (Bayerische Staatsforsten AöR); German Science Foundation for project PR 292/12-1; the European Union for funding the COST Action FP1206 EuMIXFOR; FEDER/COMPETE/POCI under Project POCI-01-0145-FEDER-006958 and FCT–Portuguese Foundation for Science and Technology under the project UID/AGR/04033/2013; Swiss National Science Foundation grant 310030B_147092; the EU H2020 PEGASUS project (no 633814), EU H2020 Simwood project (no 613762); and the European Union’s Horizon 2020 research and innovation program within the framework of the MultiFUNGtionality Marie Skłodowska-Curie Individual Fellowship (IF-EF) under grant agreement 655815. The expeditions in Cameroon to collect the data were partly funded by a grant from the Royal Society and the Natural Environment Research Council (UK) to Simon L. Lewis. Pontifica Universidad Católica del Ecuador offered working facilities and reduced station fees to implement the census protocol in Yasuni National Park. We thank the following agencies and organization for providing the data: USDA Forest Service; School of Natural Resources and Agricultural Sciences, University of Alaska Fairbanks; the Ministère des Forêts, de la Faune et des Parcs du Québec (Canada); the Alberta Department of Agriculture and Forestry, the Saskatchewan Ministry of the Environment, and Manitoba Conservation and Water Stewardship (Canada); the National Vegetation Survey Databank (New Zealand); Italian and Friuli Venezia Giulia Forest Services (Italy); Bavarian State Forest Enterprise (Bayerische Staatsforsten AöR) and the Thünen Institute of Forest Ecosystems (Germany); Queensland Herbarium (Australia); Forestry Commission of New South Wales (Australia); Instituto de Conservação da Natureza e das Florestas (Portugal). M'Baïki data were made possible and provided by the ARF Project (Appui la Recherche Forestière) and its partners: AFD (Agence Française de Développement), CIRAD (Centre de Coopération Internationale en Recherche Agronomique pour le Développement), ICRA (Institut Centrafricain de Recherche Agronomique), MEDDEFCP (Ministère de l'Environnement, du Développement Durable des Eaux, Forêts, Chasse et Pêche), SCAC/MAE (Service de Coopération et d’Actions Culturelles, Ministère des Affaires Etrangères), SCAD (Société Centrafricaine de Déroulage), and the University of Bangui. All TEAM data were provided by the Tropical Ecology Assessment and Monitoring (TEAM) Network—a collaboration between Conservation International, the Smithsonian Institute, and the Wildlife Conservation Society—and partially funded by these institutions: the Gordon and Betty Moore Foundation, the Valuing the Arc Project (Leverhulme Trust), and other donors. The Exploratory plots of FunDivEUROPE received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement 265171. The Chinese Comparative Study Plots (CSPs) were established in the framework of BEF-China, funded by the German Research Foundation (DFG FOR891); The Gabon data set was provided by the Institut de Recherche en Ecologie Tropicale (IRET)/Centre National de la Recherche Scientifique et Technologique (CENAREST); Dutch inventory data collection was done with the help of Probos, Silve, Bureau van Nierop and Wim Daamen, financed by the Dutch Ministry of Economic Affairs. Data collection in Middle Eastern countries was supported by the Spanish Agency for International Development Cooperation [Agencia Española de Cooperación Internacional para el Desarrollo (AECID)] and Fundación Biodiversidad, in cooperation with the governments of Syria and Lebanon. We are grateful to the Polish State Forest Holding for the data collected in the project “Establishment of a forest information system covering the area of the Sudetes and the West Beskids with respect to the forest condition monitoring and assessment” financed by the General Directory of State Forest National Holding. We thank two reviewers who provided constructive and helpful comments to help us further improve this paper. The data used in this manuscript are summarized in the supplementary materials (tables S1 and S2). All data needed to replicate these results are available at https://figshare.com and www.gfbinitiative.org. New Zealand data (doi:10.7931/V13W29) are available from S.W. under a materials agreement with the National Vegetation Survey Databank managed by Landcare Research, New Zealand. Access to Poland data needs additional permission from Polish State Forest National Holding, as provided to T.Z.-N.
View Abstract

Navigate This Article