Nitrogen and phosphorus constrain the CO 2 fertilization of global plant biomass

Elevated CO 2 (eCO 2 ) experiments provide critical information to quantify the effects of rising CO 2 on vegetation 1,2,3,4,5,6 . Many eCO 2 experiments suggest that nutrient limitations modulate the local magnitude of the eCO 2 effect on plant biomass 1,3,5 , but the global extent of these limitations has not been empirically quantified, complicating projections of the capacity of plants to take up CO 27,8 . Here, we present a data-driven global quantification of the eCO 2 effect on biomass based on 138 eCO 2 experiments. The strength of CO 2 fertilization is primarily driven by nitrogen (N) in ~65% of global vegetation and by phosphorus (P) in ~25% of global vegetation, with N- or P-limitation modulated by mycorrhizal association. Our approach suggests that CO 2 levels expected by 2100 can potentially enhance plant biomass by 12 ± 3% above current values, equivalent to 59 ± 13 PgC. The global-scale response to eCO 2 we derive from experiments is similar to past changes in greenness 9 and biomass 10 with rising CO 2 , suggesting that CO 2 will continue to stimulate plant biomass in the future despite the constraining effect of soil nutrients. Our research reconciles conflicting evidence on CO 2 fertilization across scales and provides an empirical estimate of the biomass sensitivity to eCO 2 that may help to constrain climate projections.


