Rising CO2 and warming reduce global canopy demand for nitrogen.

Nitrogen (N) limitation has been considered as a constraint on terrestrial carbon uptake in response to rising CO2 and climate change. By extension, it has been suggested that declining carboxylation capacity (Vcmax ) and leaf N content in enhanced-CO2 experiments and satellite records signify increasing N limitation of primary production. We predicted Vcmax using the coordination hypothesis and estimated changes in leaf-level photosynthetic N for 1982-2016 assuming proportionality with leaf-level Vcmax at 25°C. The whole-canopy photosynthetic N was derived using satellite-based leaf area index (LAI) data and an empirical extinction coefficient for Vcmax , and converted to annual N demand using estimated leaf turnover times. The predicted spatial pattern of Vcmax shares key features with an independent reconstruction from remotely sensed leaf chlorophyll content. Predicted leaf photosynthetic N declined by 0.27% yr-1 , while observed leaf (total) N declined by 0.2-0.25% yr-1 . Predicted global canopy N (and N demand) declined from 1996 onwards, despite increasing LAI. Leaf-level responses to rising CO2 , and to a lesser extent temperature, may have reduced the canopy requirement for N by more than rising LAI has increased it. This finding provides an alternative explanation for declining leaf N that does not depend on increasing N limitation.

Introduction

Atmospheric CO2 has increased from c. 280 to nearly 420 μmol mol–1 (in 2022) since the Industrial Revolution. This increase is smaller than expected from the cumulative anthropogenic emissions of CO2. The difference is partly due to terrestrial ecosystems, which have taken up around one‐third of these emissions (Friedlingstein et al., 2021). There is concern, however, that this ‘terrestrial carbon sink’ may not persist due to the increasing limitation of primary production by nutrient availability (Reich et al., 2006). Much of this literature has focused on nitrogen (N) as a limiting nutrient, hence our focus on N here – although we recognize that other plant nutrients (notably phosphorus) may be more important in constraining primary production in the tropics, and in regions with high rates of N deposition.

The incorporation of carbon–nitrogen (C–N) cycle coupling has become a priority in the development of dynamic global vegetation models (DGVMs), partly in order to address concerns about the persistence of the terrestrial carbon sink. Yet, there is little consensus as to how C–N coupling should be represented in DGVMs. The mechanisms by which plants acquire additional N in response to elevated CO2 (e.g. Finzi et al., 2007) are often not well represented in DGVMs (Zaehle et al., 2014). Meanwhile, satellite‐derived records of global changes in chlorophyll content, and independent measurements of leaf N content, have been used to quantify trends in ecosystem C–N coupling (He et al., 2017; Penuelas et al., 2020 and references therein). It has been found that leaf N is declining and inferred that N limitation (i.e. an excess of N demand over supply) is increasing. Moreover, increasing N limitation has been taken to imply a reduction in the CO2 fertilization effect on gross primary production (GPP). As this effect is considered to be a significant mechanism driving the terrestrial carbon sink (Walker et al., 2021), increasing N limitation has been taken to imply a potential reduction in the efficiency of this sink (He et al., 2017; Wang et al., 2020) – albeit that other mechanisms, whose relative magnitudes are still uncertain, may also contribute (Walker et al., 2021).

Here, we interrogate the idea of increasing N limitation by considering how leaf photosynthetic N might be expected to change with rising CO2 and warming, based on leaf‐level optimality principles. We also consider how changes in leaf photosynthetic N might translate into changes in the whole‐canopy photosynthetic N (given observed increases in leaf area index (LAI)) and the annual canopy demand for N (which also depends on the rate at which leaves are replaced, and the efficiency with which the N in senescing leaves is resorbed). Studies have shown that leaf N is influenced by environmental factors during growth, including temperature and CO2 (Ainsworth & Long, 2005; Dong et al., 2017; Smith & Keenan, 2020). Many of these responses are predicted by the coordination hypothesis, which states that the Rubisco‐ and electron transport‐limited rates of photosynthesis tend to be colimiting under typical daytime conditions (Chen et al., 1993; Dong et al., 2017). According to this hypothesis, leaf‐level V cmax should acclimate such that currently available light is fully utilized. A growing body of both experimental (Scafaro et al., 2017) and observational (Dong et al., 2017; Smith et al., 2019) evidence supports the predictions of this hypothesis with respect to both experimentally imposed changes and spatial gradients in plants’ light environment, growth temperature and atmospheric dryness vapour pressure deficit (VPD). V cmax when standardized to a common temperature (typically 25°C, written here as V cmax25 and taken to be proportional to Rubisco content) is predicted to increase in proportion to absorbed light; to decline with increasing growth temperature, even as V cmax measured at growth temperature increases; and to increase with VPD, compensating for greater stomatal closure. V cmax at growth temperature increases to counteract increased photorespiration at higher temperatures (Smith et al., 2019), while V cmax25 declines because less Rubisco is required (Scafaro et al., 2017) to achieve a given photosynthetic rate – this second effect being stronger than the first (Fig. 1).

