The enhanced vegetation productivity driven by increased concentrations of carbon dioxide (CO 2 ) [i.e., the CO 2 fertilization effect (CFE)] sustains an important negative feedback on climate warming, but the temporal dynamics of CFE remain unclear. Using multiple long-term satellite-and ground-based datasets, we showed that global CFE has declined across most terrestrial regions of the globe from 1982 to 2015, correlating well with changing nutrient concentrations and availability of soil water. Current carbon cycle models also demonstrate a declining CFE trend, albeit one substantially weaker than that from the global observations. This declining trend in the forcing of terrestrial carbon sinks by increasing amounts of atmospheric CO 2 implies a weakening negative feedback on the climatic system and increased societal dependence on future strategies to mitigate climate warming.

T errestrial ecosystems have accounted for more than half of the global carbon sink during the last six decades and have thus substantially mitigated climate warming (1). Global process-based models attribute part of the increasing land carbon sink (2) to the increase in vegetation productivity driven by the fertilization effect of increasing atmospheric CO 2 concentration (3), i.e., the CO 2 fertilization effect (CFE), a process that acts as a negative feedback in the climate system (4). First introduced by Keeling (5), the b factor generally is used to characterize the plant response to increasing CO 2 concentration. To compare the results between global analysis and experimental measurements, we used an approximation form of the original b (b ¼ @GPP @Ca ) (6), which is defined as the relative increase in gross primary production (GPP) in response to a 100-ppm increase in atmospheric CO 2 concentration (C a ).

Increases of GPP originate from the direct acceleration of photosynthesis in response to the increased supply of CO 2 (7), but they are also modified by a suite of indirect responses, including water saving because of the reduced stomatal conductance under increasing atmospheric CO 2 concentrations (8,9) and nutrient limitation (10). Enhanced GPP and net primary production (NPP) are commonly observed in field experiments, such as free-air CO 2 enrichment (FACE) experiments (7) and open-top chamber experiments (11), where ecosystems are exposed to elevated CO 2 in the range of two times ambient values, albeit with different response magnitudes across sites (7).

The fingerprint of CFE on the carbon and water cycles on global scales (3,8) is, however, more elusive, given the covariation of atmospheric CO 2 with other environmental drivers of vegetation productivity. Elucidating this fingerprint of CFE is a scientific problem involving detection and attribution methods, requiring either statistical methods and longterm data such that effects of non-CFE drivers can be removed empirically (12,13), or process-based models with which CFE can be isolated by deliberate factorial simulations (14).

Process-based models of the terrestrial carbon cycle have indicated that CFE accounts for ~70% of the increasing global trend in foliar area, i.e., global greening (15), and up to 60% of the current terrestrial carbon sink (3). These models project that CFE will induce increased land carbon storage by the end of this century, despite opposing effects from climate change (3,16). A cascade of uncertainties in projecting future land carbon storage arises from the responses of photosynthesis, NPP, and ecosystem carbon turnover times to increasing CO 2 , especially in the presence of increasing limitations from nutrient and water availability. Therefore, accurately quantifying the temporal dynamics of CFE on GPP is essential for reducing the uncertainties of future land carbon storage and climate projections derived on the basis of Earth system models (4).

A study using the long-term CO 2 concentration measurements at the Barrow station (Alaska, 71°N) found that the sensitivity of the seasonal amplitude of atmospheric CO 2 to the increase in CO 2 concentration has decreased, which suggests a declining response of GPP to CO 2 at northern high latitudes (17). If such declining trends prevail across the globe, the terrestrial vegetation carbon sink response to increasing CO 2 will decrease, with important consequences for the global carbon budget and for the effort required from mitigation policies to meet future climate targets. For this reason, a comprehensive assessment of the temporal dynamics of CFE on global GPP is timely as a first step to understand the impacts of CFE on the trends of the land carbon sinks. We used three long-term satellite datasets collected during 1982 to 2015: a recently developed vegetation index (NIR V ) (18) from the Advanced Very High Resolution Radiometer (AVHRR NIR V ), the fusion of NIR V from AVHRR and the Moderate Resolution Imaging Spectroradiometer (AVHRR+MODIS NIR V ), and the fusion of NIR V and sun-induced chlorophyll fluorescence (SIF) (19) (NIR V +SIF), which are proxies for GPP (hereafter, satellite GPP proxies) (supplementary text S1 and figs. S1 to S5) (20). We also corroborated our findings using long-term GPP time series from eddy covariance (EC) flux towers, from a light use efficiency (LUE) model (21), and from an ensemble of terrestrial carbon cycle models (TRENDY v6, trends and drivers of the regional scale sources and sinks of carbon dioxide) (14).