Introduction
Levels of eCO 2 affect the functioning and structure of terrestrial ecosystems and create a negative feedback that reduces the rate of global warming 8,9,11,12,13,14 . However, this feedback remains poorly quantified, introducing substantial uncertainty in climate change projections 7,8 . Experiments with eCO 2 simulate the response of plants to eCO 2 and thereby provide important empirical and mechanistic constraints for climate projections. Numerous eCO 2 experiments have been conducted over the last three decades and they collectively provide strong evidence for a fertilizing effect of eCO 2 on leaf-level photosynthesis 6 . At the ecosystem level, however, individual CO 2 experiments show contrasting results for the magnitude of the growth and biomass response to eCO 2 , ranging from strongly positive in some studies 2 to little or no response with N 1 , P 5 or water 3 limitations in other studies. Despite this conflicting evidence at the ecosystem scale, a global-scale carbon (C) sink in terrestrial ecosystems is robustly inferred 12 .
Here, we synthesize 1,432 observations from 138 eCO 2 studies in grassland, shrubland, cropland and forest systems (Supplementary Figs. 1 and 2 and Supplementary Table 1), encompassing free-air CO 2 enrichment (FACE) and chamber experiments. We train a random-forest meta-analysis model with this dataset and identify the underlying factors that explain variability within it. We use these relationships to estimate the global-scale change in biomass in response to an increase in atmospheric CO 2 from 375 ppm to 625 ppm, which is the increase in CO 2 expected by 2100 in an intermediate emission scenario.
We included 56 potential predictors of the CO 2 effect (Supplementary Table 2) belonging to four main categories: nutrients (N, P, mycorrhizal association; see ref. 4 ), climate (for example, precipitation and temperature), vegetation (age and type) and experimental methodology (for example, the increase in CO 2 concentration (∆CO 2 ) and the type of CO 2 fumigation technology). More details on the model selection are available in the Supplementary Discussion.
The random-forest meta-analysis indicated that the most important predictors of the CO 2 fertilization effect on biomass in our dataset were experiment type (FACE or chambers), soil C:N ratio (an indicator of N availability), soil P availability and mycorrhizal type, with different relationships for C:N and P between mycorrhizal types (y ≈ Mycorrhizal_type × N + Mycorrhizal_type × P + Fumigation_type, pseudo-R 2 = 0.94). A sensitivity test using a larger dataset of 205 studies confirmed the robustness of the relationships described by the statistical model (Supplementary Discussion). Among 56 potential predictors, mycorrhizal type was the primary modulator of above-ground biomass responses to eCO 2 (P < 0.001) ( Supplementary Fig. 3).
The eCO 2 effect in arbuscular mycorrhizal (AM) plants was best predicted by soil C:N ( Fig. 1a, P < 0.001), but not significantly by P ( Supplementary  Fig. 4a, P = 0.2830). The C:N ratio of soil organic matter is a proxy for plant N availability because it is associated with stoichiometric limitations of microbial processes in the soil 15 . Although the constraining role of N on CO 2 fertilization has been reported in many eCO 2 studies 1,3,6 , here we find that soil C:N is a powerful indicator to quantify the N-limitation on CO 2 fertilization across experiments.
In contrast, the eCO 2 effect in ectomycorrhizal (ECM) plants was best predicted by soil P (Fig. 1b, P < 0.001), but not significantly by soil C:N ( Supplementary Fig. 4b, P = 0.1141). The critical role of P on CO 2 fertilization across a large number of studies was unexpected, but consistent with an increasing body of research 5,16 .
Once the effects of mycorrhizal type, C:N, P and fumigation type were accounted for, other predictors such as climate, biome type (for example, temperate tree versus grass) or the age of the vegetation did not explain an important fraction of the variability in the effect ( Supplementary Fig. 5). Previous studies have variously attributed differences in the magnitude of the CO 2 effect to either average temperature (MAT) or precipitation (MAP), or to both 17 (see Supplementary Discussion). Using the model y ≈ MAT + MAP + Fumigation_type instead of our final model reduced explained variability (R 2 ) from 0.94 to 0.05. These results suggest that the CO 2 fertilization can only be reliably predicted when nutrient availability is considered.
We used the quantitative relationships derived from the meta-analysis to predict the global distribution of the eCO 2 effect based on maps for soil C:N, P and mycorrhizal type. Plant responses to eCO 2 were significantly higher in open top chamber and growth chamber experiments than in FACE ( Supplementary Fig. 5, P < 0.001) (see Supplementary Discussion), so we included Fumigation_type as a predictor in the scaling model to produce projections that are consistent with the response found in FACE experiments, as they allow CO 2 to be fumigated with as little disturbance as possible.
Our global projections from FACE experiments show a relative increase in biomass of 12 ± 3% ( Fig. 2a and Table 1) for the average 250 ppm ∆CO 2 across experiments. The magnitude of the global effect is less than the overall effect of ~20% found previously in meta-analyses 4,6 and the ~30% effect found in several FACE experiments 2,4 . This reduction arises in part because many CO 2 experiments were conducted in relatively fertile soils or under nutrient fertilization regimes. Thus, extrapolating nutrient relationships to areas with naturally poor soils results in a lower global effect. In absolute terms, we estimated a global increase in total biomass of 59 ± 13 PgC for a 250 ppm ∆CO 2 ( Fig. 2b and Table 1), scaled from satellite observations of current above-ground biomass 18 and region-specific total to above-ground biomass ratios from the literature (Supplementary Table 4). Global anthropogenic emissions are currently around 10 PgC annually 12 , hence the additional C-sequestration in biomass is equivalent to 5-6 years of intermediate CO 2 emissions.
Forests show the largest relative increases in biomass (Table 1 and Fig. 2a). Tropical forests are characterized by low P (Supplementary Fig. 6). However, their association with AM fungi, together with relatively high N ( Supplementary Fig. 6), support a widespread, though moderate, biomass enhancement. Our approach does not explicitly include symbiotic acquisition of atmospheric N (N 2 -fixation), which is relatively common in tropical forests 19 . Indeed, tropical N 2 -fixing species can show larger CO 2 effects than non N 2 -fixing species 20 , and thus the response in tropical forests in our model may be underestimated. Nevertheless, our dataset contains tropical N 2 -fixing species 21 , indirectly including this effect. Temperate grasslands, which are also dominated by AM plants, show the lowest relative biomass increment as a result of N-limitations. In temperate forests, some of the largest relative increases (~30%) occur in ECM-forests when P is high, but AM-forests show low relative biomass increases due to moderately high C:N ( Supplementary  Fig. 6).
The absolute eCO 2 effect is dominated by tropical forests (Table 1 and Fig. 2b), consistent with ground-based measurements showing increases in above-ground biomass in recent decades in intact tropical forests 22 , with CO 2 identified as the main driver 22,23 . To account for uncertainties, and to highlight the environmental conditions not well represented in eCO 2 experiments, we computed the standard error of the projections (Methods). Wet-tropical and boreal ecosystems show the largest uncertainties in absolute and relative terms, respectively (Fig. 2c,d), reflecting the limited number of studies in ecosystems with extreme values of climate and nutrient availability.
To assess the magnitude of the global eCO 2 effect we derive from FACE, we compared it with the increase in biomass attributed to rising CO 2 concentration (β) from 1980 to 2010 by the TRENDY ensemble of dynamic global vegetation models (DGVMs), standardized to 100 ppm ∆CO 2 . Our estimated rate of increase in total biomass is 25 ± 4 PgC 100 ppm −1 , a value within the range of DGVMs and slightly larger than the multimodel ensemble mean β (Fig. 3a). This similarity is remarkable given the independency of both approaches and reported large inconsistencies in DGVMs in partitioning total to above-ground biomass 24 .
For comparing the geographical distribution of our global eCO 2 effect, we used satellite-based observations of changes in leaf area (greening) 9 attributed to CO 2 rising in the period 1982-2009. Although changes in greenness and above-ground biomass are not necessarily correlated, we found an intriguingly strong correlation between the contemporary CO 2 -driven increase in greenness and our independently estimated biomass projections (Fig. 3b,c).
In summary, our results suggest that plant biomass responses to eCO 2 are driven primarily by interactions with N and P modulated by mycorrhizal status. N constrains the strength of CO 2 fertilization in most AM plants (Fig. 1a), which currently store ~65% of terrestrial vegetation C 25 , probably because the ability of AM fungi to supply plants with N is relatively small 26,27 . In contrast, we observed that P availability alters the biomass response to eCO 2 in ECM plants, which store ~25% of terrestrial vegetation C 25 . The sensitivity of ECM plants to P availability may be driven by the positive effect of eCO 2 on N uptake in ECM plants 27 , which, together with widespread N deposition, might reinforce the limiting role of P 28 in the ecosystem.
Although our analysis uses the most comprehensive dataset of eCO 2 observations currently available, it has several limitations. First, our data-driven approach, unlike DGVMs, is not intended to capture the complex interactions that drive long-term changes in the C cycle, such as warming, disturbance, changes in water availability or N deposition. Instead, it is aimed at the empirical quantification of net CO 2 effects, providing constraints on the attribution of modelled biomass responses to CO 2 and a better mechanistic understanding of the underlying drivers of the effect. Second, tropical and boreal ecosystems are under-represented in global eCO 2 experiments (Supplementary Fig. 1). We have accounted for this uncertainty in our estimates, which we also use to highlight the specific regions where eCO 2 experiments are urgently needed. Furthermore, it is critical that comprehensive soil data in eCO 2 experiments are reported, ideally in more long-term studies.
We observed a strong similarity between the global-level responses to eCO 2 found in FACE and past changes in biomass and greening attributed to CO 2 . The implications of this finding are threefold. First, this convergence supports our projections, indicating that empirical relationships with soil nutrients can be powerful for explaining large-scale patterns of eCO 2 responses, despite ecosystem-level uncertainties. Second, the effect attributed to rising CO 2 in past decades by DGVMs is similar in magnitude to our predicted effect of increasing CO 2 expected in the future (Fig. 3a), suggesting that the past CO 2 fertilization effect may continue at a similar magnitude for some time, despite nutrient limitations. Third, all else being equal, the same ecosystems that are currently responsible for most of the greening 9 and C uptake 11,14 are likely to remain important for future increases in biomass under eCO 2 (see Fig. 3b,c).
A key strength of our upscaling approach is that it synthesizes observational evidence at local scales and captures a global view of the eCO 2 effect on plant biomass and its drivers. DGVMs differ at the process level (including the current effects of CO 2 on biomass, see Fig. 3a), and consequently vary when projecting the future. Our data-based approach, along with new data from ongoing experiments, can be updated continuously and used to calibrate DGVMs, providing an empirical constraint for model simulations of the biomass sensitivity to CO 2 . This research accounts for the extent of nutrient limitations on the eCO 2 fertilization effect and shows that, despite local limitations, a global and positive effect, consistent with independent evidence of past CO 2 fertilization, can be inferred. This result challenges the strong and pervasive limitations on the projected eCO 2 fertilization suggested by some nutrient-enabled models 29 . For example, in the TRENDY ensemble of models in Fig. 3a, only OCN and CLM4CN take N limitations into account, and none of them to our knowledge include P limitations. While model simulations of the CO 2 effect on biomass by OCN closely match our data-driven results, CLM4CN underestimates the CO 2 fertilization effect by half and thus overestimates nutrient limitations. This may be related to the limited capacity of plant N uptake to mediate an excessively open N cycle in CLM4CN 30 .
Our results highlight the key role of terrestrial ecosystems, in particular forests, in mitigating the increase in atmospheric CO 2 resulting from anthropogenic emissions. Thus, if deforestation and land use changes continue decreasing the extent of forests, or if warming and other global changes diminish or reverse the land carbon sink, we will lose an important contribution towards limiting global warming.

