CO2 fertilization of terrestrial photosynthesis inferred from site to global scales

Edited by Christopher Field, Stanford Woods Institute for the Environment, Stanford University, Stanford, CA; received August 24, 2021; accepted January 12, 2022
March 1, 2022
119 (10) e2115627119


The magnitude of the CO2 fertilization effect on terrestrial photosynthesis is uncertain because it is not directly observed and is subject to confounding effects of climatic variability. We apply three well-established eco-evolutionary optimality theories of gas exchange and photosynthesis, constraining the main processes of CO2 fertilization using measurable variables. Using this framework, we provide robust observationally inferred evidence that a strong CO2 fertilization effect is detectable in globally distributed eddy covariance networks. Applying our method to upscale photosynthesis globally, we find that the magnitude of the CO2 fertilization effect is comparable to its in situ counterpart but highlight the potential for substantial underestimation of this effect in tropical forests for many reflectance-based satellite photosynthesis products.


Global photosynthesis is increasing with elevated atmospheric CO2 concentrations, a response known as the CO2 fertilization effect (CFE), but the key processes of CFE are not constrained and therefore remain uncertain. Here, we quantify CFE by combining observations from a globally distributed network of eddy covariance measurements with an analytical framework based on three well-established photosynthetic optimization theories. We report a strong enhancement of photosynthesis across the observational network (9.1 gC m−2 year−2) and show that the CFE is responsible for 44% of the gross primary production (GPP) enhancement since the 2000s, with additional contributions primarily from warming (28%). Soil moisture and specific humidity are the two largest contributors to GPP interannual variation through their influences on plant hydraulics. Applying our framework to satellite observations and meteorological reanalysis data, we diagnose a global CO2-induced GPP trend of 4.4 gC m−2 year−2, which is at least one-third stronger than the median trends of 13 dynamic global vegetation models and eight satellite-derived GPP products, mainly because of their differences in the magnitude of CFE in evergreen broadleaf forests. These results highlight the critical role that CFE has played in the global carbon cycle in recent decades.
Vegetation photosynthesis is responsible for the largest flux of carbon from the atmosphere into the biosphere (1). Both theory and experimental observations show that the resulting gross primary production (GPP) increases with rising atmospheric CO2, a process known as the CO2 fertilization effect (CFE) (28). The CFE plays a strong role in offsetting anthropogenic emissions by directly increasing the terrestrial carbon assimilation rate. It affects both the global carbon (6) and water budgets (9), which leads to changes in temperature and precipitation, and indirectly affects climate change through vegetation–climate feedbacks (10) by increasing plant growth (11). Understanding the CFE is thus critical to understanding the evolution of Earth’s climate.
Despite the importance of global photosynthesis, however, there are few long-term records of observationally inferred GPP in natural environments. This leads to large uncertainty in the magnitude of photosynthetic change over time and the CFE (4, 1214). Eddy covariance (EC), a measurement of gas exchange deployed at hundreds of ecosystems worldwide, provides estimates of photosynthesis that could contain information on the CFE (14). Estimating CFE from such observations is challenging, however, because of the large confounding effects from climatic variability of natural environments and the short duration of measurements at many sites. Studies that have attempted to isolate the CFE from carbon fluxes have been limited either to simple statistical approaches that preclude causative attribution (15, 16) or to process-based approaches that, despite capturing the sign of responses, report widely varying sensitivities (17, 18). This uncertainty has led to debate regarding the magnitude of the CFE, as evidenced by the large spread between satellite- and Earth system model (ESM)–based CFE estimates (1923). Indeed, because of the lack of observational evidence at the ecosystem scale, many GPP products based on satellite-derived metrics do not explicitly account for CFE (20, 24, 25), while current representations of CFE also vary markedly among ESMs (2628).
Here, we leverage globally distributed EC observations and reanalysis data to detect and attribute CFE using an eco-evolutionary optimality (EEO) framework. This framework reconciles three well-established optimization theories in combination with Fick’s law and the Farquhar-von Caemmerer-Berry (FvCB) photosynthesis model (29), constraining photosynthesis and key intermediate variables, namely stomatal conductance (g) (30, 31), intercellular leaf CO2 concentration (ci) (32, 33), and photosynthetic capacity (3438) (see Materials and Methods for general descriptions and SI Appendix for derivations). The framework can analytically diagnose GPP and its sensitivity to seven measurable variables, namely atmospheric CO2 concentration (ca), leaf area index (LAI), air temperature (Ta), volumetric soil water content (SWC), specific humidity (qa), incident shortwave radiation (SWin), and surface pressure (P), without the need to prescribe biome-specific photosynthetic properties (35). Our results show an overall increase in photosynthesis at measurement sites worldwide and attribute a large proportion of the increase to the effect of rising CO2.

Strong CO2-Induced Increase in Site-Scale Photosynthesis.