Fig. 1

(a) Ratio of V cmax (μmol m–2 s–1) to the mean V cmax over a range of growth CO2 levels. Red dots and bars, experimental data on soybean (Zheng et al., 2019), showing means and 95% confidence intervals (n = 5 for each treatment). Black dots, predicted optimal V cmax at growth temperature 25°C, photosynthetic photon flux density (PPFD) 1000 μmol m–2 s–1 and vapour pressure deficit (VPD) 1.5 kPa. (b) Predicted variations of V cmax and (c) V cmax25 with 400 ppm CO2 and PPFD and VPD as in (a), as a function of growth temperature. Both quantities are shown as ratios to their values at 25°C.

Free air carbon enrichment (FACE) experiments (Ainsworth & Long, 2005), and comparisons between plants growing in natural CO2 springs and in control plots (Saban et al., 2019), have shown that leaf‐level photosynthetic capacity (V cmax) and leaf N generally decline with enhanced CO2. This decline is predicted by the coordination hypothesis due to the different sensitivities of Rubisco‐ and electron transport‐limited photosynthesis to CO2. The Rubisco‐limited rate is more sensitive, so – all else equal – the optimal Rubisco content of leaves should decline as CO2 increases; otherwise, plants would have to bear the cost of maintaining Rubisco in excess of the amount they need in order to use the available light. Fig. 1 shows that the coordination hypothesis predicts the generally decreasing trend of experimentally induced reductions in V cmax as growth CO2 is increased (Zheng et al., 2019). Although these data show unexplained deviations from a smooth curve, the comparison (for soybean) is apposite because photosynthesis in this N‐fixing plant is unlikely to be limited by N supply. Yet, declining V cmax and leaf N with increasing CO2 are also observed in C3 plants more generally, and the coordination hypothesis supplies a simple and parsimonious explanation (Smith & Keenan, 2020).

In this Viewpoint article, we have used the coordination hypothesis to predict V cmax, yielding a global spatial pattern similar to that presented by Smith et al. (2019), and compared this predicted spatial pattern with a recently developed, remotely sensed map of leaf‐level chlorophyll content (LCC). LCC has been shown to be tightly linked to V cmax (Luo et al., 2019, 2021). We then predicted temporal variation in V cmax at growth temperature from the beginning of the satellite era (1982–2016) using the coordination hypothesis; converted predicted V cmax to V cmax25 using the Arrhenius equation; and converted predicted changes in V cmax25 to changes in leaf‐level photosynthetic (Rubisco) N. We further estimated changes in the whole‐canopy photosynthetic N, taking into account satellite‐derived changes in LAI, with the help of an empirical extinction coefficient describing the decline in leaf‐level V cmax with canopy depth. Finally, we estimated how the annual leaf‐ and canopy‐level demand for N uptake has changed over time, with the help of leaf turnover times separately estimated for different phenological types (evergreen vs deciduous plants).

Methods