Overview
The goal of this paper is to scale the effects of eCO 2 on biomass globally. This scaling requires a quantification of 'current' plant biomass and its distribution worldwide together with a model based on the environmental drivers (predictors) that statistically best explain the observations derived from eCO 2 studies. We collected data on above-ground biomass (Supplementary Figs. 1 and 2 and Supplementary Table 1) because (1) above-ground biomass is the metric most commonly reported in eCO 2 studies and (2) satellites can only detect above-ground biomass; thus, upscaling the effects of eCO 2 on above-ground biomass avoids some of the uncertainties related to modelled products of plant productivity or total (above-ground and below-ground) biomass.
From an initial pool of 56 potential predictors, we selected the most important predictors based on variable importance metrics from randomforest meta-analysis. We built a mixed-effects meta-regression model with the most important predictors of the effect, and applied this model with global maps to scale the effects of eCO 2 on above-ground biomass.
Finally, our results were evaluated in terms of distribution and magnitude. For the distribution of the effect, we compared the latitudinal distribution of our estimates with the latitudinal effects of CO 2 on changes in greenness (LAI) in the past three decades 9 . For the magnitude of the effect, we compared our sensitivity of biomass changes to eCO 2 with the sensitivity of biomass changes to the historical increase in atmospheric CO 2 (β) derived from the TRENDY ensemble of global vegetation models 10 .