Three metrics are used to quantify and assess the response of GPP to CO2.
βCO2 (gC m−2 year−1 ppmv−1) is the partial differential sensitivity of GPP to CO2 for given environmental conditions (i.e., GPPca in Eq. 1). It is directly diagnosed by the EEO framework, which decouples the confounding effects on GPP from other environmental conditions.
βappln is a dimensionless, apparent logarithmic GPP response ratio to CO2 following Ref. 8 for comparison with other studies (Eq. 3). Changes in GPP include both CO2 and non-CO2 effects.
βdirln is a dimensionless, direct logarithmic GPP response ratio to CO2, which means changes in GPP include CO2 effects only. Confounding effects are neglected or analytically decoupled through ideal experimental manipulations such as free-air CO2 enrichment (FACE) experiments or partial derivative approaches. It shares the same equation as βappln, and is only used for comparison with other studies.
We find a strong increasing trend of GPP measured across 632 site-years of observations in the EC network of 9.1 gC m−2 year−2 between 2001 and 2014 [interquartile range (IQR) ∈ [8.5, 11.5], P < 0.01, βappln = 1.24]. This EC-inferred trend is captured by our EEO framework (8.3 gC m−2 year−2, IQR ∈ [7.0, 9.4], P < 0.01, βappln = 1.12) (Figs. 1A and 2A). Both βappln values are relatively high compared with expectations from measured direct responses to CO2 (i.e., βdirln) but are similar to several published indirect apparent GPP proxies (8): carbonyl sulfide records (βappln = 0.95) (39), ice-core measurements of atmospheric O2 isotopes (βappln = 1.3) (40), and satellite monitoring of water-use efficiency (βappln = 1.1) (41). Since βappln is derived from natural environments and influenced by CO2 and climate variability, it differs from the direct βdirln discussed later. The detected overall trends are robust to the uncertainty caused by filtering of low-quality data, and the uneven distribution of sites and site-years (SI Appendix, Fig. S2). Variations in the EC-inferred GPP due to partitioning methods and friction velocity filtering methods have a minimal effect on the EEO-inferred GPP (shaded area in Fig. 1A), where the former is used to calibrate the canopy upscaling of the latter. In general, the EEO-inferred GPP reproduces the interannual variability (IAV) of the EC-inferred GPP (Pearson’s r = 0.91) (Fig. 1A; SI Appendix, Fig. S1). The high level of agreement from annual (Fig. 1B) to monthly (SI Appendix, Fig. S3) scales supports the use of the EEO framework for attribution analysis.
Fig. 1.
The EEO-inferred GPP reproduces the trend and IAV of site-scale EC-inferred GPP. (A) Trends and interannual variations. GPP anomaly is defined as the deviation from climatological mean during the study period, further averaged across all 68 sites (see anomalies for individual sites in SI Appendix, Fig. S1). Blue numbers indicate the number of sites with available GPP. The shaded area represents the range of variation due to partitioning and friction velocity filtering methods in EC-inferred GPP (gray), and their corresponding variability propagated to the EEO-inferred GPP through canopy upscaling calibration (light red). Solid lines are the ensemble mean of EC-inferred GPP based on the above methods (black) and the resulting EEO-inferred GPP (red). (B) Scatter plot of EEO-inferred and EC-inferred GPP. Dashed lines indicate best-fit lines.
To quantify the GPP trends contributed by different factors, we use the EEO framework to perform a univariate sensitivity analysis (SI Appendix, Fig. S4). The aggregated trend from individual factors (10.3 gC m−2 year−2, IQR ∈ [7.7, 11.0]) is similar to the EEO-inferred and EC-inferred trends (Fig. 2A). We find a strong direct contribution from ca at the site level, which accounts for 44% (4.5 gC m−2 year−2, P < 0.001) of the aggregated GPP trend (Fig. 2A). We also diagnosed another overall CO2-induced GPP trend using the partial differential approach (4.9 gC m−2 year−2, the first term on the right-hand side of Eq. 1), resulting in a small difference (0.4 gC m−2 year−2) compared with the univariate analysis. We then convert the estimate from the univariate analysis, which corresponds to a direct response ratio (βdirln) of 0.61 for GPP to CO2. This βdirln is lower than those light-saturated leaf-level analyses in FACE experiments (βdirln = 0.79) (4, 8) and those using deuterium isotopomers (βdirln = 1.0) (8, 42) but is comparable to a theoretical βdirln of 0.6 in Ref. 8 that considers canopy radiative transfer, suggesting a coexisting of light-saturated and light-limited photosynthesis over our study sites. CFE is stronger under light-saturated conditions but is present regardless of light conditions (27, 28). An increase in ca eventually elevates leaf ci through optimization (SI Appendix, Eq. S16), which in the FvCB scheme leads to higher light use efficiency and carbon assimilation rates (7, 32).
Fig. 2.
GPP trend and IAV are dominated by different factors at site-scale. (A) Probability density of GPP trends by 50% resampling of site-years for each site for 200,000 times (see illustration in SI Appendix, Fig. S3D). Attribution of GPP trends is computed by univariate sensitivity analysis. The other variables are kept as the climatological mean of their study period. Blue * and × indicate statistically significant (P < 0.05) and insignificant trends in the estimated GPP, respectively. The hollow circles, thick black vertical lines, and gray horizontal lines are the median, IQR, and mean of trends. “EC-inferred”: GPP trends estimated from the EC network. “EEO-inferred”: GPP trends directly estimated by the EEO framework. “Aggregated”: sum of GPP trends for each factor from the sensitivity analysis. (B) GPP IAV is proxied by the product of partial derivatives of GPP to control factors and IAVs of control factors. Partial derivatives are diagnosed using climatological means of the control factors and IAVs are defined as 1 SD of the data with long-term trends removed. (C) GPP IAV induced by temperature IAV through different pathways. K is the Michaelis-Menten coefficient, Γ* is the CO2 compensation point, and η* is the ratio of water viscosity under ambient air temperature to a reference temperature at 298.15 K. In B and C, each violin consists of values from 68 sites.
In addition to ca, Ta plays the second most important role, accounting for 28% (2.9 gC m−2 year−2, P < 0.05) of the aggregated trend (Fig. 2A), mainly because of the effect of warming on photosynthetic capacity, i.e., the maximum rate of Rubisco activity, Vcmax, and the maximum rate of electron transport, Jmax (Fig. 2C) (38). Although increasing temperature negatively affects GPP through, for example, vapor pressure deficit (VPD) and some biochemical reaction pathways (Fig. 2C), these negative effects are weaker than the positive effects of temperature (Fig. 2 A and B). SWC and specific humidity, as proxies for plant water supply and atmospheric water demand, together contribute 14% of the aggregated trend (0.91 and 0.50 gC m−2 year−2, P > 0.05) but are statistically insignificant. These two factors (i.e., SWC and qa), however, are the largest contributors to GPP IAV, and their contributions (median values are 48 and 33 gC m−2 year−1, respectively; Fig. 2B) are an order of magnitude larger than those of the GPP trends (Fig. 2A). These results are consistent with previous data and ESM-driven studies on plant responses to water (4347) and are further supported by observations that GPP anomalies are regulated by stomatal conductance and leaf ci, as plants attempt to maximize their efficiency of carbon gain per unit water loss (30, 32, 48, 49).
Other factors, such as long-term changes in LAI and radiation, could also lead to changes in GPP across the EC network. However, we find a small and statistically insignificant trend in GPP due to increased LAI (0.88 gC m−2 year−2; Fig. 2A), which is consistent with a previous study that solely considered the response of GPP to changes in LAI (50). This small LAI-induced trend in GPP is comparable to the LAI trend, both of which are on the order of ‰ year−1 (51), although our framework diagnoses a large underlying sensitivity of GPP to LAI (mean GPPLAI = 119, intersite SD = 178, unit: gC m−2 year−1). In addition, light saturation in canopy radiative transfer (52) may lead to a low LAI-induced GPP trend as leaf area increases do not effectively translate to increases in the fraction of absorbed photosynthetically active radiation at high LAI (3, 20). For the changes in incoming shortwave radiation, their contribution to the aggregated GPP is also small (0.43 gC m−2 year−2; Fig. 2A) since there is no large trend in radiation over the study sites at the annual scale. Therefore, we conclude that neither changes in LAI nor radiation significantly contribute to the observed upward trend in GPP. The effect of changes in surface pressure is also negligible (Fig. 2 A and B).

Effects of CO2 Fertilization at the Global Scale.