LCC has been estimated using radiative transfer theory to downscale remotely sensed canopy reflectance to leaf reflectance, as input to a leaf optical model (Luo et al., 2019). Croft et al. (2020) thereby retrieved global monthly LCC from land surface spectral reflectance data provided by the Medium Resolution Imaging Spectrometer (MERIS) instrument aboard Envisat. MERIS provided red, near‐infrared and red‐edge bands at 300‐m spatial resolution for 7‐d intervals during 2003–2012. The correlation between chlorophyll and V cmax was used to map global V cmax based on LCC, using empirical relationships between V cmax25 and LCC for different plant functional types (Luo et al., 2019). However, whereas Luo et al. (2019) labelled their mapped quantity as V cmax25, our analysis presupposes that LCC is more closely related to V cmax at growing‐season temperature than to V cmax25. The empirical relationships used to relate V cmax and LCC were based on measurements made in a narrow range (17–20°C) before adjustment to 25°C (Teklemariam et al., 2009), so relationships of V cmax or V cmax25 to LCC could not be distinguished in those measurements.

Optimal V cmax values according to the coordination hypothesis were estimated using the following expression (Dong et al., 2017) for C3 plants acclimated to their growth conditions: where φ 0 is the intrinsic quantum efficiency of photosynthesis (mol mol–1), I 0 is the incident photosynthetic photon flux density (PPFD) (μmol m–2 s–1), c i is the intercellular partial pressure of CO2 (Pa) (for simplicity, and in line with most of the current plant functional ecology and carbon cycle modelling literature, we have disregarded the limitation imposed by finite mesophyll conductance), K is the effective Michaelis–Menten coefficient of Rubisco (Pa), and Γ* is the photorespiratory compensation point (Pa). Eqn 1 is a generic expression, independent of the precise formulation of J max limitation effects (models used by Wang et al. (2017) and Smith et al. (2019), for example, used slightly differing empirical formulations of the response of electron transport to light). Reference values and temperature dependencies of φ 0, K and Γ* in Eqn 1 were derived from Bernacchi et al. (2001). Values of c i in Eqn 1 were estimated by the least‐cost hypothesis as described in Prentice et al. (2014): where

Eqn (Eqn 2(a)), (Eqn 2(b)) yields the c i value (for a given ambient CO2 partial pressure, c a) that minimizes the sum of the costs (per unit assimilation) of maintaining the capacities for transpiration and carboxylation. The parameter ξ (Pa½) is inversely related to the sensitivity of stomata to VPD (D 0). The scaling factor β was estimated as 146 from a global analysis of leaf δ13C measurements (Wang et al., 2017). η* is the viscosity of water at growth temperature, relative to its value at 25°C. V cmax was converted to V cmax25 by inverse application of the Arrhenius equation with activation energy as given in Bernacchi et al. (2001).

Rubisco is the most abundant photosynthetic protein in the leaf. Photosynthetic proteins together account for 50–60% of leaf N, and Rubisco alone, 25–30%. We used the N in Rubisco as an index of total photosynthetic N (Evans & Seemann, 1989), implicitly assuming that other proteins – including those involved in electron transport and the Calvin cycle – covary in proportion with Rubisco. The leaf N required for Rubisco (N rubiscoL, g m–2 leaf area) was estimated from V cmax25 as follows (Dong et al., 2017): where M n = 14 g mol–1 is the molecular mass of N; M r = 0.55 g μmol–1 is the molecular mass of Rubisco; [N r] = 0.0114 mol g–1 is the N concentration of Rubisco; n = 8 is the number of catalytic sites per mole of Rubisco; and k cat = 3.5 s–1 is the catalytic turnover number at 25°C. To scale up leaf Rubisco N to the canopy, we used a big‐leaf approximation: where N rubiscoC (g m–2 ground area) is the canopy Rubisco N; L is the LAI; and k is the extinction coefficient for V cmax, assigned the value of 0.175 based on its median value in studies summarized by Lloyd et al. (2010). Eqn 4 accounts for the declining leaf‐level V cmax and Rubisco N with depth in the canopy, and results in a nonlinear response of Rubisco N to L. In our main analysis, L was estimated from monthly Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) LAI3g from 1982 to 2016. There are substantial differences in temporal variations among satellite‐based LAI products (Jiang et al., 2017). We therefore also tried other LAI products: TCDR LAI during 1982–2016; GLASS LAI during 1982–2014; and GLOBMAP LAI during 1982–2011 (Jiang et al., 2017). Differences among products are especially apparent before 2000, the first year after the launch of the Terra satellite (Supporting Information Fig. S1), and significant differences persist (Jiang et al., 2017; Cortés et al., 2021). However, with continuous improvement, most LAI products have become consistent in showing long‐term positive trends in many regions (Chen et al., 2019; Cortés et al., 2021).

