Forest production efficiency increases with growth temperature
Definitions of terms
GPP is defined here as ‘the sum of gross carbon fixation (carboxylation minus photorespiration) by autotrophic carbon-fixing tissues per unit area and time54. GPP is expressed as mass of organic carbon produced per unit area and time, over at least one year. NPP consists of all organic carbon that is fixed, but not respired over a given time period54:
$${mathrm{NPP}} = {mathrm{GPP}}-R_{mathrm{a}} = {Delta}B + L + F + H + O = {mathrm{BP}} + O$$
(3)
with all terms expressed in unit of mass of carbon per unit area and time. Ra is autotrophic respiration (composed of growth and maintenance respiration components); ΔB is the annual change in standing biomass carbon; litter production (roots, leaves and woody debris) is L; fruit production is F; the loss to herbivores is H, which was not accounted here because of the very limited number of observations available. BP is biomass production4. Symbol O represents occult, carbon flows, i.e. all other allocations of assimilated carbon, including changes in the nonstructural carbohydrate pool, root exudates, carbon subsidies to symbiotic fungi (mycorrhizae) or bacteria (e.g. nitrogen fixers), and BVOCs emissions (Supplementary Fig. 1). These ‘occult’ components are often ignored or unaccounted when estimating NPP, hence this bias is necessarily propagated into the Ra estimate when Ra is calculated as the difference between GPP and NPP55.
Estimation methods
We grouped the ‘methods’ into four categories:
biometric: direct tree stock measurements, or proxy data together with biomass expansion factors, allometric equations and the stock change as a BP component. If not otherwise stated, we assumed that the values included both above- and below-ground plant parts (n = 13 for GPP; n = 200 for NPP or BP).
micrometeorological: micrometeorological flux measurements using the eddy-covariance technique to measure CO2 flux and partitioning methods to estimate ecosystem respiration and GPP (n = 98 for GPP; n = 4 for NPP or BP).
model: model applications ranging from single mathematical equations (for canopy photosynthesis and whole-tree respiration) to more complex mechanistic process-based models to estimate GPP and Ra, with NPP as the net difference between them (n = 53 for GPP; n = 24 for NPP or BP).
scaling: upscaling of chamber-based measurements of assimilation and respiration (GPP and Ra) fluxes at the organ scale, or the entire stand (n = 73 for GPP; n = 9 for NPP or BP).
The difference between ‘scaling’ and ‘modelling’ lies in the data used. In the case of ‘scaling’ the data were derived from measurements at the site. ‘Model’ means that a dynamic process-based model was used, but with parameters calibrated and optimized at the site, based on either biometric or micrometeorological measurements.
Data selection
The data were obtained from more than 300 peer-reviewed articles (see also ref. 5), adding, merging and extending published works worldwide on CUE or BPE4,9,11,23,25,56,57. Data were extracted from the text, Tables or directly from Figures using the Unix software g3data (version 1.5.2, Jonas Frantz). In most studies, NPP, BP and GPP were estimated for the tree stand only. However, GPP estimated from CO2 flux by micrometeorological methods applies to the entire stand including ground vegetation. We therefore included only those micrometeorological studies where the forest stand was the dominant primary producer. The database contains 244 records (197 for BPE and 47 for CUE) from >100 forest sites (including planted, managed, recently burned, N-fertilized, irrigated and artificially CO2-fertilized forests; Supplementary Information, Supplementary Fig. 3 and online Materials; https://doi.org/10.5281/zenodo.3953478), representing 89 different tree species. Globally, 170 records out of the total data are from temperate sites, 51 from boreal, and 23 for tropical sites, corresponding to 79 deciduous broad-leaf (DBF), 14 evergreen broad-leaf (EBF), 132 evergreen needle-leaf (ENF) and 19 mixed-forests records (MX). The majority of the data (∼93%) cover the time-span from 1995 to 2015. We assume that when productivity data came from biometric measurements the reported NPP would have to be considered as BP because ‘occult’, nonstructural and secondary carbon compounds (e.g. BVOCs or exudates) are not included. In some cases, multiple datasets from the same site were included, covering different years or published by different authors. We considered only those values where either NPP (or BP) and GPP referred to the same year. From studies where data were available from more than 1 year, mean values across years were calculated. When the same reference for data was found in different papers or collected in different databases, where possible, we used data from the original source. When different authors described the same values for the same site, one single reference (and value) was used (in principle the oldest one). By using only commonly available environmental drivers to analyse the spatial variability in CUE and BPE, we were able to include almost all of the data that we found in the literature. We examined as potential predictors site-level effects of: average stand age (n = 204; range from 5 to ∼500 years), mean annual temperature (MAT; n = 230; range −6.5 to 27.1 °C) and total annual precipitation (TAP; n = 232; range from ∼125 to ∼3500 mm yr−1), method of determination (n = 237), geographic location (latitude and longitude; n = 241, 64°07′N to −42°52′S and 155°70′W to −173°28′E), elevation (n = 217; 5–2800 m, above sea level), leaf area index (LAI, n = 117; range from 0.4 to 13 m2 m−2), treatment (e.g.: ambient or artificially increased atmospheric CO2 concentration; n = 34), disturbance type (e.g.: fire n = 6; management n = 55), and the International Geosphere-Biosphere Programme (IGBP) vegetation classification and biomes (n = 244), as reported in the published articles (online Materials). The methods by which GPP, NPP, BP (and Ra) were determined were included as random effects in a number of possible mixed-effects linear regression models (Supplementary Table 4).
We excluded from statistical analysis all data where GPP and NPP were determined based on assumptions (e.g. data obtained using fixed fractions of NPP or Ra of GPP). In just one case GPP was estimated as the sum of upscaled Ra and NPP58; however, this study was excluded from the statistical analysis. NPP or Ra estimates obtained by process-based models (n = 23) were also not included in the statistical analysis. No information was available on prior natural disturbance events (biotic and abiotic, e.g. insect herbivore and pathogen outbreaks, and drought) that could in principle modify production efficiency, apart from fire. The occurrence of fire was reported by only a few studies59,60,61. These data were included in the database but fire, as an explanatory factor, was not considered due to the small number of samples in which it was reported (n = 6).
Data uncertainty
Uncertainties of GPP, NPP and BP data were all computed following the method based on expert judgment as described in Luyssaert et al.55. First, ‘gross’ uncertainty in GPP (gC m−2 yr−1) was calculated as 500 + 7.1 × (70−|lat|) gC m−2 yr−1 and gross uncertainties in NPP and BP (gC m−2 yr−1) were calculated as 350 + 2.9 × (70−|lat|). The absolute value of uncertainty thus decreases linearly with increasing latitude for GPP and for NPP and BP, because we assumed that the uncertainty is relative to the magnitude of the flux, which also decreases with increasing |lat|. Subsequently, as in Luyssaert et al.55, uncertainty was further reduced considering the methodology used to obtain each variable, by a method-specific factor (from 0 to 1, final uncertainty (δ) = gross uncertainty × method-specific factor). Luyssaert et al.55 reported for GPP-Micromet a method-specific factor of 0.3 (i.e. gross uncertainty is reduced by 70% for micrometeorological measurements); and for GPP-Model, 0.6. GPP-Scaling and GPP-Biometric were not explicitly considered in ref. 55 for GPP. We we used values of 0.8 and 0.3, respectively. For BP-Biometric and NPP-Micromet we used a reduction factor of 0.3; for NPP-Model, 0.6; and for NPP-Scaling (as obtained from chamber-based Ra measurements), 0.8. When GPP and/or NPP or BP methods were not known (n = 7), a factor of 1 (i.e. no reduction of uncertainty for methods used, hence maximum uncertainty) was used. The absolute uncertainties on CUE (δCUE) and BPE (δBPE) were considered as the weighted means62 by error propagation of each single variable (δNPP or δBP and δGPP) as follows:
$$delta {mathrm{CUE}} = sqrt {left( {frac{{delta {mathrm{NPP}}}}{{{mathrm{GPP}}}}} right)^2 + left( {delta {mathrm{GPP}}frac{{{mathrm{NPP}}}}{{{mathrm{GPP}}^2}}} right)^2}$$
(4)
and similarly for δBPE, by substituting NPP with BP and CUE with BPE.
Data and model selection
The CUE and BPE data were combined into a single variable, as sites for which both types of estimates existed did not show any significant differences between these entities (Supplementary Fig. 2). CUE values based on modelling were excluded (in our database we do not have BPE data from modelling). Tests showed that the CUE value was systematically higher when GPP was estimated with micrometeorological methods, compared to values based on biometric or scaling methods. Only data with complete information on CUE, MAT, age, TAP, and latitude were used. Altogether, 142 observations were selected.
In order to use the most complete information possible, a full additive model was constructed first (Eq. (1)). The method used for estimation of GPP (GPPmeth) was specified as a random effect on the intercept, as visual inspection suggested that CUE values were smaller where ‘scaling’ was used to estimate GPP compared to cases where ‘micromet’ was used to estimate GPP.
In Eq. (1) the variable ‘age’ represents the development status of the vegetation, i.e. either average age of the canopy forming trees or the period since the last major disturbance. The other three parameters represent different aspects of the climate. The absolute latitude, |lat|, was chosen as a proxy of radiation climate, i.e. day length and the seasonality of daily radiation. The term ηZGPPMeth represents the random effect on the intercept due to the different methods of estimating GPP.
These variables were not independent (Supplementary Table 1). If the different driver variables contain information that is not included in any of the other driver variables, multiple linear regression is nonetheless able to separate the individual effects. If, on the contrary, two variables exert essentially the same effect on the response variable (CUE) this can be seen in an ANOVA based model comparison. These considerations led us to the selection procedure in which we started with the full model (Eq. (1)) and compared it with all possible reduced models (Supplementary Table 2). The result of this analysis is the model with the smallest number of parameters that does not significantly differ from the full model.
We also examined, whether there were any significant interactions of predictor variables. There were not.
We used the R function lmer from the R-package lme463 to fit the mixed and ordinary multiple linear models to the data. We checked for potential problems of multicollinearity using the variance inflation factor (VIF)64. All predictors had VIF More