To examine global CFE patterns, we apply our framework to Moderate Resolution Imaging Spectroradiometer (MODIS) LAI and European Centre for Medium-Range Weather Forecasts Reanalysis Version 5 – Land data (hereafter, ERA5-land). Here, our EEO framework is calibrated by the ensemble mean of all eight satellite-derived GPP products. We then assess both the GPP trend directly estimated by the EEO framework and the underlying differential sensitivities, which themselves can be used to reconstruct a composite trend estimate (Eq. 1). The two approaches show almost identical apparent GPP trends of about 4.7% decade−1 (P < 0.001) relative to their mean GPP throughout 2001 and 2016 (A0 and A1 in Fig. 3A). The resulting CO2-induced GPP trend is about 4.1% decade−1 (4.4 gC m−2 year−2, P < 0.001) (A2 in Fig. 3A; SI Appendix, Fig. S5). We note that the fractional contribution of ca to the apparent GPP trend at the global-scale analysis may be uncertain, and we do not report it because of the uncertainties in apparent GPP trend caused by the uncertain trends in climate forcings. However, estimates of the CO2-induced GPP trends (A2 in Fig. 3, calculated as the product of βCO2 and Δca) are independent of the uncertainties in the climate forcing trends because of the nature of our partial differential approach: 1) The spatial variation of βCO2 is rational because the climate forcings well capture their spatial variations at the global scale (5355), and 2) Δca is derived from the National Oceanic and Atmospheric Administration (NOAA)’s observations with high confidence. As a result, after excluding the spatial heterogeneity introduced by LAI and canopy structure, we show a good agreement when comparing the leaf-level CO2-induced GPP trend at the global-scale analysis with their EC site counterparts (r = 0.88; SI Appendix, Fig. S6A).
Fig. 3.
Trends in EEO-inferred, DGVM, and satellite-derived GPP at the global and ecosystem scales relative to their respective mean GPP from 2001 to 2016. (A) All biomes. (B) Individual biomes, including EBFs, other forests (OF), short woody vegetation (SW), grasslands (GRA), croplands (CRO), and C4 vegetation (C4). “A0”: estimated directly from our EEO framework. “A1”: diagnostic trends composited from partial derivatives using the EEO framework. “A2”: diagnostic CO2-induced trends through the partial derivative approach using the EEO framework. “B1”: 13 DGVM models, all forcing time-varying. “B2”: 13 DGVM models, CO2 only. “C1”: eight satellite-derived GPP products. For B1, B2, and C1, the bars represent their median. The gray bars represent GPP trends considering all effects. The red bars represent GPP trends caused by CO2 only. BEPS, Boreal Ecosystem Productivity Simulator; BESS, Breathing Earth System Simulator.
We compare our CO2-induced GPP trends (A2 in Fig. 3) with simulations from 13 dynamic global vegetation models (DGVMs) and eight satellite-derived GPP products during their overlap period (2001 to 2016). Our estimates are at the upper end of their trend distribution, at least one-third stronger than their median (Fig. 3A). The lower GPP trends in the median of DGVMs (B1) and satellite-derived products (C1) is likely due to their smaller CFE in evergreen broadleaf forests (EBFs) (Fig. 3B). Our EEO framework suggests similar relative GPP trends caused by CO2 among all C3 species (including EBFs), ranging from 4.3 to 5.1% decade−1 (A2 in Fig. 3B). We further replace the ensemble mean of all satellite-derived GPP with each individual product to calibrate the EEO framework, finding a small spread in the diagnosed CO2-induced relative GPP trends (SI Appendix, Table S2), despite the fact that these individual satellite GPP products used for calibration vary markedly (global GPP in 2001 with a minimum of MODIS Collection 5.5 [MOD-C55] at 104 Pg C year−1 and a maximum of Pmodel-s0 at 131 Pg C year−1). In EBFs, it is noteworthy that the satellite GPP products that only consider the indirect CFE on satellite reflectance, but not the direct CFE on gas exchange [i.e., FluxCom, MOD-C55, MOD-C6, and vegetation photosynthesis model (VPM)], do not report an increasing GPP trend (Fig. 3B). Their GPP trends are not a proxy of CFE because they model GPP solely based on climate forcings and leaf greenness or its surrogates. The climate forcings in reanalysis data may show negative impacts on GPP because of excessive temperature and drought stress in the tropics (see the EBF column in SI Appendix, Table S3), and the weak trends in leaf greenness cannot fully extrapolate the direct CFE effect because of the saturation of reflectance in EBFs (5662). While the choice of climate forcing could introduce some uncertainties to our CFE estimates, we find similar diagnostic results with another forcing dataset [i.e., Climatic Research Unit and Japanese 55-year Reanalysis (CRU-JRA55) in SI Appendix, Fig. S7]. Furthermore, we calculate βappln and βdirln for the global-scale analysis (SI Appendix, Table S3), which results in the same conclusions as those using the relative GPP trend metric.

The Partial Differential Sensitivity of GPP to CO2.

To further characterize the underlying CFE, we show the spatial distribution of βCO2, which disentangles the confounding effects from other factors. The magnitude of βCO2 is largest in hot and humid regions at the annual scale (Fig. 4A). Since the CO2 trend is assumed to be constant worldwide, the spatial pattern of βCO2 represents that of CO2-induced GPP trends (A2 in SI Appendix, Fig. S5). At the ecosystem level, βCO2 is highest in EBFs at about 5.8 gC m−2 year−1 per ppmv CO2, followed by croplands and other forests from the temperate to boreal regions at about 2.2 gC m−2 year−1 per ppmv CO2 (Fig. 4B). EBFs have the longest growing season, and their high temperature and radiation determine the high photosynthetic capacity (Vcmax and Jmax) and thus βCO2 (4, 28, 35). On the other hand, croplands cultivated in mild regions also have high photosynthetic capacity (63) and additional resource inputs from human management, including multiple cropping and irrigation (51), leading to a relatively high annual βCO2. Compared with C3 plants, our results suggest a positive but much weaker βCO2 in C4 plants (3, 4, 13), on the order of 0.6 gC m−2 year−1 per ppmv CO2, because the C4 CO2 pump inhibits the benefit of CFE (49).
Fig. 4.
Average sensitivity of annual GPP to CO2 over 2001 to 2018. (A) Spatially distributed sensitivity of GPP to CO2. (B) Ecosystem-level sensitivity of GPP to CO2. Black error bars indicate 1 SD spatial variation. Red error bars indicate 1 SD temporal variation. The EEO framework is calibrated by the ensemble mean of eight satellite-derived GPP products.
According to plant photosynthesis optimization theories and gauged by the GPP magnitude estimated from an EC network or satellites, our EEO framework analytically constrains the sensitivity of GPP to CO2 (i.e., βCO2) from the gas exchange perspective. Our results suggest that the analytical solution of βCO2 is insensitive to the GPP uncertainty used for calibration of the framework (Fig. 1A; SI Appendix, Table S2) and largely unaffected by interannual fluctuations of climatic conditions in recent years (Fig. 4B). This highlights the advantage of the EEO framework in that the partial derivatives decouple the confounding effects of climate on CFE. A recent study using regression methods showed a 40% decline in CFE from 1982 to 2015 (16). However, their key findings are subject to the large temporal uncertainty in the satellite data and limited by the causal ambiguity of their regression (6466). In addition, our EEO framework assumes that plants coordinate nutrient investment to optimize photosynthetic capacity (Vcmax and Jmax) with environmental constraints (3437), rather than prescribing specific values for plant functional types. Photosynthetic capacity not only controls the slope of the GPP response to CO2 but also strongly influences it by determining whether photosynthesis is light saturated or light limited (27, 28). We explore the influence of coordination at different timescales, comparing the results with the EC-inferred GPP (Fig. 1A; SI Appendix, Fig. S8), and adopt a decadal coordination to the environments in the peak LAI month for nutrient investment (67, 68) while adjusting the enzyme kinetics monthly with temperature (38). These assumptions result in a fixed reference Jmax/Vcmax ratio for each site (or grid) over the study period while still allowing both to vary spatially to account for different growing climates. The decadal coordination allows for the presence of light-saturated and light-limited photosynthesis because the coordination is set at a much longer timescale than the stomatal and leaf ci optimization (35). Under these conditions, the EEO framework diagnoses a large proportion of marginally light-saturated photosynthesis, which responds more strongly to changes in ca than does light-limited photosynthesis (27, 28). Specifically, 63 to 69% of C3 canopy photosynthesis originates from light-saturated conditions for all seasons using the ERA5-land forcing. Together, these features explain the high βCO2 predicted by our framework and also lead to a good agreement with the overall trend in EC-inferred GPP and supported by multiple independent studies (3942).