We estimated the annual leaf‐level demand for photosynthetic N by: (with a parallel expression for the canopy level relating N demandC to N rubiscoC) where c is the N resorption efficiency, set at 0.5 for woody vegetation and zero for herbaceous vegetation (Aerts, 1996), and τ is the leaf replacement time (years). For deciduous vegetation, τ was set to 1 yr. For evergreen vegetation, τ = τ ev was estimated from the theoretical framework developed by Wang et al. (2021) for the leaf economics spectrum (Wright et al., 2004), which relates leaf mass per area (LMA, g m–2) to leaf lifespan: where Here, u = 768 (dimensionless) is a constant that relates the leaf ageing rate to LMA and photosynthetic capacity, and C = 23 (dimensionless) is a multiplier that accounts for the total costs of constructing leaves and other tissues required to support them. Both parameters were estimated from data by Wang et al. (2021). f is the growing‐season length, as a fraction of the year; k = 30 g (dry biomass) mol–1 (C) is a conversion factor between carbon and dry biomass; ΣI 0 is the daily total average photosynthetic photon flux density (PPFD) (mol m–2 d–1); and h T is the Arrhenius function for the response of V cmax to temperature, relative to its value at 25°C (Wang et al., 2021). LMA for evergreen plants was derived from a scaled‐up global trait data set (Fig. S2c). Herbaceous vegetation was treated as deciduous in climates with at least one month’s mean temperature < 0°C, otherwise as evergreen. This is a rough approximation that accounts for the fact that herbaceous biomass commonly turns over annually in cold‐winter climates. It is unrealistic for semi‐arid vegetation, but we assume the consequences are minor for global totals. Evergreen, deciduous and herbaceous vegetation fractions were assigned based on ESA CCI land cover products at 0.5° resolution (Li et al., 2018). The thermal growing season was defined as the period with daily mean temperatures (linearly interpolated between months) > 0°C. Fig. S2 shows values of τ and LMA used in these calculations.

Mean values of air temperature (T, °C), VPD (D 0, Pa) and total (ΣI 0, mol m–2 d–1) and average (I 0, μmol m–2 s–1) PPFD during the thermal growing season were calculated based on the CRU TS4.3 climate data for 1982–2016 at 0.5° resolution. Mean daytime air temperature was estimated from daily temperature maxima and minima by sinusoidal interpolation. Three historical simulations were performed: (1) all factors (CO2, climate) varied; (2) climate varied, with CO2 fixed at 340 ppm (its value in 1982); and (3) CO2 varied, with climate variables fixed at their mean values from 1982 to 2016. The Theil–Sen regression slopes were fitted to indicate the direction and magnitude of the trends. The Theil–Sen regression slope is the median slope of all straight lines joining pairs of data points and provides a robust estimate that is less sensitive to outliers than ordinary least‐squares linear regression. The Theil–Sen regression was implemented using the ‘spatialEco’ package in R. Segmented regressions (‘segmented’ package) were applied to assess the timing of breaks in the time series of predicted canopy‐level N demand based on each of the LAI data products. We also mapped the spatial pattern of the historical simulation based on CO2, climate and LAI. All analyses and graphics were developed in R.

Results

The global pattern of predicted V cmax at growing‐season temperature (Fig. 2a) shows good general agreement with satellite LCC‐derived V cmax (Fig. 2b) (r = 0.56, P < 0.001). There is underestimation in predominantly cropland areas in interior North America, Europe and East Asia. The LCC‐derived V cmax for croplands depends on a cropland‐specific conversion factor between LCC and V cmax; our model did not distinguish crops from other plants and therefore might underestimate V cmax in croplands (which can be influenced by fertilization and irrigation, as well as varietal selection). The spatial correlation increased to r = 0.63 after excluding croplands. There are a number of specific differences between the predicted and observed maps that we do not explore here. However, there are notable, large‐scale geographic features in common, including a belt of exceptionally high values in the Sahel (consistent with recent measurements by Sibret et al., 2021), high values in north‐western India and steep declines north of 50°N. All of these features were already shown in the global map of predicted optimal V cmax presented by Smith et al. (2019). Fig. 2 demonstrates that they are present in nature and observable from space.