Data collection
We collected 1,432 above-ground biomass observations from 205 studies that met our criteria (below), of which 138 had data for all predictors considered and were therefore included in our analysis. Repeated measurements over time within the same plots (that is, annual or seasonal measurements) were considered non-independent, and were thus aggregated so that only one synthetic measurement per study was included in the meta-analysis. Different species or treatments within the same site were considered independent, but we included 'site' as a random effect in the mixed-effects meta-analysis to account for this potential source of nonindependency (see Meta-analysis). We consulted the list of CO 2 experiments from INTERFACE (https://www.bio.purdue.edu/INTERFACE/experiments.php), the Global List of FACE Experiments from the Oak Ridge National Laboratory (http://facedata.ornl.gov/global_face.html), the ClimMani database on manipulation experiments (www.climmani.org) and the databases described by Dieleman et al. 31 ,Baig et al. 32 and Terrer et al. 4,27,33 . We used Google Scholar to locate the most recent publications for each of the previously listed databases.
We included as many observations as possible for our analysis. Criteria for exclusion from the main analysis were: (1) soil C:N and N content data for the specific soils in which the plants were grown were not reported-for example, studies that included a N fertilization treatment were only included when C:N was measured in situ, and not in unfertilized plots; (2) species did not form associations with either AM or ECM-only species in two studies were non-mycorrhizal, insufficient to identify the drivers of the eCO 2 response in this group; and (3) the duration of the experiment was less than 2 months.
We considered the inclusion of factorial CO 2 × warming or CO 2 × irrigation studies when specific soil data for those additional treatments were measured and reported. These treatments were treated as independent and were included in the dataset using the specific MAT and MAP for the warming and irrigation treatments, respectively. Approximately a quarter of the studies were irrigated, with irrigation more common in cropland studies. In those cases, and if the total amount of water used in irrigation was not indicated, we assigned the historical maximum value of MAP extracted from the coordinates of the site in the period 1900-2017 from ref. 34 . Although in some studies we found soil data for several soil depth profiles, soil data were most commonly reported for a depth of 0-10 cm. We thus collected soil data at 0-10 cm, and scaled CO 2 effects using global gridded datasets for this depth increment.
Data for MAP, MAT, soil C:N, soil N content, pH, available P and vegetative and experimental predictors were reported in the literature. Data for the rest of the predictors were not commonly reported, so we extracted these data from global gridded datasets (Supplementary Table 2).
We used the check-lists in refs. 35,36 , with additional classifications derived from the literature, to classify plant species as ECM, AM or non-mycorrhizal. Species that form associations with both ECM and AM fungi (for example, Populus spp.) were classified as ECM because these species can potentially benefit from increased N availability due to the presence of ECM fungi 4,27 , as hypothesized. Overall, CO 2 responses from species associated with AM and ECM were similar to strictly ECM species, and their exclusion did not alter the results of the meta-analysis, as found previously 4 .
Where possible, data were collected at the species level, and different species from the same site were considered independent when grown in monoculture with sufficient replication (that is, multiple plots of the same species and multiple individuals of the same species in the same plot).
Using these criteria, we found a total of 205 studies with data on aboveground biomass, with 138 of them including data for all the predictors considered, and thus included in the main analysis. Additionally, we ran a sensitivity test including data from our full dataset of 205 studies, estimating missing soil N and P data from proxies, in the following order of preference: (1) from studies that, due to proximity, used similar soils; (2) from gridded datasets (Supplementary Table 2) in the case of non-fertilized soils; and (3) using the mean values in the dataset for fertilized and non-fertilized studies within ecosystem types. For example, if a study comprised of temperate trees in a fertilized soil did not report soil data, and the characteristics of these soils could not be estimated from similar known soils, we assigned missing data as the average values in the dataset for temperate trees in fertilized soils. Supplementary Table 1, data included in the meta-analysis in Supplementary  Fig. 2 and location of the studies in Supplementary Fig. 1. An overview of the studies excluded from the main analysis is given in Supplementary Table 3, and included in a sensitivity test.

Model selection and relative importance
We used random-forest model selection in the context of meta-analysis to identify the most important predictors of the CO 2 effect in the dataset. This method has the advantage over maximum likelihood model-selection approaches that can handle many potential predictors and their interactions, and considers nonlinear relationships.
Some of the 56 potential predictors included in the analysis were extracted from global datasets using the coordinates of the experiments (Supplementary Table 2), and thus included missing values. Because random-forest and meta-analysis require complete data, and no methods for multiple imputation are currently available, we applied single imputation using the missForest 37 algorithm. Like any random forests-based technique, the main advantage of this method is that it does not make any distributional assumptions, which means it easily handles (multivariate) non-normal data and complex interactions and nonlinear relations amongst the data.
Some of the potential predictors provided redundant and potentially correlated information (that is, multiple methods to measure soil P and multiple climate predictors) (see Supplementary Table 2). We used principal component analysis (PCA) for dimensional reduction, extracting components from map-based, potentially redundant predictors.
We included all field-based predictors, together with PCA map-based predictors, in a bootstrapped random-forest meta-analysis recursive preselection with the metaforest 38 R package. We trained a random-forest meta-analysis with preselected predictors and calculated variable importance with metaforest 38 (Supplementary Fig. 3). Based on partial dependence plots ( Supplementary Fig. 5), we used reciprocal transformations for nonlinear predictors showing ceiling/floor effects. We included the ten most important predictors in a mixed-effects metaregression model with the metafor 39 R package, including reciprocal transformations for nonlinear predictors and potential interactions. Finally, we pruned the model once, keeping only significant predictors.
As a sensitivity test, we ran an alternative model-selection procedure using maximum likelihood estimation. For this purpose, we used the rma.mv() function from the metafor R package 39 and the glmulti() function from the glmulti R package 40 to automate fitting of all possible models containing the seven most important predictors and their interactions. Model selection was based on Akaike Information Criterion corrected for small samples as criterion, using a genetic algorithm for faster fitting of all potential models. The relative importance value for a particular predictor was equal to the sum of the Akaike weights (probability that a model is the most plausible model) for the models in which the predictor appears. A cut-off of 0.8 was set to differentiate between important and redundant predictors, so that predictors with relative importance near or less than 0.8 are considered unimportant.

Meta-analysis
We used the response ratio (mean response in elevated-to-ambient CO 2 plots) to measure effect sizes 41 . We calculated the natural logarithm of the response ratio (logR) and its variance for each experimental unit to obtain a single response metric in a weighted, mixed-effects model using the rma.mv function in the R package 'metafor' 39 . We included 'site' as a random effect (because several sites contributed more than one effect size and assuming different species or treatments within one site are not fully independent), and weighting effect size measurements from individual studies by the inverse of the variance 42 . Some 5% of studies did not report standard deviations, which were thus imputed using Rubin and Schenker's 43 resampling approach from studies with similar means and performed using the R package metagear 44 .
Measurements across different time-points (that is, over several years or harvests) were considered non-independent, and we computed a combined effect across multiple outcomes (for example, time-points) so that only one effect size was analysed per study. The combined variance that accounts for the correlation among the different time-point measurements was calculated following the method described in Borenstein et al. 45 , using a conservative approach by assuming non-independency of multiple outcomes (r = 1) and performed using the MAd package in R 46 .
We considered nonlinear mixed-effects meta-regression models, which were fitted using reciprocal transformations (1/variable).

Quantification of uncertainties
Extrapolating the empirical relationships that drive biomass responses to eCO 2 (for example, y ≈ C:N) in the dataset to the globe has an error associated with the mixed-effects meta-regressions. For the case of soil C:N, for example, this error is large for high C:N values, as the representativeness of soils with high C:N values in the dataset is lower, increasing uncertainty in the regression. For the projections of the eCO 2 effect, we limited the maps of C:N and P to be constrained by the minimum and maximum values in the dataset of eCO 2 studies, thus assuming saturating responses to avoid extremely high or low (negative) effects that are not representative of the observed effects. For the projection of uncertainties in Fig. 2, however, we aimed at representing not only the uncertainty associated with the representativeness of the most important predictors (C:N and P), but also the uncertainty associated with the lower sampling effort in areas with extreme climate (for example, very dry and warm-deserts-or cold and dry-boreal -or wet and warm-tropical). We therefore ran an alternative model that included temperature and precipitation in addition to C:N and P. We extrapolated the standard error of this alternative model using unconstrained maps of temperature, precipitation, C:N and P to account for the higher level of uncertainty in areas with climate and soil values that are not well represented by eCO 2 experiments. Thus, uncertainties in our projections represent the unconstrained standard error of the mixed-effects meta-regression, with larger values under soil and climate conditions that are not adequately studied due to low sample size.
Global estimates of N and P availability N can be limiting for plants (1) if there is little total N content or (2) because N is bound in organic matter with a high C:N ratio. In the latter case, soil microbes that degrade the organic matter become N-limited, resulting in low amounts of free N available for plant uptake. Therefore, soil N content and C:N ratio were included as potential predictors of the CO 2 effect. Other potential predictors, such as nitrate and ammonium contents and N mineralization, were not generally available and were therefore not included in the analysis. Because soil C:N ratio was an important predictor of the CO 2 -driven increase in biomass in our dataset ( Fig. 1a and Supplementary Fig. 3), we used a global dataset of soil C:N ratio from ISRIC-WISE on a 30 × 30 arcsec grid 47 to upscale this effect. The range of C:N values covered by eCO 2 experiments is representative of the range of C:N values represented in the C:N map 47 .
Arid regions typically have very low soil C:N ratios as a result of a small organic C pool and also low N content 48,49 . Therefore, soil C:N is not a good indicator of N availability in arid soils, and the model would overestimate the CO 2 effect in these areas because it would assume relatively high N availability. To avoid the overestimation of the CO 2 effect in arid areas with low C:N, yet low N availability, we followed the approach of Wang et al. 50 , who found a threshold of 0.32 in aridity index (ratio of precipitation to mean temperature) below which plant N uptake is limited by water availability, and characterized by low soil C:N despite extremely low soil N availability. We converted areas with aridity index <0.32 to null values in the map of soil C:N, thereby treating these areas as missing data for analyses including soil C:N. We used the aridity data from the CGIAR-CSI Global-Aridity Database 51 . In our dataset of CO 2 experiments, the Nevada Desert FACE fell within this category, with low soil C:N, but low total N 52 , and no CO 2 effect on biomass 53 , supporting this assumption. Running the model strictly in areas with aridity index >0.32 resulted in 0.4 PgC less than by running the model globally. This small difference was the result of the extremely low above-ground biomass in arid regions ( Supplementary Fig. 7), rendering small absolute increases in biomass when incorporated in the analysis. Nevertheless, these areas were not included in the final analysis because it is not likely they could increase their biomass under elevated CO 2 due to extremely low N availability. In areas outside this maximum aridity threshold limiting nitrogen uptake, we studied the impact of climatic and water availability predictors in explaining the magnitude of the CO 2 effect.
The amount of P in the soil estimated by the Bray method was one of the important predictors of the biomass responses to eCO 2 . We constrained the map of P amount by the minimum and maximum values of P in the dataset of eCO 2 studies, 2-64 ppm, assuming these values are representative of the conditions at <2 and >64 ppm.

Climate data
For the model selection analysis (Fig. 2) we used MAT and MAP data for the individual studies reported in the papers. As an alternative climatic predictors to MAT and MAP to account for the effect of temperature and water availability, we tested additional predictors not commonly reported in the papers, calculated using temperatures and precipitation values from CRU or extracted from other gridded datasets (Supplementary Table 2).

Current above-ground biomass
As global estimates of current above-ground biomass carbon we used passive microwave-based global above-ground biomass carbon from Liu et al. 18 (v.1.0) at 0.25° resolution and available online for the period 1993-2012 (http://www.wenfo.org/wald/global-biomass/).

Land cover types
Calculations of changes in biomass in response to CO 2 across biomes were performed through zonal statistics with the land cover maps from ESA (http://maps.elie.ucl.ac.be/CCI/viewer/download.php) at 300 m resolution (Table 1) and MODIS IGBP (http://glcf.umd.edu/data/lc/) at 5´ resolution (Supplementary Table 4). Both maps were aggregated by dominant classes. The indication of climatic region (that is, temperate, boreal, tropical) within forest land cover types was based on the classification by Pan et al. 54 .

Changes in LAI
In order to evaluate the geographical patterns of our predictions, we compared the latitudinal distribution of the effects of elevated CO 2 on aboveground biomass with changes in LAI attributed to CO 2 in the period 1982-2009 (ref. 9 ). We used LAI data from three different satellite records and averaged them, as described in ref. 9 . The attribution of the relative and absolute effects of CO 2 on LAI was estimated through vegetation models, as described in Zhu et al. 9 .
For the calculation of the effects of elevated CO 2 on biomass, regions where water availability limits N uptake (aridity index < 0.32) were excluded from the analysis (see Global estimates of N and P availability). Thereby, for the comparison of biomass and LAI changes, these arid regions were excluded from both maps.

Global vegetation models
In order to evaluate the magnitude of the sensitivity of plant biomass to eCO 2 derived from our analysis, we analysed biomass β for the historical increase in atmospheric CO 2 derived from the DGVMs considered in the TRENDY intercomparison project (http://dgvm.ceh.ac.uk/node/9). We used TRENDY-v1, which includes nine DGVMs with common input forcing data, varying CO 2 only from 1980 to 2010 (S1) and calculated biomass β as the change in biomass relative to the change in atmospheric CO 2 . For more details on the TRENDY model simulations see Sitch et al. 10 .

Calculation of total biomass carbon
The TRENDY models considered here output total biomass (above ground + below ground), whereas our results refer to above-ground biomass only. In order to compare the magnitude of the eCO 2 effect derived from models and our approach, we have estimated the potential effect of eCO 2 on total biomass using region-specific ratios of total biomass and above-ground biomass reported in the literature (Supplementary Table 4).

Data availability
The biomass data from CO 2 experiments summarized in Supplementary  Fig. 2 supporting the findings of this study are available in published papers, and soil and climate data required to upscale CO 2 effects are available in published datasets (Supplementary Table 2). Raw data can be obtained from the corresponding author on reasonable request.

Code availability
The R code used in the analysis presented in this paper is available online and can be accessed at https://github.com/cesarterrer/CO2_Upscaling.