We detect a large increase in GPP from a globally distributed EC network. Our EEO framework successfully captures the trends and IAVs in the EC-inferred GPP. The site-scale analysis suggests that CO2 is the dominant contributor to the GPP trends, while SWC and specific humidity control the IAVs. Using this framework to scale the GPP globally, the estimated absolute CO2-caused GPP trend is comparable to its EC-inferred counterpart and translates this CO2 fertilization effect to a global increase in photosynthesis of 4.1% decade−1 since the 2000s (or 5.1 PgC decade−2 for global sum, 4.4 gC m−2 year−2 for global average). However, our EEO framework diagnoses a high CFE compared with the median values of the DGVMs and the satellite-derived GPP products, with the largest disagreement in tropical forests. These tropical forests are responsible for about one-third of global GPP, and thus, an accurate estimation of CFE in tropical forests is critical to global carbon cycle modeling. Finally, we urge expansion of the currently sparse EC observing network in tropical ecosystems to improve understanding of their physiological responses to CO2 and climate change.

Materials and Methods

EEO Framework for GPP and CFE.

We derive a parsimonious EEO framework that constrains plant photosynthesis and water loss at a monthly scale. In addition to attributing changes in GPP through factorial simulations, the framework allows for partial differential sensitivity of GPP to atmospheric CO2 and various measurable environmental conditions. This is particularly important as neither GPP nor CFE can be measured directly in the natural environment, and the apparent changes in GPP are due to the confounding effect of changes in CO2 and other environmental conditions. We use Fick’s law and FvCB photosynthesis model as the governing equations for leaf-level photosynthesis. Photosynthesis is then constrained by three well-established optimization theories, which are widely supported by theoretical and empirical evidence (3, 4, 30, 32, 3437, 44, 48, 6973). Reconciling these optimization theories aims to constrain key variables in photosynthesis (33), such as the marginal water-use efficiency (mWUE) (74), stomatal conductance (g), leaf intercellular CO2 concentrations (ci), and photosynthetic capacity. Details of the EEO framework are derived in the SI Appendix, Text S1 to S7. We briefly summarize the concept of this framework below.
According to Fick’s law of mass transfer, neglecting the leaf boundary layer and mesophyll resistances, leaf-level photosynthesis is the product of stomatal conductance (g) and the gradient between atmospheric (ca) and leaf intercellular CO2 concentrations (ci). For a first intuitive guess, an increase in ca can lead to increased photosynthesis. However, plants actively optimize the tradeoff between carbon gain and water loss. Optimization partially offsets the direct photosynthetic benefits of elevated ca by indirectly reducing g and increasing ci to increase the carbon gain per unit water loss eventually.
Theory 1 states that vegetation maximizes the sum of carbon gain and water loss through optimizing g in response to changing environment (30).
Theory 2 states that vegetation minimizes the summed cost of transpiration and carboxylation per unit carbon assimilation through optimizing ci (32).
Specifically, theory 1 constrains g, and theory 2 constrains ci. Combining theories 1 and 2 provides a direct constraint on the mWUE. The optimization of g and ci should be sufficient and necessary because not only can the variation in g affect ci by controlling the inflow rate of ambient CO2 (30), but also the variation in ci can regulate stomatal aperture by adjusting guard cell membrane potential (4, 75, 76). To further consider the water stress on photosynthesis, we use a logistic function to approximate the change of mWUE with SWC (SI Appendix, Fig. S9) (33, 49, 72, 77). This directly allows the framework to incorporate the impact of SWC into the photosynthetic optimization processes (SI Appendix, Fig. S10) by calibrating a site-specific constant (i.e., ζo; SI Appendix, Text S4).
In terms of the biogeochemical response, we assume that plants optimize nutrient allocation in response to the elevated CO2 and their growing environment. We use the photosynthetic coordination theory (theory 3) to constrain how plants optimize their nutrient investment, which allows for the estimation of photosynthetic capacity (i.e., Vcmax and Jmax) without prescribing biome types (3437).
Theory 3 states that at a longer timescale than theories 1 and 2, plants adjust reference photosynthetic capacity to a nearly equal limitation of photosynthesis under average daytime conditions by ribulose-1,5-bisphosphate carboxylase (Rubisco) activity and ribulose-1,5-bisphosphate (RuBP) regeneration (interchangeable with light saturated and light limited, respectively) (3437).
In other words, the coordination theory states that leaf nutrient content is regulated to reflect photosynthetic capacity, which is opposite to the hypothesis that photosynthetic capacity is limited by leaf nutrient content (78). Although the theory is well supported by multiple studies (3437), determining the causation between the photosynthetic capacity and leaf nutrient content, or indeed a third alternative of finding common factors that limit the two, warrants further investigation (78). Nevertheless, the coordination theory captures the correlation between photosynthetic capacity and leaf nutrient content and implies that nutrient limitation should be implicit in LAI (37). We tested different timescales of the photosynthetic capacity coordination with respect to the monthly optimization of g and ci and compared our results with the EC-derived GPP (SI Appendix, Fig. S8). We use the timescale that best captures the trend and IAV of the EC-derived GPP (Fig. 1A) as the basis of our analysis. That is, the reference rates of Vcmax and Jmax acclimate to the average environment of the peak LAI month during the study period (67, 68). For the other months, Vcmax and Jmax are a function of temperature and their reference values (38). Therefore, our framework allows the photosynthetic capacity to be optimized to the changing atmospheric CO2 and meteorological conditions on a decadal timescale, and the photosynthesis is roughly colimited by Rubisco and RuBP. Other timescales explored include reference Vcmax and Jmax values optimized according to 1) year-to-year variation of peak month CO2 + average of peak month meteorological conditions during the study period, 2) year-to-year variation of peak month CO2 and meteorological conditions, 3) month-to-month variation of CO2 + average of peak month meteorological conditions during the study period, and 4) month-to-month variation of CO2 and meteorological conditions. These represent a range of acclimation timescales, from short-term monthly and annual acclimation to longer-term decadal acclimation.
In order to preserve the parsimonious and analytical nature of the EEO framework, leaf-level photosynthesis is upscaled by LAI using a big-leaf approach. Changes in LAI implicitly account for canopy-scale effects of growing season extensions, recovery from disturbance, and succession. The upscaling factor for light extinction, denoted as G(μ)μ¯, represents the monthly average canopy shape (G) and solar zenith angle (μ) of a particular location (SI Appendix, Text S5). This factor is obtained by calibration with the EC-inferred GPP for site-scale analysis and satellite-derived GPP for global-scale analysis. G(μ)μ¯ for each month is then kept constant without year-to-year variation. We stress that the statistical significance of annual trends and interannual variations in EEO-inferred GPP are not caused by calibration but by the seven forcing variables: ambient CO2 concentration (ca), LAI, air temperature (Ta), volumetric SWC, specific humidity (qa), incident shortwave radiation (SWin), and surface pressure (P). We finally evaluate the EEO framework qualitatively, which meets six of the seven criteria proposed by Ref. 72 (SI Appendix, Text S6).
A key benefit of the analytical framework, in addition to allowing a direct estimate of GPP, is that it can be used to assess the partial derivatives of GPP to different variables, which can further be used to reconstruct a composite estimate of GPP changes for attribution analysis. The following equation can express the change in GPP and sensitivities of GPP to multiple factors:
where ΔGPP is the GPP change; GPPcaΔca is the direct CFE; GPPca is the partial derivative of GPP to CO2, which reflects the absolute response of GPP to unit change in CO2; and so on. We treat these variables as independent forcings to different physiological processes within the EEO framework, but these physiological processes eventually interact to influence the gas exchange. For those larger-scale land–atmosphere interactions, the covariances between the forcing variables are implicit in their measurements (or reanalysis data). We denote GPPca as βCO2 for simplicity. The symbol Δ can either be the interannual variation or long-term trend. To decompose the sensitivity of GPP to Ta into different pathways (i.e., Fig. 2C), we additionally use the chain rule of calculus, as follows.
where D is VPD, Vcmax should be replaced with Jmax if the photosynthesis is light limited, K is the Michaelis-Menten coefficient, Γ* is the CO2 compensation point, and η* is the ratio of water viscosity under ambient air temperature to a reference temperature at 298.15 K. The full analytical form of the partial derivatives can be retrieved from the MATLAB script in the SI Appendix.
To compare with other studies, we compute the changes in GPP relativized by changes in CO2 with a metric used in Ref. 8:
where GPPs and ca,s are of values in the starting year and GPPe and ca,e are of values in the ending year. A βln value of a unity indicates direct proportionality GPP response to CO2. We remind that this βln factor is only for intercomparison purposes; it is not a partial derivative and does not, by default, isolate the response of GPP to factors other than CO2 as βCO2 does. Since our framework can analytically attribute the trends to individual factors, we distinguish direct CO2 responses, βdirln, from apparent "CO2 responses," βappln.

