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 time⁵⁴. 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 period⁵⁴:

with all terms expressed in unit of mass of carbon per unit area and time. R_a 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 production⁴. 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 R_a estimate when R_a is calculated as the difference between GPP and NPP⁵⁵.

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 CO₂ 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 R_a, 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 R_a) 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. ⁵), adding, merging and extending published works worldwide on CUE or BPE^{4,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 CO₂ 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 CO₂-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⁻¹), 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 m² m⁻²), treatment (e.g.: ambient or artificially increased atmospheric CO₂ 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 R_a) 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 R_a of GPP). In just one case GPP was estimated as the sum of upscaled R_a and NPP⁵⁸; however, this study was excluded from the statistical analysis. NPP or R_a 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 studies^59,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.⁵⁵. First, ‘gross’ uncertainty in GPP (gC m⁻² yr⁻¹) was calculated as 500 + 7.1 × (70−|lat|) gC m⁻² yr⁻¹ and gross uncertainties in NPP and BP (gC m⁻² yr⁻¹) 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.⁵⁵, 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.⁵⁵ 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. ⁵⁵ 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 R_a 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 means⁶² by error propagation of each single variable (δNPP or δBP and δGPP) as follows:

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 ηZ_GPPMeth 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 lme4⁶³ 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)⁶⁴. All predictors had VIF < 5 (between 1.1 and 3.8).The model residuals were also tested for normality (using the Anderson-Darling test of non-normality, in the R-package nortest⁶⁵). For models that did not take a random intercept regarding ‘GPPmeth’ into account (16–30 in Supplementary Table 2) the Anderson-Darling test found significant deviation from normality of the model residuals, hence these models were excluded from the analysis. The remaining models were compared with one another using the function ANOVA of the R-package lmerTest⁶⁶. This resulted in a 15 × 15 matrix of model comparisons in which the full model turned out to be significantly different from all other models.

The same analysis was also performed with a log-transformed version of Eq. (1):

where 7.5 °C was added to MAT in order to make its minimum 1 °C. Note that the linear model from the log-transformed variables differs from the untransformed linear model. The coefficients, here noted with a prime, can be interpreted on the basis of the back-transformed model. Contrary to the untransformed linear model where effects are additive, the back-transformed model is a multiplicative effect model, with the slope parameters as exponents for each variable and the intercept ((beta _0^prime) as power of e). As with the untransformed model, negative slope parameter values lower CUE, positive increase it with increasing driver variable values.

The results from this analysis were, as with the original additive model (Eq. (1)), (i) the full model could not be reduced any further and (ii) the directions of the effects were the same as with the additive model, i.e. the predicted CUE increased with increasing MAT, TAP and |lat| but decreased with increasing age.

The AIC and BIC values were lower for the log-transformed model compared to the untransformed model, with AIC values of −169.7 and −157.2 and BIC values of −149.0 and −136.5 for the log-transformed and untransformed models, respectively. The coefficients and model performance parameters of the untransformed and the log-transformed models are shown in Table 1 and Supplementary Table 4. The adjusted squared correlation coefficients were similar: 0.306 for the untransformed and 0.321 for the log-transformed model. Despite considerable uncertainty of the CUE values, it was possible to derive significant, systematic, linear relationships between the four driver variables and CUE or ln (CUE). Both model variants showed the same direction and similar magnitudes of the effects. It can be concluded that CUE (or ln (CUE)) from a global dataset of a large variety of forests is significantly positively affected by MAT, TAP and |lat|, and significantly negatively affected by age. Even excluding from the analysis the five tropical forest data with |lat| < 20 degrees did not alter significantly the empirical relationship (Supplementary Table 5).

Because the parameters of the untransformed, additive model are much easier to interpret, we use the additive model in the main text and use the log-transformed model only as a confirmation of trends found in the additive model.

Outputs from TRENDY v.7

We used the simulations from eight Dynamic Global Vegetation Models (DGVMs) performed in the framework of the TRENDY v.7 project^2,67 (http://dgvm.ceh.ac.uk/node/9; data downloaded 27 November 2019). Models that did not provide NPP and GPP at plant functional type level were excluded because of the need to analyse CUE in forests without significant contributions from shrubs, grassland or crops. The selection comprises the following models: ISAM, JULES, LPJ-GUESS, CABLE-POP, ORCHIDEE, ORCHIDE-CNP, JSBACH and SDGVM (for references on models see refs. ^2,67 and Supplementary Table 6). All the models represent the surface fluxes of CO₂, water and the dynamics of carbon pools in response to changes in climate, atmospheric CO₂ concentration, and land-use change across a global grid. However, processes underlying the exchanges of water and carbon are based on different formulations in different models.

In the TRENDY protocol all DGVMs were forced with common historical climate fields and atmospheric CO₂ concentrations over the period from 1700 to 2017. Climate fields were taken from the CRU-JRA55 dataset², whereas the time series of atmospheric CO₂ concentrations were derived from the combination of ice core records and atmospheric observations. Land-use change was taken into account in the simulations (S3). However, similar simulations without land-use change (S2) were also tested, showing no differences. CUE was estimated as NPP/GPP (where NPP is commonly obtained in models by subtracting R_a from GPP) for the forest plant functional types simulated to be present in each grid cell. The model outputs refer to the mean from 1995 to 2015 for comparability with the records used when showing global land analysis (Fig. 5 and Supplementary Fig. 4). At site level, the same dates as the observations were chosen from the model outputs.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com