Linear and nonlinear multiple regression approaches were used to estimate b from satellite GPP proxies at each pixel across the globe (12,13,20). Pixels with large land-cover changes were excluded from the analyses (supplementary text S2 and fig. S6). The global median b during 1982 to 2015 was 16.1 ± 11.5% 100 ppm -1 (fig. S7), consistent with the FACE experiments (7), which on average suggested a 15.5% increase in light-saturated photosynthesis per 100 ppm CO 2 . Results showed that the estimated b values for different vegetation physiology and biome types were also aligned with experimental results (7), e.g., the b values of C 4 plants were much smaller than those of C 3 plants (supplementary text S3 and fig. S8). We next calculated the time series of b with 15-year moving windows during 1982 to 2015 and found that b significantly decreased at a rate of -0.92 ± 0.12% 100 ppm -1 year -1 (p < 0.01) (Fig. 1A). This decrease was also evident after excluding the crop areas (fig. S9 A B C D Fig. 2. Global declining trend of b. Spatial patterns of trends in b derived from three datasets using 15-year moving windows, reported as the percentage per 100 ppm per year. (A) Means from three satellite GPP proxies. (B) The revised EC-LUE GPP. (C) Multimodel mean GPP determined from the TRENDY v6 ensemble. All datasets are from 1982 to 2015. The Mann-Kendall test was used to estimate the trend of b pixel by pixel, and the regions with black dots indicate significant trends (p < 0.05). The pixel size is 1°. (D) Mean declining trend of b from different datasets. The error bars represent SEs.

significantly declined from 22.1 ± 13.3% 100 ppm -1 during 1982 to 1996 to 12.9 ± 11.3% 100 ppm -1 during 2001 to 2015 (Fig. 1B).

We also found a robust decrease in b when we used another satellite NIRv time series (fig. S11), various land cover change masks (fig. S12), different window lengths (supplementary text S4 and fig. S13), various combinations of explanatory variables (supplementary text S5 and figs. S14 and S15), and different definitions of growing seasons (fig. S16) and when we included the consideration of seasonal precipitation (supplementary text S6 and fig. S17). After considering the uncertainties from original satellite data and from the methods, the global b estimated from satellite GPP proxies during 2001 to 2015 was still significantly lower than that during 1982 to 1996 (supplementary text S7 and figs. S18 and S19). To verify whether the trend of b was an artefact of our regression method, we also estimated the temporal dynamics of b using the optimal fingerprint attribution method (supplementary text S8) and again found that the declining trend of b was significant (fig. S20). Altogether, these results based on satellite observations of GPP proxies suggest a significant declining trend in the response of GPP to increasing atmospheric CO 2 .

To further assess the decreases in b inferred from satellite GPP proxies, we used other independent datasets: (i) the GPP time series derived from the revised EC-LUE model (21) (which also accounted for the direct CO 2 effects on LUE); (ii) satellite-based leaf area index (LAI) time series; (iii) multiyear GPP estimations from 22 EC flux sites (table S4); and (iv) a global long-term GPP dataset from the Carbon Cycle Data Assimilation System (CCDAS) (supplementary text S9 to S12). We also found a significant declining trend of b on the basis of the revised EC-LUE GPP, with a rate of -0.82 ± 0.18% 100 ppm -1 year -1 (fig. S21), similar to the estimates from satellite GPP proxies. However, the result from satellitebased LAI showed a smaller declining rate of b (-0.59 ± 0.11% 100 ppm -1 year -1 ) (fig. S22). This divergence suggested that the decrease in b was likely due to two factors: the direct effect on foliar physiology and the indirect effect on LAI. The former refers to the CO 2induced stimulation of carboxylation per unit leaf area, whereas the latter depicts the increased carbon sequestration used for leaf area expansion. The smaller declining b trend from satellite LAI suggested that these two effects were both relevant for the full description of the global decreases in the CO 2 effect on GPP. Ground GPP estimates from EC flux sites also confirmed our findings, by indicating an average declining rate of b at -0.70% 100 ppm -1 year -1 , which was comparable to the results from satellite GPP proxies around these sites (fig. S23). On the basis of the CCDAS GPP, a data-model fusion product constrained by atmospheric CO 2 observations and which is independent from satellite data, we also found a declining b trend at a rate of -0.62% 100 ppm -1 year -1 (fig. S24). The overall consistency between various remote-sensing data, ground measurements, and a carbon assimilation system suggests that the global decline of b is robust and coherent across multiple, independent observations.

The geographic distribution of the temporal trends of b from satellite GPP proxies (Fig. 2A) revealed that b decreased in ~86% of global lands (fig. S25, A to C). Areas with declining b spanned over most of the globe. By contrast, increasing b was found in some limited regions of Southeast Asia, east Australia, and North America (Fig. 2, A and B). These increasing b trends are likely driven by an intensification of management in croplands (i.e., irrigation and fertilization) or related to in-creasing atmospheric nutrients deposition in recent years (fig. S33). We also observed similar spatial patterns for b trends calculated using both 10-year and 17-year moving windows (fig. S26). Overall, the decreases of b in tropical areas were smaller, whereas cold regions had a slightly stronger declining trend (figs. S10A and S25, A to C). The latter may be because the declining signals of b in northern high latitudes include both the indirect effects on LAI and the direct effect on foliar physiology, whereas those in warmer climates include only the physiological effect, given that LAI in some tropical humid regions may already be close to saturation. Results of the satellite LAI confirmed this interpretation by suggesting a larger reduction of b in colder northern regions than in warmer climates (fig. S25, E and F). Similar declines in b were obtained from the revised EC-LUE GPP, which covered ~74% of the global terrestrial area (Fig. 2B). In general, all datasets indicated that b has been declining in most global lands and across various vegetation types and climate zones.

We also investigated whether an ensemble of state-of-the-art carbon cycle models ( 14) was able to reproduce these observed global declines of b on GPP. We used results from 12 models that contributed to the TRENDY v6 ensemble to calculate b (20) (supplementary text S13). For each model, we extracted the differences in simulated GPP between S1 (timevarying CO 2 only) and S0 (constant CO 2 ) scenarios, which represented the CO 2 effect on vegetation productivity. These models predicted a negative trend in b, both for the multimodel mean and for the individual one, but the declining rate was clearly lower than that derived from satellite GPP proxies and with no evident spatial variations (Fig. 2C and fig. S27). Grouping the estimates of CFE trends, we found that the global b declined at a rate of -0.92 ± 0.12% 100 ppm -1 year -1 for satellite GPP proxies and only -0.12 ± 0.01% 100 ppm -1 3 of 6

year -1 for TRENDY GPP (Fig. 2D). b trends of individual models were variable, from -0.06 to -0.21% 100 ppm -1 year -1 (fig. S28), but all of them were lower than that from satellite GPP proxies. Given that we defined b as the percent increase in GPP per 100-ppm increase in CO 2 and that the S1 simulations solely considered the CO 2 effect, the smaller decreases in b from TRENDY GPP were likely caused by the saturating physiological response of GPP to CO 2 (22), without adequately capturing the concurrent emergence of other limiting factors driven by the changing environmental conditions. We also used GPP from the TRENDY S2 and S3 scenarios to estimate the trends in b by using the regression approaches (20).

The temporal trends of b remained unaltered, with a value of approximately -0.12 ± 0.12% 100 ppm -1 year -1 for both TRENDY S2 and S3 simulations (fig. S29). These results highlight that the ongoing strong decrease in b inferred from satellite datasets is probably underestimated by TRENDY models.

Two possible non-mutually exclusive hypotheses were proposed to account for the declining b and to explain why TRENDY models failed to adequately replicate the magnitude of this decline: (i) the increasing constraints on vegetation productivity from emerging nutrient limitations and imbalances that were not adequately represented in models and (ii) current models underestimated the sensitivity of terrestrial GPP to changes in water availability.

The first hypothesis relates to the possible effect of the growing limitation of key nutrients, including nitrogen (N) and phosphorus (P). Using 410 groups of ground-based foliar N and P observations (supplementary text S14 and fig. S30), we found a general decrease in foliar N and P concentrations, with mean values of -0.24 ± 0.06% year -1 and -0.55 ± 0.06% year -1 , respectively (fig. S31, A and B). Our findings are consistent with a recent study suggesting a general global pattern of decreasing foliar N concentration (23) and with many examples of local to regional decreases in foliar N and P concentrations (24,25). Enhanced GPP from the increasing atmospheric CO 2 concentrations, leading to larger NPP and higher nutrient demands by plants, may partly explain the observed declines of foliar N and P concentrations. Concurrently, these decreases in key foliar nutrients may impose limitations on GPP and thereby limit b. To test this hypothesis, we applied a linear spatial mixedeffects model (table S5) to investigate the relationship between the trend of b and foliar N or P concentrations (supplementary text S14 and fig. S32). After accounting for the trends in b explained statistically by mean annual air temperature and mean annual precipitation, we found clear positive correlations between the residual declining trends of b with both foliar N and P concentrations across European forests (Fig. 3). These results suggested that vegetation with lower foliar nutrient concentrations generally showed larger declines in b, therefore supporting our hypothesis of a role for nutrient limitation in the temporal dynamics of b.

This phenomenon was also verified from global trends in atmospheric N and P depositions, which suggested that areas with decreasing atmospheric nutrient supplies probably had larger decreases in b (such as Europe and Siberia) and vice versa (such as East Asia) (fig. S33). This may be because vegetation in areas with declining nutrient supplies from atmospheric depositions tend to have larger N and P limitations on GPP. Moreover, the increases of N and P depositions in some regions of East Asia may explain the increasing trend of b in these areas (Figs. 2A and fig. S33). The ongoing decreases in foliar nutrients might constrain the plant photosynthetic capacity and result in the decline of b, which might not have been adequately represented in current models (sup-plementary text S14). Regarding this aspect, the models that included C-N cycle interactions to emulate nutrient constraints exhibited a larger declining rate of b (fig. S34), partly confirming this interpretation. The role of N limitation on b has been widely suggested by experimental evidence (10,26,27), model analyses (28), and synthesis reviews (29,30), all of which have been consistent with our analyses. Our finding regarding P limitation was also consistent with FACE experiments, which demonstrated that insufficient P availability generally had negative impacts on b (31). Furthermore, as foliar P:N ratios positively correlated with plant net photosynthesis and growth (32), the decreasing foliar P:N ratios in European forests (fig. S31C) suggested a worsening nutrient imbalance that may partly account for the observed decline in b.

The second hypothesis to explain the divergence of the b trend between the TRENDY models and the satellite-derived estimates stated that these models underestimated the sensitivity of GPP to water availability, because the coupling between water and carbon in models is underestimated (33). To test this hypothesis, we analyzed the sensitivity of satellite GPP A B C RESEARCH | RESEARCH ARTICLE Erratum 5 February 2020. See Erratum. on June 29, 2021 http://science.sciencemag.org/ Downloaded from

proxies and TRENDY GPP to water availability using a moving window of 15 years (supplementary text S15). We used the terrestrial water storage (TWS) data to represent the availability of water to plants and selected arid, temperate, and dry tropical climate zones as the research areas (fig. S35), as the GPP in these zones was found to be highly sensitive to TWS (fig. S36). We observed that the sensitivity of GPP to TWS was relatively constant across these three zones for the multimodel mean of the TRENDY GPP, whereas the TWS sensitivities derived from satellite GPP proxies exhibited significantly increasing trends for arid (2.93 ± 0.38% year -1 ), temperate (2.12 ± 0.39% year -1 ), and dry tropical (2.30 ± 0.57% year -1 ) ecosystems (Fig. 4). In temperate and dry tropical areas, the majority of models exhibited divergent results compared to satellite GPP proxies, and almost all of them largely underestimated the GPP sensitivities to TWS in arid areas (fig. S37). The decreases in b were larger in the regions where the increases of TWS sensitivities were higher (fig. S38). These strongly divergent sensitivities between satellite GPP proxies and TRENDY GPP were also confirmed when using various TWS products (20) (supplementary text S15 and fig. S39) and when using a shorter moving window (fig. S40).

This finding implies that GPP and consequently land carbon uptake are more sensitive to the variations in water availability than assumed by the TRENDY models, as previously suggested by several recent studies (34,35). Moreover, a recent study showed that the CFE of grasslands could be reduced under drier conditions (36). The relationship between drought stress and CFE is complex and ecosystem specific and may be affected by the total annual rainfall as well as rainfall seasonality. It was shown, for instance, in 19 grasslands experiments that the annual CFE was negated when spring precipitation became too low (37). Field experiments also support our results by suggesting that the CFE on vegetation productivity is at least partly limited by water availability (38). The significant increases in the GPP sensitivities to TWS may thus partly explain the decreases of b in arid, temperate, and dry tropical climate zones.

From a theoretical perspective, the global declines of CFE may result from several factors. First, given that the CO 2 -induced photosynthesis stimulation at the leaf level is scaled up to the canopy level through LAI (39,40), the declining b on the basis of satellite LAI (fig. S22) could partly explain the global decreases of the GPP response to CO 2 . Consistent with our findings, a recent study using the FACE experiments on mature forests found a relatively low CO 2 effect on GPP (41), possibly because the LAI of these forests did not change much (39). Second, the CO 2 effect on photo-synthesis at the leaf level involves both the stimulation of carboxylation and increases in water use efficiency (WUE). The former may possibly be regulated by the foliar key nutrients (N and P), and the latter is related to the water availability. According to the progressive N limitation theory, the soil N availability for plant growth may be expected to diminish over time (42), possibly leading to the observed global decreases in b. A recent study using carbon isotope measurements revealed diminishing CO 2 -induced WUE gains across global forests (43), which supported our findings of the effects of water supply limitation on CFE and could partly explain the global declines in b. Moreover, complex interactive effects between nutrients and water supply may also have impacts on CFE. For example, low soil water supply could possibly strengthen the nutrient constraints on CFE through limiting nutrient decomposition and diffusion in soils (38). Our further analysis using model simulations showed that the GPP trends have clear reductions in arid areas, once the interactions between N limitation and climate constraints are considered, supporting this hypothesis (fig. S41B and supplementary text S16). Nevertheless, an excess of N fertilizer may result in a reduction in the soil water and therefore possibly lead to drought stress (44). Regarding this aspect, our analysis also highlights the need for future efforts to better understand the complex interactions between nutrients, water, and CO 2 effects on vegetation from a climate change perspective. Additional mechanisms (e.g., plant acclimation or changes in plant species over time) may also explain the observed global decrease in CFE (supplementary text S17).

Our analyses showed a significant and spatially extensive decline in b, which implies a substantial reduction of the positive effects of increasing atmospheric CO 2 on terrestrial carbon uptake. A recent study suggested that the CO 2 effect on the carbon cycle in tropical regions (3) might be counteracted by impacts from climate-driven changes (45), in agreement with our findings. Although still under debate, the possible increasing trend of the airborne fraction of anthropogenic CO 2 may imply a saturation of the CO 2 sinks from land and oceans (46)(47)(48)(49), which may be partly caused by the global decline in CFE. Current carbon cycle models also exhibit such global decreases in b but fail to adequately detect the sharp declining trend that we identified from satellite data. This divergence between observations and process-based models possibly originates from the models' limitations in adequately representing the emerging decline in key foliar nutrient concentrations and the increasing constraints of water limitations on vegetation productivity. Ultimately, these results indicate that terrestrial photosynthesis may not increase as much as models project, possibly reducing the potential of land-based climate mitigation, further accelerating global warming and exacerbating the efforts required for meeting climate targets. Our findings also highlight the need for better characterizations of the biogeochemical and hydrological effects on vegetation in current carbon cycle models to produce more robust projections of the terrestrial carbon budget for the next decades.