FLUXNET Dataset.

In this study, site-scale EC-inferred GPP and meteorological forcings are provided by the FLUXNET2015 Tier 1 dataset, and the uncertainties associated with this dataset have been previously reviewed (14). To calculate the EEO-inferred GPP (and its sensitivities to various factors), we use air temperature “TA_F,” incoming shortwave radiation “SW_IN_F,” surface pressure “PA_F,” VPD “VPD_F” (convertible to specific humidity), and volumetric SWC “SWC_F_MDS_1” in the dataset as inputs. The EEO framework performs calculations on a monthly basis for daytime averages.
Filtering of meteorological data. We only use data in which air temperature, incoming shortwave radiation, and VPD are measured or gap filled with good quality [i.e., quality flag (QF) is 0 or 1]. The same filtering is also applied to the EC-inferred GPP (i.e., FLUXNET2015 GPP).
Calculation of monthly daytime averages. We first compute the monthly average of the variables at a given half-hour/hour and aggregate them to monthly daytime averages. Daytime is determined by the nighttime flag (QF = 0) provided by the dataset and further checked with the half-hourly/hourly shortwave radiation (>0 Wm−2).
Calculation of reference EC-inferred GPP. In order to compare our results with EC-inferred GPP, we use the mean of “GPP_NT_VUT_MEAN,” “GPP_DT_VUT_MEAN,” “GPP_NT_CUT_MEAN,” and “GPP_DT_CUT_MEAN” as the main reference of EC-inferred GPP. In addition, we also use the individual EC-inferred GPP based on different partitioning and friction velocity filtering methods to analyze the corresponding GPP uncertainty propagation (shaded area in Fig. 1A).
Exclusion of low-quality monthly data. We exclude all the monthly averages with less than 50% good-quality timestamps or negative GPP.
Aggregation of annual GPP. Since the quality filtering from steps 1 to 4 may lead to data gaps, in annual aggregations, we ensure that no sites have more than one missing value of EEO-inferred GPP and EC-inferred GPP during the growing season (defined as monthly climatological mean EC-inferred GPP > 30 gC m−2 mo−1). If there is one and only one such missing value during a year, EEO-inferred and EC-inferred GPP will be filled by their corresponding low-quality data without the filtering of step 4. Otherwise, that particular year will be excluded.
Calculation of GPP annual anomalies. We compute the annual anomalies of the EEO-inferred GPP and EC-inferred GPP for each site independently across the observational network (SI Appendix, Fig. S1) and aggregate the anomalies across sites (Fig. 1A).
Last, our evaluations are performed for those sites with more than 5 y of good-quality data at the annual scale.
After the preprocessing, there are 68 qualified flux sites (632 site-years) from 2001 to 2014, covering 10 different biome types, including croplands (10), closed shrublands (1), deciduous broadleaf forests (11), EBFs (4), evergreen needleleaf forests (18), grasslands (12), mixed forest (6), open shrublands (2), savannas (1), and woody savannas (3) (SI Appendix, Table S1). In the site-scale analysis, all biome types are treated as C3 species because of lack of information. The GPP trends of site-scale analysis are examined by the Mann-Kendall test (R package: An overall uncertainty of the site-scale GPP trends is characterized by the IQR of trends through resampling of site-years (Fig. 2A; SI Appendix, Fig. S3D). In addition, we calculated the IQRs due to site heterogeneity by resampling of sites (SI Appendix, Fig. S3 B and E) and due to the uneven distribution of sites and site-years over the study period by randomly shuffling the first-year data (SI Appendix, Fig. S3 C and F). The shaded area in Fig. 1A shows the variations in the EC-inferred GPP due to the partitioning and frictional velocity filtering methods and their propagation to the EEO-inferred GPP. On the other hand, we analyze the dataset without the quality filtering in steps 1, 2, and 4. This allows a comparison of our EEO-inferred GPP with the fully gap-filled FLUXNET data (SI Appendix, Fig. S3A), ensuring that the resulting findings are not due to data filtering.

Global Meteorological Forcing Datasets.

Our main global-scale results are computed using the ERA5-land hourly data (2001 to 2018, 0.1° × 0.1°,!/dataset/10.24381/cds.e2161bac). We use the following variables: “2m temperature,” “surface solar radiation downward,” “surface pressure,” “2m dew point temperature” (converted to specific humidity), and “volumetric soil water layers 1 and 2” (0 to 28 cm). We convert them to a monthly daytime average by examining the downward surface solar radiation (>20 W m−2) and resampling to 0.5° × 0.5°.
We estimate another EEO-inferred GPP using CRU-JRA55 forcing (2001 to 2018, 0.5° × 0.5°, for comparison. We use the following variables: “temperature at 2m,” “maximum temperature at 2m,” “downward solar radiation flux,” “pressure,” and “specific humidity.” Since there is no volumetric SWC in CRU-JRA55, we adopt it from ERA5-land. Daytime temperature is averaged from “temperature at 2m” and “maximum temperature at 2m.”

MODIS Leaf Area Index.

We use Collection 6 MODIS LAI as our LAI inputs (52). We refine the MODIS Terra and Aqua LAI (8-d frequency, 500 m) by checking their quality flags following a previous study (51). The quality of this dataset has been extensively validated (55) and reported elsewhere (51, 79). We convert the original 8-d LAI into monthly frequencies. For the flux-site analysis, we match the 500-m LAI to each site’s longitude and latitude and average over in a 3 × 3 window (1.5 km × 1.5 km). For the global-scale analysis, we convert the 500-m LAI into 0.5° × 0.5°.

Atmospheric CO2.

We use the monthly average atmospheric CO2 concentration from the Mauna Loa Observatory and the South Pole Observatory provided by NOAA's Earth System Research Laboratory (

DGVM Gross Primary Production.

DGVM outputs from the Trends in Net Land-Atmosphere Exchange (TRENDY-v6) project are used in this study ( (1). Here, we use two model experiments from TRENDY-v6, either varying CO2 only (time-invariant “preindustrial” climate and land use mask, S1) or varying CO2, climate, and land use (all forcing time-varying, S3). We exclude Sheffield DGVM (SDGVM) as it has an unrealistic sudden drop in GPP after 2007. More details are documented in SI Appendix, Table S4.

Satellite-Derived GPP Products.

The following GPP products are used: 1) Boreal Ecosystem Productivity Simulator (80), 2) Breathing Earth System Simulator (81), 3) FluxCom (82), 4) and 5) MODIS Collection 5.5 (MOD-C55) and Collection 6 (MOD-C6) (83), 6) photosynthesis model (Pmodel-s0) (44), 7) photosynthesis–respiration model (PRmodel) (7), and 8) VPM (84). The ensemble global mean GPP of all these eight products through 2001 to 2016 is about 119 PgC y−1. More details about these products are documented in SI Appendix, Table S4.