Fig. 2

Global distributions of (a) satellite‐derived V cmax (V cmax,chl, μmol m–2 s–1) from leaf chlorophyll content and (b) predicted V cmax at growth temperature during the period from 2003 to 2012.

During the period from 1982 to 2016, CO2 increased by 58.5 ppm and global mean land temperature by 0.5°C. Predicted V cmax at growing‐season temperature generally increased, while V cmax25 decreased, in response to climate change (Fig. 3c,d). These responses are consistent with experimentally observed leaf‐level responses to increasing temperature (Scafaro et al., 2017). Both quantities declined in response to rising CO2 (Fig. 3e,f). The net effect of warming and CO2 increase was variable in sign for V cmax at growing‐season temperature, but showed an almost consistently negative trend for V cmax25 (Fig. 3a,b). The rate of decline in leaf‐level photosynthetic N over the period studied was c. 0.27% yr–1 (Fig. 4a). For comparison, measured rates of decline in total leaf N have been in the range from 0.2 to 0.25% yr–1 (Craine et al., 2018; Penuelas et al., 2020; Wang et al., 2020).

Fig. 3

Global temporal trends (Theil–Sen slope coefficients) in predicted V cmax (μmol m–2 s–1) at growing‐season temperature (left panel) and at standard temperature, 25°C (V cmax25, μmol m–2 s–1) during 1982–2016 for three scenarios: (a, b) all factors; (c, d) varying climate only, with CO2 fixed at 340 ppm; and (e, f) varying CO2 only, with climate fixed at multi‐year mean values. During this period, mean annual temperatures on land increased by 0.5°C.

Fig. 4

Time series of modelled global (a) leaf‐level N rubisco (g m–2 leaf area), declining by 0.27% yr–1 after 1996; (b) canopy‐level N rubisco (g m–2 ground area), declining by 0.36% yr–1 after 1996; (c) annual leaf‐level N rubisco demand (g m–2 leaf area yr–1), declining by 0.28% yr–1; and (d) annual canopy‐level N rubisco demand (g m–2 ground area yr–1), declining by 0.35% yr–1 after 1996. Leaf area index variations were based on the LAI3g data set. Black lines, all factors varied. Blue lines, varying climate only, with CO2 fixed at 340 ppm. Red lines, varying CO2 only, with climate variables fixed at their multi‐year mean values. The vertical grey lines in (b) and (d) indicate the breakpoint (1996) identified by segmented regression.

Predicted leaf‐level photosynthetic N was lowest in the tropics and highest at high latitudes and elevations (Fig. 5a). This first‐order pattern is predicted by the coordination hypothesis as a consequence of enzyme kinetics. As temperature increases, less Rubisco is required to match the light‐limited rate of photosynthesis. This gradient was partially offset by long leaf turnover times (Fig. S2b) in boreal evergreen forests, leading to a more even spatial distribution of the annual leaf‐level requirement for photosynthetic N (Fig. 5c). The estimated whole‐canopy photosynthetic N (Fig. 5b) was highest in forests, especially in the boreal zone. Similar to the pattern of leaf‐level N, this was offset by long leaf turnover times, but canopy‐level demand for photosynthetic N was nonetheless greatest in high latitudes (Fig. 5d).

Fig. 5

Predicted average global distributions of (a) leaf‐level N rubisco (g m–2 leaf area); (b) canopy‐level N rubisco (g m–2 ground area); (c) annual leaf‐level N rubisco demand (g m–2 leaf area yr–1); and (d) annual canopy level of N rubisco demand (g m–2 ground area yr–1) over the period from 1982 to 2016. The colour scale for all maps is in square root transformation.

Predicted trends with increasing CO2 in both leaf‐level photosynthetic N and N demand (Fig. 4a,c) were negative. At the canopy level, warming and (especially) increasing CO2 imply decreasing demand, while LAI increases imply increasing demand (Fig. 4b,d). Since the late 1990s, however, the modelled global trend in total canopy N (in CO2‐only simulations and all‐factor simulations) has been negative, despite widespread increases in LAI. Simulations based on alternative vegetation cover data sets showed greater variations, and differences from one another, in the early part of the record (Fig. S1). However from 1997, which was the estimated breakpoint in the LAI3g time series, all‐factor simulations using the four data sets consistently showed declining trends in the whole‐canopy photosynthetic N demand, with P < 0.01 for the all data sets except GLASS (P = 0.2) (Fig. S3).