Land Cover Maps.

In order to analyze the spatial pattern of our diagnosed results, we use Collection 6 MODIS yearly product as the reference land cover map (2001 to 2018, 0.05° × 0.05°) (85). We refine the yearly maps for the study period into a single map and resample to 0.5° × 0.5° by taking the mode class of each grid cell. We aggregate the International Geosphere-Biosphere Program classification types into five C3 biomes (EBFs, other forests, short woody vegetation, grasslands, croplands) and one C4 biome (SI Appendix, Fig. S11). There is no biome aggregation for EBFs. Other forests are aggregated from deciduous broadleaf forests, mixed forests, evergreen needleleaf forests, deciduous needleleaf forests, and woody savannas. Short woody vegetation includes closed shrublands, open shrublands, savannas, and permanent wetlands. Croplands include croplands and cropland/natural vegetation mosaics. Finally, grids with greater than 50% C4 vegetation are labeled as C4 according to the ISLSCP II C4 vegetation percentage (

Data Availability

The following data were used in this work: FLUXNET2015:; MODIS LAI: and; ERA5-land:!/dataset/10.24381/cds.e2161bac; CRU-JRA55:; TRENDY-v6:; MODIS land cover:; ISLSCP II C4 vegetation percentage:; NOAA Atmospheric CO2:; BESS GPP:; MOD-C55 GPP:; MOD-C6 GPP:; FluxCom GPP:; Pmodel-s0:; VPM GPP: All other study data are referenced in the article. Codes are included in the supporting information.


C.C., T.F.K., and W.J.R. are supported by the US Department of Energy Office of Science Biological and Environmental Research, Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation Scientific Focus Area. T.F.K. acknowledges support from NASA Terrestrial Ecology Program Awards NNH17AE86I and 80NSSC21K1705. T.F.K. and I.C.P. acknowledge additional support from the LEMONTREE (Land Ecosystem Models based On New Theory, obseRvations and ExperimEnts) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures programme.

Supporting Information

Materials/Methods, Supplementary Text, Tables, Figures, and/or References

Appendix 01 (PDF)
Dataset S01 (TXT)


P. Friedlingstein et al., Global carbon budget 2020. Earth Syst. Sci. Data 12, 3269–3340 (2020).
B. E. Medlyn et al., Effects of elevated [CO2] on photosynthesis in European forest species: A meta-analysis of model parameters. Plant Cell Environ. 22, 1475–1495 (2002).
E. A. Ainsworth, S. P. Long, What have we learned from 15 years of free-air CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2. New Phytol. 165, 351–371 (2005).
E. A. Ainsworth, A. Rogers, The response of photosynthesis and stomatal conductance to rising [CO2]: Mechanisms and environmental interactions. Plant Cell Environ. 30, 258–270 (2007).
P. J. Franks et al., Sensitivity of plants to changing atmospheric CO2 concentration: From the geological past to the next century. New Phytol. 197, 1077–1094 (2013).
D. Schimel, B. B. Stephens, J. B. Fisher, Effect of increasing CO2 on the terrestrial carbon cycle. Proc. Natl. Acad. Sci. U.S.A. 112, 436–441 (2015).
T. F. Keenan et al., Recent pause in the growth rate of atmospheric CO2 due to enhanced terrestrial carbon uptake. Nat. Commun. 7, 13428–10 (2016).
A. P. Walker et al., Integrating the evidence for a terrestrial carbon sink caused by increasing atmospheric CO2. New Phytol. 229, 2413–2445 (2021).
G. G. Katul, R. Oren, S. Manzoni, C. Higgins, M. B. Parlange, Evapotranspiration: A process driving mass transport and energy exchange in the soil‐plant‐atmosphere‐climate system. Rev. Geophys. 50, 1083–25 (2012).
C. Chen et al., Biophysical impacts of Earth greening largely controlled by aerodynamic resistance. Sci. Adv. 6, eabb1981 (2020).
Z. Zhu et al., Greening of the earth and its drivers. Nat. Clim. Chang. 6, 791–795 (2016).
R. J. Norby, D. R. Zak, Ecological lessons from free-air CO2 enrichment (FACE) experiments. Ecol Evol Syst 42, 181–203 (2011).
P. B. Reich, S. E. Hobbie, T. D. Lee, M. A. Pastore, Unexpected reversal of C3 versus C4 grass response to elevated CO2 during a 20-year field experiment. Science 360, 317–320 (2018).
G. Pastorello et al., The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. Sci. Data 7, 225 (2020).
M. Fernández-Martínez et al., Global trends in carbon sinks and their relationships with CO2 and temperature. Nat. Clim. Chang. 9, 1–9 (2018).
S. Wang et al., Recent global decline of CO2 fertilization effects on vegetation photosynthesis. Science 370, 1295–1300 (2020).
M. Ueyama et al., Inferring CO2 fertilization effect based on global monitoring land-atmosphere exchange with a theoretical model. Environ. Res. Lett. 15, 084009–084012 (2020).
W. Cai, I. C. Prentice, Recent trends in gross primary production and their drivers: Analysis and modelling at flux-site and global scales. Environ. Res. Lett. 15, 124050 (2020).
W. K. Smith et al., Large divergence of satellite and Earth system model estimates of global terrestrial CO2 fertilization. Nat. Clim. Chang. 6, 306–310 (2015).
M. G. D. Kauwe, T. F. Keenan, B. E. Medlyn, I. C. Prentice, C. Terrer, Satellite based estimates underestimate the effect of CO2 fertilization on net primary productivity. Nat. Clim. Chang. 6, 892–893 (2016).
S. Wenzel, P. M. Cox, V. Eyring, P. Friedlingstein, Projected land photosynthesis constrained by changes in the seasonal cycle of atmospheric CO2. Nature 538, 499–501 (2016).
A. J. Winkler, R. B. Myneni, G. A. Alexandrov, V. Brovkin, Earth system models underestimate carbon fixation by plants in the high latitudes. Nat. Commun. 10, 885 (2019).
T. F. Keenan et al., Terrestrial biosphere model performance for inter‐annual variability of land‐atmosphere CO2 exchange. Glob. Change Biol. 18, 1971–1987 (2012).
A. Anav et al., Spatiotemporal patterns of terrestrial gross primary production: A review. Rev. Geophys. 53, 785–818 (2015).
Y. Ryu, J. A. Berry, D. D. Baldocchi, What is global photosynthesis? History, uncertainties and opportunities. Remote Sens. Environ. 223, 95–114 (2019).
A. Rogers et al., A roadmap for improving the representation of photosynthesis in Earth system models. New Phytol. 213, 22–42 (2017).
B. E. Medlyn et al., Using ecosystem experiments to improve vegetation models. Nat. Clim. Chang. 5, 528–534 (2015).
V. Haverd et al., Higher than expected CO2 fertilization inferred from leaf to global observations. Glob. Change Biol. 26, 2390–2402 (2020).
G. D. Farquhar, S. von Caemmerer, J. A. Berry, A biochemical model of photosynthetic CO2 assimilation in leaves of C 3 species. Planta 149, 78–90 (1980).
G. Katul, S. Manzoni, S. Palmroth, R. Oren, A stomatal optimization theory to describe the effects of atmospheric CO2 on leaf photosynthesis and transpiration. Ann. Bot. 105, 431–442 (2010).
T. F. Keenan et al., Increase in forest water-use efficiency as atmospheric carbon dioxide concentrations rise. Nature 499, 324–327 (2013).
I. C. Prentice, N. Dong, S. M. Gleason, V. Maire, I. J. Wright, Balancing the costs of carbon gain and water transport: Testing a new theoretical framework for plant functional ecology. Ecol. Lett. 17, 82–91 (2014).
S. P. Harrison et al., Eco-evolutionary optimality as a means to improve vegetation and land-surface models. New Phytol. 231, 2125–2141 (2021).
H. Wang et al., Towards a universal model for carbon dioxide uptake by plants. Nat. Plants 3, 734–741 (2017).
N. G. Smith et al., Global photosynthetic capacity is optimized to the environment. Ecol. Lett. 22, 506–517 (2019).
J.-L. Chen, J. F. Reynolds, P. C. Harley, J. D. Tenhunen, Coordination theory of leaf nitrogen distribution in a canopy. Oecologia 93, 63–69 (1993).
V. Maire et al., The coordination of leaf photosynthesis links C and N fluxes in C3 plant species. PLoS One 7, e38345–e15 (2012).
J. Kattge, W. Knorr, Temperature acclimation in a biochemical model of photosynthesis: A reanalysis of data from 36 species. Plant Cell Environ. 30, 1176–1190 (2007).
J. E. Campbell et al., Large historical growth in global terrestrial gross primary production. Nature 544, 84–87 (2017).
P. Ciais et al., Large inert carbon pool in the terrestrial biosphere during the Last Glacial Maximum. Nat. Geosci. 5, 74–79 (2011).
L. Cheng et al., Recent increases in terrestrial carbon uptake at little cost to the water cycle. Nat. Commun. 8, 110 (2017).
I. Ehlers et al., Detecting long-term metabolic shifts using isotopomers: CO2-driven suppression of photorespiration in C3 plants over the 20th century. Proc. Natl. Acad. Sci. U.S.A. 112, 15585–15590 (2015).
C. A. Gunderson et al., Environmental and stomatal control of photosynthetic enhancement in the canopy of a sweetgum (Liquidambar styraciflua L.) plantation during 3 years of CO2 enrichment. Plant Cell Environ. 25, 379–393 (2002).
B. D. Stocker et al., Drought impacts on terrestrial primary production underestimated by satellite monitoring. Nat. Geosci. 12, 264–270 (2019).
P. B. Reich, S. E. Hobbie, T. D. Lee, Plant growth enhancement by elevated CO2 eliminated by joint water and nitrogen limitation. Nat. Geosci. 7, 920–924 (2014).
A. G. Konings, A. P. Williams, P. Gentine, Sensitivity of grassland productivity to aridity controlled by stomatal and xylem regulation. Nat. Geosci. 10, 284–288 (2017).
V. Humphrey et al., Soil moisture-atmosphere feedback dominates land carbon uptake variability. Nature 592, 65–69 (2021).
B. E. Medlyn et al., Reconciling the optimal and empirical approaches to modelling stomatal conductance. Glob. Change Biol. 17, 2134–2144 (2011).
S. Manzoni et al., Optimizing stomatal conductance for maximum carbon gain under water stress: A meta-analysis across plant functional types and climates. Funct. Ecol. 25, 456–467 (2011).
Y. Zhang, C. Song, L. E. Band, G. Sun, No proportional increase of terrestrial gross carbon sequestration from the greening earth. J. Geophys. Res. Biogeosci. 124, 2540–2553 (2019).
C. Chen et al., China and India lead in greening of the world through land-use management. Nat. Sustain. 2, 122–129 (2019).
R. B. Myneni et al., Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sens. Environ. 83, 214–231 (2002).
H. Hersbach et al., The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146, 1999–2049 (2020).
J. Muñoz-Sabater et al., ERA5-land: A state-of-the-art global reanalysis dataset for land applications. Earth Syst Sci Data Discuss 2021, 1–50 (2021).
K. Yan et al., Evaluation of MODIS LAI/FPAR product collection 6. Part 2: Validation and intercomparison. Remote Sens. 8, 460–26 (2016).
M. Zhao, S. W. Running, Drought-induced reduction in global terrestrial net primary production from 2000 through 2009. Science 329, 940–943 (2010).
A. Samanta et al., Comment on “Drought-induced reduction in global terrestrial net primary production from 2000 through 2009”. Science 333, 1093, author reply 1093 (2011).
B. E. Medlyn, Comment on “Drought-induced reduction in global terrestrial net primary production from 2000 through 2009”. Science 333, 1093, author reply 1093 (2011).
M. Zhao, S. W. Running, Response to comments on “Drought-induced reduction in global terrestrial net primary production from 2000 through 2009”. Science 333, 1093 (2011).
C. Chen et al., Prototyping of LAI and FPAR retrievals from MODIS multi-angle implementation of atmospheric correction (MAIAC) data. Remote Sens. 9, 370 (2017).
J. Peñuelas, J. G. Canadell, R. Ogaya, Increased water‐use efficiency during the 20th century did not translate into enhanced tree growth. Glob. Ecol. Biogeogr. 20, 597–608 (2011).
P. van der Sleen et al., No growth stimulation of tropical trees by 150 years of CO2 fertilization but water-use efficiency increased. Nat. Geosci. 8, 24–28 (2015).
K. Huang et al., Enhanced peak growth of global vegetation and its key mechanisms. Nat. Ecol. Evol. 2, 1897–1905 (2019).
Z. Zhu et al., Comment on “Recent global decline of CO2 fertilization effects on vegetation photosynthesis”. Science 373, eabg5673 (2021).
C. Frankenberg, Y. Yin, B. Byrne, L. He, P. Gentine, Comment on “Recent global decline of CO2 fertilization effects on vegetation photosynthesis”. Science 373, eabg2947 (2021).
Y. Sang et al., Comment on “Recent global decline of CO2 fertilization effects on vegetation photosynthesis”. Science 373, eabg4420 (2021).
T. F. Keenan, Ü. Niinemets, Global leaf trait estimates biased due to plasticity in the shade. Nat. Plants 3, 16201 (2016).
Ü. Niinemets, T. F. Keenan, L. Hallik, A worldwide analysis of within-canopy variations in leaf structural, chemical and physiological traits across plant functional types. New Phytol. 205, 973–993 (2015).
P. Hari, A. Mäkelä, E. Korpilahti, M. Holmberg, Optimal control of gas exchange. Tree Physiol. 2, 169–175 (1986).
G. G. Katul, S. Palmroth, R. Oren, Leaf stomatal responses to vapour pressure deficit under current and CO(2)-enriched atmosphere explained by the economics of gas exchange. Plant Cell Environ. 32, 968–979 (2009).
T. N. Buckley, K. A. Mott, Modelling stomatal conductance in response to environmental factors. Plant Cell Environ. 36, 1691–1699 (2013).
Y. Wang, J. S. Sperry, W. R. L. Anderegg, M. D. Venturas, A. T. Trugman, A theoretical and empirical assessment of stomatal optimization modeling. New Phytol. 227, 311–325 (2020).
S. von Caemmerer, G. D. Farquhar, Some relationships between the biochemistry of photosynthesis and the gas exchange of leaves. Planta 153, 376–387 (1981).
B. E. Medlyn, R. A. Duursma, M. G. D. Kauwe, I. C. Prentice, The optimal stomatal response to atmospheric CO2 concentration: Alternative solutions, alternative interpretations. Agric. For. Meteorol. 182, 200–203 (2013).
K. A. Mott, Do stomata respond to CO2 concentrations other than intercellular? Plant Physiol. 86, 200–203 (1988).
S. M. Assmann, The cellular basis of guard cell sensing of rising CO2. Plant Cell Environ. 22, 629–637 (1999).
Y.-S. Lin et al., Optimal stomatal behaviour around the world. Nat. Clim. Chang. 5, 459–464 (2015).
C. Field, H. A. Mooney, “The photosynthesis-nitrogen relationship in wild plants” in On the Economy of Plant Form and Function, T. J. Givnish, Ed. (Cambridge University Press, 1986), pp. 25–56.
S. Piao et al., Characteristics, drivers and feedbacks of global greening. Nat. Rev. Earth Environ. 1, 14–27 (2020).
J. M. Chen et al., Vegetation structural change since 1981 significantly enhanced the terrestrial carbon sink. Nat. Commun. 10, 4259 (2019).
C. Jiang, Y. Ryu, Multi-scale evaluation of global gross primary productivity and evapotranspiration products derived from Breathing Earth System Simulator (BESS). Remote Sens. Environ. 186, 528–547 (2016).
M. Jung et al., Scaling carbon fluxes from eddy covariance sites to globe: Synthesis and evaluation of the FLUXCOM approach. Biogeosciences 17, 1343–1365 (2020).
S. W. Running et al., A continuous satellite-derived measure of global terrestrial primary production. Bioscience 54, 547–560 (2004).
Y. Zhang et al., A global moderate resolution dataset of gross primary production of vegetation for 2000-2016. Sci. Data 4, 170165 (2017).
M. A. Friedl, D. K. McIver, J. Hodges, X. Y. Zhang, Global land cover mapping from MODIS: Algorithms and early results. Remote Sens. Environ. 83, 287–302 (2002).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 119 | No. 10
March 8, 2022
PubMed: 35238668


Data Availability

The following data were used in this work: FLUXNET2015:; MODIS LAI: and; ERA5-land:!/dataset/10.24381/cds.e2161bac; CRU-JRA55:; TRENDY-v6:; MODIS land cover:; ISLSCP II C4 vegetation percentage:; NOAA Atmospheric CO2:; BESS GPP:; MOD-C55 GPP:; MOD-C6 GPP:; FluxCom GPP:; Pmodel-s0:; VPM GPP: All other study data are referenced in the article. Codes are included in the supporting information.

Submission history

Received: August 24, 2021
Accepted: January 12, 2022
Published online: March 1, 2022
Published in issue: March 8, 2022


  1. CO2 fertilization effect
  2. photosynthesis
  3. GPP
  4. optimization theory
  5. carbon and water coupling


C.C., T.F.K., and W.J.R. are supported by the US Department of Energy Office of Science Biological and Environmental Research, Reducing Uncertainties in Biogeochemical Interactions through Synthesis and Computation Scientific Focus Area. T.F.K. acknowledges support from NASA Terrestrial Ecology Program Awards NNH17AE86I and 80NSSC21K1705. T.F.K. and I.C.P. acknowledge additional support from the LEMONTREE (Land Ecosystem Models based On New Theory, obseRvations and ExperimEnts) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures programme.


This article is a PNAS Direct Submission.



Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720
Department of Environmental Science, Policy and Management, University of California, Berkeley, CA 94720
Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720
I. Colin Prentice
Department of Life Sciences, Imperial College London, Ascot SL5 7PY, United Kingdom
Trevor F. Keenan1 [email protected]
Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720
Department of Environmental Science, Policy and Management, University of California, Berkeley, CA 94720


To whom correspondence may be addressed. Email: [email protected] or [email protected].
Author contributions: C.C. and T.F.K. designed research; C.C. performed research and analyzed data; and C.C., W.J.R., I.C.P., and T.F.K interpreted the results and wrote the paper.

Competing Interests

The authors declare no competing interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    CO2 fertilization of terrestrial photosynthesis inferred from site to global scales
    Proceedings of the National Academy of Sciences
    • Vol. 119
    • No. 10







    Share article link

    Share on social media