Discussion

Despite many approximations, the general findings of our analysis are clear. Leaf photosynthetic N is expected to decline as CO2 and temperatures rise, on optimality grounds, and the magnitude of this effect is compatible with observed declines in total leaf N. The comparison is not exact because our calculations refer only to the photosynthetic component of leaf N. However, the magnitude is similar. Moreover, the effect was large enough to counteract the potential increase in canopy N demand implied by observed, increasing trends in LAI.

It follows that declining leaf N does not necessarily imply increasing N limitation of primary production. In other words, declining leaf N is the predicted outcome of an acclimation process, by which plants avoid incurring maintenance costs in excess of requirements set by their growth conditions.

The terrestrial carbon sink does not solely depend on rising GPP. Terrestrial carbon uptake may be limited by factors other than N, especially in the tropics and in regions subject to anthropogenic N deposition, and influenced by changes in the residence time of carbon in ecosystems. We are not arguing for the absence of constraints on the carbon sink, today or in future. Our analysis nonetheless suggests that declining leaf N should not, of itself, be taken to imply that N limitation of GPP has increased.

Rising CO2 and temperatures affect land ecosystems in multiple ways. Establishing the dominant processes is a challenge. This applies particularly to N‐cycle processes, because rates of N input (Wieder et al., 2015), turnover and loss via different pathways (Fang et al., 2015; Meyerholt et al., 2020) are all challenging to determine. Moreover, the biotic and environmental controls on more readily observable quantities, including leaf N, stable isotope ratios (δ15N) and nitrous oxide (N2O) emissions, are incompletely understood. Global environmental changes induce competing effects, including increased N mineralization due to warming and increased N demand due to greening. Greening has been variously attributed (depending on the region) to rising CO2, lengthening growing seasons and land‐use change (Zhu et al., 2016). In addition to these opposing forces, our results suggest that the increase in N demand due to greening has been outweighed in recent decades by a reduction in N demand due to the down‐regulation of leaf‐level photosynthetic Ν in response to rising CO2 and warming.

Several studies have shown regional or global declines in plant N isotope signatures (δ15N) over recent decades (McLauchlan et al., 2017; Craine et al., 2018). These too have been attributed to increasing N limitation of plant growth (Caldararu et al., 2022). However, plant δ15N is not an unambiguous indicator of N limitation. Soil (and therefore plant) δ15N values are influenced by the partitioning of ecosystem N losses between gaseous (N2, N2O, NO, NH3) emissions and leaching (Houlton & Bai, 2009). Effects of elevated CO2 on gaseous N emissions can be positive or negative, with increases in wetter environments and decreases in drier environments (Phillips et al., 2001). This pattern appears consistent with one shown by McLauchlan et al. (2017): warmer and drier forested regions of the USA showed increasing δ15N, in contrast with cooler and wetter regions.

Thus, although declining leaf N has been interpreted as an indicator of increasing N limitation on plant growth (He et al., 2017), this attribution neglects the role of acclimation processes that are predicted to reduce canopy demand for N. Our calculations here indicate that these processes are of comparable magnitude to observed trends in leaf and canopy N, and large enough to counteract the increase in demand caused by increases in LAI.

Author contributions

ND and ICP designed the study. ND carried out the data analysis and graphic and drafted the first version with a significant input from ICP. ND, ICP and IJW interpreted the analysis. JMC, XZL and TFK provided the data. HW provided help on a new LES model. NGS provided codes. All authors contributed to the subsequent manuscripts.

Fig. S1 Time series of modelled global leaf‐level N rubisco, canopy‐level N rubisco, annual leaf‐level N rubisco demand and annual canopy‐level N rubisco demand with four different LAI datasets.

Fig. S2 Maps of vegetation turnover rate (yr–1), evergreen longevity (yr) and leaf mass per area (g m–2).

Fig. S3 Bivariate plots of annual canopy‐level N rubisco demand against CO2 after 1996, based on different LAI datasets.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.