Stomatal conductance limited the CO2 response of grassland in the last century
BMC Biology volume 19, Article number: 50 (2021)
The anthropogenic increase of atmospheric CO2 concentration (ca) is impacting carbon (C), water, and nitrogen (N) cycles in grassland and other terrestrial biomes. Plant canopy stomatal conductance is a key player in these coupled cycles: it is a physiological control of vegetation water use efficiency (the ratio of C gain by photosynthesis to water loss by transpiration), and it responds to photosynthetic activity, which is influenced by vegetation N status. It is unknown if the ca-increase and climate change over the last century have already affected canopy stomatal conductance and its links with C and N processes in grassland.
Here, we assessed two independent proxies of (growing season-integrating canopy-scale) stomatal conductance changes over the last century: trends of δ18O in cellulose (δ18Ocellulose) in archived herbage from a wide range of grassland communities on the Park Grass Experiment at Rothamsted (U.K.) and changes of the ratio of yields to the CO2 concentration gradient between the atmosphere and the leaf internal gas space (ca – ci). The two proxies correlated closely (R2 = 0.70), in agreement with the hypothesis. In addition, the sensitivity of δ18Ocellulose changes to estimated stomatal conductance changes agreed broadly with published sensitivities across a range of contemporary field and controlled environment studies, further supporting the utility of δ18Ocellulose changes for historical reconstruction of stomatal conductance changes at Park Grass. Trends of δ18Ocellulose differed strongly between plots and indicated much greater reductions of stomatal conductance in grass-rich than dicot-rich communities. Reductions of stomatal conductance were connected with reductions of yield trends, nitrogen acquisition, and nitrogen nutrition index. Although all plots were nitrogen-limited or phosphorus- and nitrogen-co-limited to different degrees, long-term reductions of stomatal conductance were largely independent of fertilizer regimes and soil pH, except for nitrogen fertilizer supply which promoted the abundance of grasses.
Our data indicate that some types of temperate grassland may have attained saturation of C sink activity more than one century ago. Increasing N fertilizer supply may not be an effective climate change mitigation strategy in many grasslands, as it promotes the expansion of grasses at the disadvantage of the more CO2 responsive forbs and N-fixing legumes.
Atmospheric CO2 increase and related climate change are altering fundamentally the terrestrial biogeochemical cycles of carbon (C), water, and nitrogen (N) across grasslands, forests, and croplands [1,2,3,4,5]. At local scales, these cycles are linked by processes at plant surfaces of leaves [6,7,8,9,10,11,12] and roots [3, 13, 14], critical interfaces in the soil-plant-atmosphere continuum. The uptake of CO2 by photosynthesis and evaporative loss of water vapor in transpiration occur through the same path, the stomata . These small pores in the leaf epidermis adjust their conductance in response to environmental conditions [9, 10], such as atmospheric CO2 concentration (ca) [7, 11], and internal cues that include photosynthetic activity [7, 15,16,17] which correlates with leaf nitrogen status [7, 15, 16, 18]. When stomata open for CO2 uptake, they simultaneously expose the leaf internal moisture to a comparatively dry atmosphere. This leaf-to-air vapor pressure gradient drives transpiration . The ratio of canopy-scale photosynthesis (A) and transpiration (E) determines vegetation water use efficiency (W) (Eq. 1) [12, 19]. The control function of stomatal conductance, and its implication for changes of water use efficiency (W) in the face of increasing ca, becomes evident when A and E are expressed as the products of stomatal conductance (gs) and the concentration gradients of the respective gases [12, 19]:
where ca and va are the CO2 and water vapor mole fractions in air, ci and vi that in the substomatal cavity, and 1.6 is the ratio of the diffusivities of CO2 and water vapor in air. If transpiration is altered by a change of stomatal conductance, the resulting change of the advective flux of soil water to roots may modify mass flow-dependent nutrient uptake [13, 14]. Such a change may alter the nitrogen nutrition status of vegetation . Also, reduced transpiration has been identified as one of the controls limiting N acquisition in elevated CO2 studies in the field . Understanding relationships between stomatal conductance, photosynthesis, transpiration and N acquisition are fundamental for explaining why for some temperate humid grasslands the CO2 fertilization responses in aboveground biomass production have been small or absent [4, 21]. That such interactions may have already operated in the last century is indicated by an investigation of herbage yields in 1891–1992 on several unlimed plots of the Park Grass Continuous Hay Experiment (Rothamsted, U.K.), which found no significant trend , although water use efficiency has increased . It has been hypothesized from elevated CO2 studies in the field that the limited fertilization effect of elevated ca in temperate grassland communities is connected with stronger reduction of stomatal conductance in the grasses, and greater photosynthetic acclimation, in comparison with forbs and legumes, or by nutrient limitation, especially for N [4, 7, 21]. Here, we examine evidence from the Park Grass Experiment to see in how far these mechanisms have already operated to shape grassland communities’ responses to the increase of ca and associated climate change in the last century.
The Park Grass Experiment  (hereafter referred to as Park Grass) is the oldest permanent grassland experiment in the world. As a field resource, it enables not only exploration of the relationships between grassland community composition, nutrient status, and herbage yields, but also canopy stomatal conductance, water use efficiency, and N acquisition by ways of elemental and isotopic analysis of samples stored in the Rothamsted Sample Archive (the “Methods” section). In this analysis, we also expand on the temporal scope of earlier studies on yield trends  and intrinsic water use efficiency , by investigating changes during the last 100 years by using two periods in the early twentieth (1917–1931) and twenty-first century (2004–2018), each of 15 years. Over that time span (1917 to 2018), atmospheric CO2 concentration (ca) increased by c. 30% or 100 μmol mol−1 [25, 26].
Results and discussion
Between periods, daily mean temperature increased by c. 1.5 °C at Rothamsted, but atmospheric vapor pressure deficit and rainfall during the main spring growth period (1 April to 30 June) did not change significantly (Additional file 1: Fig. S1).
We analyzed hay samples (the “Methods” section) from selected plots that had received different annual fertilizer treatments (n = 12) for over a century, starting in 1856. Treatments included limed and unlimed “control” plots given no fertilizer and plots receiving nitrate- or ammonia-nitrogen (N) at different rates amended with or without lime plus phosphorus and potassium (PK) fertilizer. These treatments have caused a great diversity of plant species richness [24, 27] (2–43 species) in associated grass- (80–100% grass biomass, n = 6 treatments) and dicot-rich communities (32–52% dicots, with 2–25% legumes, n = 6), and variable total herbage yields [22, 27] (100–560 g dry biomass m−2 year−1) in the 1st annual yield cut taken around mid-June each year (Table 1). High yields (yields > 80% of the highest yielding treatment) were obtained over the full range of N fertilizer application (0–14.4 g N m−2 year−1) in the limed treatments, providing that plots also received P and K fertilizers (Fig. 1).
ci/ca was determined from 13C discrimination (the “Methods” section), in accounting for effects of mesophyll conductance and photorespiration , factors that were not previously considered . Equations 1 and 2 show that water use efficiency depends principally on photosynthesis and stomatal conductance, if atmospheric water demand is constant, as was approximately the case at Park Grass (Additional file 1: Fig. S1). Additionally, water use efficiency increases in proportion to ca if 1 − ci/ca, the relative CO2 concentration gradient between the atmosphere and the leaf intercellular space, does not vary. Our results demonstrated that ci/ca did not change significantly between the start and the end of the century in the dicot-rich communities, but decreased by 0.01 in grass-rich communities (P < 0.01). Besides, ci/ca was slightly higher (and varied less within periods) in dicot-rich relative to grass-rich treatments, both at the beginning (+ 0.01, P < 0.001) and at the end of the century (+ 0.02, P < 0.001) (Fig. 2). The modest observed variation of ci/ca across the century agrees with other 13C-based observations over geological timescales with greatly varying ca  and with plant responses from seasonal- to century-scale changes of climatic, edaphic, and nutritional stresses in forests and grasslands, including at Park Grass [5, 23, 32]. Modest variation of ci/ca is an expected result of coordinated reciprocal adjustments of photosynthesis and stomatal conductance that serve to optimize photosynthetic C gain relative to water loss by transpiration [5, 12, 18, 33]. On average, intrinsic water use efficiency at Park Grass increased by + 9.6 μmol mol−1 or c. 31% compared to the beginning of the century (Fig. 3a, P < 0.001), mainly due to the increase of ca. However, the increase of intrinsic water use efficiency was c. 33% greater in the grass-rich swards than dicot-rich plant communities, due to the small decrease of ci/ca.
Yield trends over the century diverged for the dicot- and grass-rich treatments (Fig. 3b, P < 0.001). While the yields of the grass-rich treatments decreased by 16% on average, yield trends of the dicot-rich plots and the average of all plots were not significant, aligning with the previous observation that annual yields of the Park Grass plots have not increased systematically in the twentieth century . Changes of yield (Y) are linked to changes of growing season-integrated water use efficiency via canopy photosynthesis as Y = A (1 –ϕ) (1 – r) , with ϕ the proportion of carbon respired and r that allocated to roots and non-harvested or senesced shoot biomass. Growing season-integrated canopy photosynthesis of the different treatments may have changed over the century in response to rising ca, warming (and associated earlier spring growth) and precipitation patterns and weather dynamics in intricate but unknown ways. Also, there is some uncertainty about long-term and treatment effects on ϕ and r (but see the “Methods” section). Interestingly, however, yield trends were closely negatively related to intrinsic water use efficiency (R2 = 0.53; P < 0.01). That negative relationship points to stomatal conductance (integrated over the growing season and canopy) as the likely primary control of treatment-dependent yield trends in the last century.
Indeed, treatment effects on yield trends also correlated negatively with trends of δ18O of cellulose (δ18Ocellulose, the “Methods” section) (Fig. 4b, R2 = 0.52, P < 0.01). Increases of δ18Ocellulose in plant biomass are indicative of decreases of stomatal conductance, under otherwise equal environmental conditions [34, 35], including rainfall, atmospheric humidity, soil conditions affecting plant-water relations, and isotopic inputs of rain (δ18Orain) and atmospheric moisture (δ18Ovapour) (the “Methods” section). Climate warming should increase δ18Orain  and has caused an increase of vapor pressure deficit in Europe, which contributes to reduce stomatal conductance . Yet, vapor pressure deficit and rainfall during spring growth have not changed significantly at Park Grass. Secondly, measurements at a 55 km-distant station (Wallingford) and predictions from global circulation models provide no indication for significant trends of annual δ18Orain (Additional file 1: Fig. S2). And thirdly, the long-term trends of δ18Ocellulose were independent of interannual variations of meteorological variables (such as vapor pressure deficit), which were very similar at the beginning of the 20th and the twenty-first century (Fig. 5). Most importantly, the century-scale divergence of δ18Ocellulose between treatments was not attributable to eventual changes of δ18Orain or environmental conditions, as all treatments experienced the same site and weather conditions.
On average of all treatments, δ18Ocellulose increased over the century (Fig. 3c, average + 0.7‰, P < 0.001), consistent with a general long-term decrease of stomatal conductance. The increase of δ18Ocellulose was much stronger in grass-rich (+ 1.1‰; P < 0.001) than dicot-rich (+ 0.3‰; P = 0.06) plots. The trends of δ18Ocellulose were closely proportional (R2 = 0.70; P < 0.001) to trends in the ratio of hay yield to the CO2 concentration gradient between the atmosphere and the leaf internal gas space, ca – ci (Fig. 6a), another proxy of stomatal conductance integrated over the growing season and canopy (the “Methods” section). The proportionality indicated that a long-term change in δ18Ocellulose only occurred when stomatal conductance changed. Additionally, the sensitivity of δ18Ocellulose changes to changes of stomatal conductance estimated for Park Grass agreed with sensitivities observed by others (the “Methods” section) in leaf- or canopy-scale stomatal conductance studies with a range of plant species or genotypes in different environmental conditions (Fig. 6b), instilling further confidence in δ18Ocellulose change as a well-grounded measure of long-term canopy-scale stomatal conductance variation at the Park Grass site in the past century. It is also notable in that context, that interannual variation of δ18Ocellulose correlated negatively with ci/ca, both at the beginning and end of the century (P < 0.001), with greater ranges of variation in grass-rich than dicot-rich plots (Fig. 2b). This relationship also points to stomatal conductance-variation as an important cause for interannual variation of intrinsic WUE at Park Grass, particularly in the grass-rich plots.
Negative yield trends in grass-rich communities over the century were connected with decreased net rates of N acquisition (Fig. 3d, − 1.9 g N m−2 year−1, P < 0.001). Similar decreases of N acquisition were generally not observed in dicot-rich treatments (− 0.3 g m−2 year−1, P = 0.4), that received little or no N fertilizer. These observations match observations in elevated CO2 studies, which found consistently decreased rates of N acquisition in grassland, forests and crops under conditions where elevated CO2 failed to enhance yields . The Park Grass data indicate that such ca- and climate change-induced reductions of N acquisition may have occurred mainly in grass-rich grassland communities in the last century. At Park Grass, these decreases in N acquisition were closely correlated with increases of δ18Ocellulose (Fig. 4c, R2 = 0.56; P < 0.01), in agreement with stomatal conductance-mediated reductions of transpiration and associated mass flow-dependent N acquisition [3, 13, 14]. Absence of significant reductions of N acquisition in the non-N-fertilized (dicot-rich) communities agrees with much smaller (if any) reductions of stomatal conductance, but may also be influenced by the fact that N acquisition in N-poor soils is more closely controlled by diffusion , and atmospheric deposition of N and biological N fixation account for a larger proportion of total N acquisition in these plots on average . But, regardless of the mechanism underlying the reduced N acquisition  in grass-rich plots, any impairment of N status would reduce photosynthetic capacity and feedback on stomatal conductance to adjust ci/ca for optimization of C gain per unit water loss by transpiration [6, 7, 10, 15, 16, 18].
N fertilizer addition was the strongest experimental determinant of grass dominance at Park Grass (P < 0.001). Meanwhile, PK fertilization had a much greater effect on yields, although all plots remained N limited or N and P co-limited (Table 1). N limitation [2, 20, 28] is a well-known factor limiting the CO2 fertilizer effect in grassland [4, 21]. During the century, the reduction of N acquisition aggravated the N limitation of the grass-rich communities, as evidenced by a decreasing nitrogen nutrition index  (− 0.08, P < 0.001). Conversely, the nitrogen nutrition index decreased much less in the dicot-rich communities (− 0.04, P = 0.01). That difference must have contributed to the convergence of yields of grass- and dicot-rich plots over the century. For example, among the high-yielding, limed and PK-fertilized plots, the dicot-rich, non-N-fertilized plot yielded 22% less than its grass-rich, N-fertilized counterparts at the beginning of the twentieth century, but only 6% less at the end of the twenty-first century. Thus, the yield enhancement generated by N fertilizer addition decreased by 78% over the century. At the same time, the apparent efficiency of N fertilizer uptake (net N acquisition per unit of added N fertilizer) decreased by 16% in the grass-rich (P < 0.001) but not in the dicot-rich communities (P = 0.2). Although the cutting frequency (two cuts per year) at Park Grass does not represent a very intensive utilization practice, the decreasing N fertilizer uptake efficiency of grass-rich communities does raise questions concerning the fate of N in the ecosystem and the future sustainability of high inputs of N fertilizer in such grasslands . Much of the Western European grassland receives very high N inputs .
The present work provides evidence for an interaction between plant functional groups and the ca- and related climate change-response of temperate-humid permanent grassland in the last century, managed under a hay-cutting regime with a range of fertilizer inputs. While part of that evidence is correlative, its interpretation is underpinned by empirical knowledge and mechanistic understanding obtained in more controlled experimental conditions, although these were mainly concerned with the effects of future elevated CO2 and changed climate [3, 6, 7, 21]. The observed interactions between plant functional group composition and climate change responses in the last century appeared to be associated primarily with the much greater ca-sensitivity of grasses in terms of their stomatal conductance compared with that of forbs and legumes, resulting in effects on grassland vegetation water use efficiency , N acquisition and aboveground biomass production, central features of the role of vegetation in biogeochemical cycles and feedbacks to the climate system. Our work also shows that key predictions drawn from free-air CO2 enrichment experiments concerning ca-saturation of the CO2 fertilizer effect on aboveground biomass production [4, 6, 21], N acquisition  and the lesser competitive ability of grasses relative to forbs and legumes [6, 21] have already taken hold in some temperate grasslands in the last century or before.
The Park Grass Experiment
The Park Grass Experiment (hereafter referred to as Park Grass), located at the Rothamsted Research Station in Hertfordshire, U.K., approx. 40 km north of London (0°22′ West, 51°48′ North), is the oldest grassland experiment in the world. It was started in 1856 on c. 2.8 ha of old grassland, and its original purpose was to investigate the effect of fertilizers and organic manures on hay yields of permanent grassland . The experiment comprises 20 main plots with different fertilizer inputs; most plots have been divided into subplots “a”–“d” and receive different lime inputs. Since the beginning of the experiment, herbage has been cut in mid-June and made into hay (first cut). Herbage regrowth of the sward was grazed by sheep during the first 15 years, but in each year since 1875, a second cut has been taken and removed from the plot. Originally, total herbage yields were determined in situ by weighing the dried material (hay) for whole plots. Since 1960, one or two sample strips, depending on the plot size, were cut on each plot using a forage harvester. The fresh material was weighed and yields were calculated after drying. Since 1856, representative hay samples from each harvest have been stored in the Rothamsted Sample Archive .
For this study, we used dried hay or herbage samples from the first cut (spring growth) of 13 plots (12 different treatments), with contrasting fertilizer and lime inputs (Table 1). The treatments included four nitrogen levels (0, 4.8, 9.6, and 14.4 g N m−2 year−1, with nitrogen added either as ammonium sulfate or sodium nitrate), lime application (chalk applied as necessary to maintain soil at pH 7, 6, and 5 on sub-plots “a,” “b,” and “c”), one sodium/potassium level (3.5 g P m−2 year−1 as triple superphosphate and 22.5 g K m−2 year−1 as potassium sulfate; P was applied at 1.7 g m−2 year−1 since 2017), and the control plots (no fertilizer application). The plots that received P and K also received 1.5 g Na m−2 year−1 as sodium sulfate and 1 g Mg m−2 year−1 as magnesium sulfate. N was applied in spring and P and K in autumn.
Functional group composition
Plant functional group (PFG) composition was strongly altered by fertilization and lime inputs on Park Grass. The original vegetation was classified post-hoc as dicotyledon-rich Cynosurus cristatus-Centaurea nigra grassland , but it changed rapidly following the introduction of the different treatments, creating very contrasting community compositions, which have reached a dynamic equilibrium since 1910 . We consider the functional group composition (%-grasses, %-forbs and %-legumes in standing dry biomass) of each treatment to have remained relatively constant during the past 100 years. The average functional group composition was calculated from individual datasets sampled between 1915 and 1976, 1991 and 2000, and 2010 and 2012 (e-RA, http://www.era.rothamsted.ac.uk/). Given the high but variable proportion of grasses (48–100% of standing dry biomass, depending on treatment), we used this characteristic to represent the functional composition heterogeneity among treatments. Treatments with a grass proportion > 80% were termed grass-rich, and the other treatments (grass proportion 50–68%) dicot-rich (Table 1).
Sample preparation and analysis
Representative subsamples (2–3 g) of plant material from the first cut were taken from the Rothamsted sample archive. Two 15-year periods were selected, corresponding to the beginning (1917–1931, period 1) and the end (2004–2018, period 2) of the last 100 years. The period length was selected to account for the interannual variability in δ18O of α-cellulose (δ18Ocellulose), which is associated with interannual variation of meteorological conditions (see the “Stomatal conductance” section and Fig. 5). The subsamples were used to determine the δ13C of bulk plant material (δ13Cp) and δ18Ocellulose. Carbon and oxygen isotope composition were expressed as:
with Rsample the 13C/12C (or 18O/16O) ratio of the sample and Rstandard that in the international standard (Vienna Pee Dee Belemnite, V-PDB for carbon or Vienna Standard Mean Ocean Water, V-SMOW for oxygen isotope ratio).
α-cellulose was extracted from 50 mg of dry sample material by following the Brendel et al.  protocol as modified by Gaudinski et al. . The extracted cellulose was re-dried at 40 °C for 24 h; 0.7 ± 0.05 mg aliquots were weighed into silver cups (size 3.3 × 5 mm, IVA Analysentechnik e.K., Meerbusch, Germany) and stored in an exsiccator vessel containing Silica Gel Orange (ThoMar OHG, Lütau, Germany) prior to analysis. Samples were pyrolyzed at 1400 °C in a pyrolysis oven (HTO, HEKAtech, Wegberg, Germany), equipped with a helium-flushed zero blank auto-sampler (Costech Analytical technologies, Valencia, CA, USA) and connected via an interface (ConFlo III, Finnigan MAT, Bremen, Germany) with a continuous-flow isotope ratio mass spectrometer (Delta Plus, Finnigan MAT). Solid internal laboratory standards (SILS, cotton powder) were run each time after the measurement of four samples. Both cellulose samples and SILS were measured against a laboratory working standard carbon monoxide gas, previously calibrated against a secondary isotope standard (IAEA-601, accuracy of calibration ± 0.25‰ SD). The precision for the laboratory standard was < 0.3‰ (SD for repeated measurements).
For measurement of δ13Cp, the plant material was dried at 40 °C for 48 h, ball milled to a fine powder, re-dried during 24 h at 60 °C, and aliquots of 0.7 ± 0.05 weighed into tin cups (size 3.3 × 5 mm, LüdiSwiss, Flawil, Switzerland). Samples were combusted in an elemental analyzer (NA 1110; Carlo Erba, Milan, Italy) interfaced (Conflo III; Finnigan MAT, Bremen, Germany) with an isotope ratio mass spectrometer (Delta Plus; Finnigan MAT) and measured against a laboratory working standard CO2 gas, previously calibrated against a secondary isotope standard (IAEA-CH6 for δ13C, accuracy of calibration 0.06% SD). As a control, SILS which were previously calibrated against the secondary standard and had a similar C/N ratio as the samples (a fine wheat flour) were run after every tenth sample. The long-term precision for the SILS was < 0.2‰. C and N elemental concentration (%C and %N in dry biomass) were also measured in the same sequence. Additionally, P elemental concentration (%P in dry biomass) was determined on 20–25 mg of dry sample material for a 10-year period (2000–2009) using a modified phosphovanado-molybdate colorimetric method following acid digestion .
Meteorological and yield data
Meteorology and yield datasets were obtained from the electronic Rothamsted Archive (e-RA, http://www.era.rothamsted.ac.uk/). Rainfall data have been collected continuously since 1853 at the Rothamsted Weather Station, temperature since 1873 and additional meteorological variables have been added gradually. We calculated a yearly spring value for the last 100 years (1917–2018) by averaging all daily values between 1 April and 30 June for the following meteorological variables: maximum air temperature (Tmax, °C), average air temperature (Tmean, °C), average relative humidity of air (RH, %), and average vapor pressure deficit (VPD, kPa). In the case of precipitation (rain, mm), we used the cumulative rain amount for the same period. Yearly spring values were used for the analysis of long-term climatic changes.
Total herbage yield data from the first cut between 1917 and 2018 (expressed in g dry biomass m−2 year−1) were used to determine long-term yield trends during spring growth and for the calculation of (1) the yearly nitrogen acquisition (N acquisition, g m−2 year−1), calculated as dry matter yield × N content of dry biomass, and (2) a yield-weighted canopy stomatal conductance (see the “Stomatal conductance” section), for each of the 12 selected treatments. Given the change in the harvest method since 1960 [22, 46], we calculated an offset-correction for yield estimates pre-1960 to permit estimation of yield changes over the century on the same methodological basis. The offset-correction was based on the relationship between atmospheric CO2 concentration (as a long-term proxy for climate change) and measured yields from 1917 to 2018. For this, we fitted the following multiple linear regression to the data of every treatment:
with Y denoting measured yield, β and γ the coefficients to be determined, period a dummy variable coded with 0 for yields before 1960 and 1 for yields after 1960 and ε the random error. The treatment-specific estimated γ parameter was used as the offset-correction for spring yield determinations before 1960 (Additional file 1: Fig. S3 and Table S3).
Intrinsic water use efficiency
Intrinsic water use efficiency (Wi) is a physiological efficiency  that represents the ratio of net photosynthesis rate (A) and stomatal conductance to water vapor (gw):
That relationship is controlled by the CO2 concentration gradient between the atmosphere (ca) and the leaf internal gas space (ci), ca - ci, and the ratio of the diffusivities of water vapor and CO2 (1.6) [19, 47]. The CO2 concentration gradient can also be expressed as the product of ca and the relative CO2 concentration gradient between the atmosphere and the leaf internal gas space (1 – ci/ca) (Eq. 5). Thus, solving the equation required two parameters, ca, known from measurements of free air and air bubbles separated from ice cores , and ci/ca, that we estimated according to Ma et al.  (the bracketed term in Eq. 5). In that term, a is the 13C discrimination during diffusion of CO2 in air through the stomatal pore (4.4‰), am (1.8‰) the fractionation associated with CO2 dissolution and diffusion in the mesophyll, and b (29‰) and f (11‰) the fractionations due to carboxylation and photorespiration, respectively; Г* is the CO2 compensation point in the absence of mitochondrial respiration calculated following Brooks and Farquhar , and gs/gm the ratio of stomatal and mesophyll conductance . Relying on theory developed by Farquhar et al.  and Farquhar and Richards , this model accounts for the effects of mesophyll conductance and photorespiration on Δ13C-based estimations of intrinsic water use efficiency. Traditionally, estimations of Wi from Δ13C have used a more abbreviated version of the Farquhar model of photosynthetic carbon isotope discrimination (Δ13C) in C3 plants. That version neglected the photorespiration term and was based on the simplifying assumption that mesophyll conductance is infinite. However, Wi estimated from Δ13C using the abbreviated model systematically overestimates gas exchange-based estimates of Wi, an error that can be corrected by using a constant gs/gm ratio (0.79) based on measurements on a wide range of plant species from different functional groups (including grasses and herbaceous legumes), in moist and dry conditions . Also, global scale effects of mesophyll conductance and photorespiration of C3 vegetation on 13C discrimination of the terrestrial biosphere have been evident for the last four decades of atmospheric CO2 increase and concurrent change of the 13CO2/12CO2 ratio . Estimates of intrinsic water use efficiency at Park Grass calculated following Ma et al.  were closely correlated (R2 = 0.97) with those presented by Köhler et al.  using the abbreviated Wi model for the same treatments and time spans.
Δ13C was obtained from the measured δ13Cp of samples and the δ13C of atmospheric CO2 (δ13Ca), as
Estimates of ca and δ13Ca were obtained as in Wittmer et al.  and Köhler et al. . The time series generated by Köhler et al.  with ca and δ13Ca yearly average values during spring (April–June, 1917–2009) was updated until 2018. Monthly values between 2010 and 2018 from the NOAA ESRL atmospheric stations Mauna Loa, Mace Head, and Storhofdi Vestmannaeyjar [26, 53] were used for ca and δ13Ca estimation. Calculations of Δ13C and Wi were performed for each treatment, individually.
As intrinsic water use efficiency was estimated from a representative sample of a hay harvest, it represents a growing-season assimilation- and allocation-weighted measure.
Extrapolating from elevated CO2 studies , we hypothesized that the atmospheric CO2 increase had already caused partial stomatal closure in grassland vegetation during the last century. This would also contribute to explain why intrinsic water use efficiency has increased at Park Grass  although herbage dry matter yields have not, generally, shown a similar increase. We used two independent approaches to estimate changes of stomatal conductance at Park Grass: one based on changes of yields, atmospheric CO2 concentration and 13C discrimination, and the other on the concurrent changes of the oxygen isotope composition of cellulose, as we explain below.
Following reasoning in Farquhar and Richards , canopy-scale, growing season-integrated stomatal conductance to CO2 (gs, in mol m−2 s−1) can be estimated as:
with Y denoting yield (in moles of C in harvested biomass per ground area and per year), ϕ the proportion of carbon respired and r that allocated to roots and non-harvested shoot biomass. As measurements of A were unavailable, we used the yield data as a proxy and assumed that ϕ (0.35) and r (0.4) were constants across treatments and during the century. From estimates of stomatal conductance for every year, we then calculated the percent change of stomatal conductance (%gs-change) between the beginning (period 1, subscript 1) and end (period 2, subscript 2) of the century for every treatment as:
Note that this procedure eliminated the constants for ϕ and r, so that the estimated %gs-changes were independent of ϕ and r, except if they had changed during the century.
For the 18O-based inference, theory [34, 54] and observations [35, 55,56,57] demonstrate that changes of stomatal conductance (gs) cause inverse changes of 18O-enrichment of cellulose above source water (Δ18Ocellulose, with Δ18Ocellulose ≈ δ18Ocellulose − δ18Osource). If the δ18O of source water (δ18Osource) – the water taken up by the root system – is invariant, then an increase of Δ18Ocellulose implies a parallel increase of δ18Ocellulose. Given that δ18Osource is determined by the δ18O of precipitation (δ18Orain) and plots at Park Grass receive the same precipitation, we presumed that any divergence in the δ18Ocellulose changes was independent of δ18Osource and hence attributable to changes of Δ18Ocellulose. Indeed, no long-term changes were detected for δ18Orain at Park Grass, estimated using yearly averages from the nearest monitoring station (Wallingford GNIP, c. 55 km west-southwest of Park Grass, data since 1982)  and outputs from the ECHAM5 global circulation model  for the location of the Park Grass Experiment (pixel resolution 1.125° × 1.12°, data since 1958) (Additional File 1: Fig. S2).
The effect of stomatal conductance on Δ18Ocellulose is primarily determined by three factors: (1) all oxygen in cellulose originates from water [60, 61], (2) a high proportion of that water is leaf water, which is evaporatively 18O-enriched during daytime and imprints its 18O-signal on photosynthetic products used in cellulose synthesis , and (3) evaporative 18O-enrichment of leaf water decreases with increasing stomatal conductance . The δ18Ocellulose of a representative hay sample is a canopy-scale, growing season-integrated signal reflecting assimilation over the total time span that contributed substrate for cellulose synthesis in the tissues that comprise the sample . However, the relation between stomatal conductance and Δ18Ocellulose is influenced by multiple morpho-physiological vegetation characteristics and environmental variables , so that stomatal conductance cannot be simply and quantitatively estimated from Δ18Ocellulose or δ18Ocellulose at present . As no direct measurements of stomatal conductance were available at Park Grass during the last century, we compared δ18Ocellulose changes during the century with concurrent changes of stomatal conductance as estimated by Eqs. (7) and (8), to obtain an estimate of the stomatal conductance sensitivity of a unit change of δ18Ocellulose. Further, we compared the δ18Ocellulose change versus %gs-change relation with published data [35, 55,56,57] and one unpublished data set. The published data covered a wide range of plant species or genotypes in different environmental conditions in controlled environments and in the field. The data from the graphic displays of the different studies were digitized using ImageJ . In this analysis, we also included (1) the relationship between Δ18Ocellulose and canopy-scale stomatal conductance predicted by a process-based 18O-enabled soil-plant-atmosphere model for a multi-seasonal data set of grassland  and (2) an unpublished data set from monocultures of Lolium perenne grown at CO2 concentrations of 200, 400, or 800 μmol mol−1 in controlled environments with air temperature controlled at 20 °C:16 °C and relative humidity at 50%:75% during the 16:8 h day to night cycle, and a photosynthetic photon flux density of 800 μmol m−2 s−1 during the photoperiod .
Individual sensitivities for each study were calculated as the slopes of ordinary least-squares linear regressions. We then calculated a mean sensitivity ± CI95% (expressed as the decrease of stomatal conductance per unit increase in δ18Ocellulose) based on the sensitivities in the individual studies (n = 6, Additional file 1: Table S2). Additionally, we used an analogous procedure to calculate the stomatal conductance-sensitivity in terms of %gs-decrease (Eq. 8) per unit increase in δ18Ocellulose. In that, we took the mean stomatal conductance value of each study as the reference for calculating the %-decreases of stomatal conductance using the stomatal conductance versus δ18Ocellulose relationship of that study. On average of all the studies, stomatal conductance decreased by 39% per 1‰ increase of δ18Ocellulose.
We tested the long-term response of intrinsic water use efficiency, yield, δ18Ocellulose, and N acquisition by comparing the average values of 15-year periods at the beginning (1917–1931, period 1) and the end (2004–2018, period 2) of the last 100 years. The long-term change was calculated as the difference between the averages from period 1 and 2, for (1) every treatment individually (n = 26–30), (2) grouping the data into dicot-rich (n = 164–173) and grass-rich (n = 161–174) treatments, and (3) combining all data together (n = 330–348) (see Additional file 1: Table S1 for detailed information about the sample size). t tests were used for testing the differences between periods. Additional t tests were performed to determine whether the long-term response of the analyzed variables differed between treatments with different plant functional group composition (dicot-rich vs. grass-rich communities). For this, data corresponding to period 2 were normalized by subtracting the average value from period 1, for every individual treatment (thus obtaining the net long-term change). Additionally, the relationship between the long-term changes of different variables (e.g., change in intrinsic water use efficiency vs. change in δ18Ocellulose) was tested using ordinary least-squares linear regressions. Regression analysis was also used to test the effect of multiple parameters on yield. All statistical analyses were conducted in R v.4.0.2 . The R package ggplot2  was used for data plotting.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Craine JM, Elmore AJ, Wang L, Aranibar J, Bauters M, Boeckx P, et al. Isotopic evidence for oligotrophication of terrestrial ecosystems. Nat Ecol Evol. 2018;2:1735–44.
Du E, Terrer C, Pellegrini AFA, Ahlström A, van Lissa CJ, Zhao X, et al. Global patterns of terrestrial nitrogen and phosphorus limitation. Nat Geosci. 2020;13:221–6.
Feng Z, Rütting T, Pleijel H, Wallin G, Reich PB, Kammann CI, et al. Constraints to nitrogen acquisition of terrestrial plants under elevated CO2. Glob Chang Biol. 2015;21:3152–68.
Körner C. Plant CO2 responses: an issue of definition, time and resource supply. New Phytol. 2006;172:393–411.
Lavergne A, Voelker S, Csank A, Graven H, de Boer HJ, Daux V, et al. Historical changes in the stomatal limitation of photosynthesis: empirical support for an optimality principle. New Phytol. 2020;225:2484–97.
Ainsworth EA, Long SP. 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. 2005;165:351–72.
Ainsworth EA, Rogers A. The response of photosynthesis and stomatal conductance to rising [CO2]: mechanisms and environmental interactions. Plant Cell Environ. 2007;30:258–70.
Berry JA, Beerling DJ, Franks PJ. Stomata: key players in the earth system, past and present. Curr Opin Plant Biol. 2010;13:232–9.
Hetherington AM, Woodward FI. The role of stomata in sensing and driving environmental change. Nature. 2003;424:901–8.
Medlyn BE, Duursma RA, Eamus D, Ellsworth DS, Prentice IC, Barton CVM, et al. Reconciling the optimal and empirical approaches to modelling stomatal conductance. Glob Chang Biol. 2011;17:2134–44.
Franks PJ, Adams MA, Amthor JS, Barbour MM, Berry JA, Ellsworth DS, et al. Sensitivity of plants to changing atmospheric CO2 concentration: from the geological past to the next century. New Phytol. 2013;197:1077–94.
Cernusak LA, Ubierna N, Winter K, Holtum JAM, Marshall JD, Farquhar GD. Environmental and physiological determinants of carbon isotope discrimination in terrestrial plants. New Phytol. 2013;200:950–65.
McGrath JM, Lobell DB. Reduction of transpiration and altered nutrient allocation contribute to nutrient decline of crops grown in elevated CO2 concentrations. Plant Cell Environ. 2013;36:697–705.
Oyewole OA, Inselsbacher E, Näsholm T. Direct estimation of mass flow and diffusion of nitrogen compounds in solution and soil. New Phytol. 2014;201:1056–64.
Kattge J, Knorr W, Raddatz T, Wirth C. Quantifying photosynthetic capacity and its relationship to leaf nitrogen content for global-scale terrestrial biosphere models. Glob Chang Biol. 2009;15:976–91.
Smith NG, Keenan TF. Mechanisms underlying leaf photosynthetic acclimation to warming and elevated CO2 as inferred from least-cost optimality theory. Glob Chang Biol. 2020;26:5202–16.
Wright IJ, Reich PB, Westoby M. Least-cost input mixtures of water and nitrogen for photosynthesis. Am Nat. 2003;161:98–111.
Buckley TN. Modeling stomatal conductance. Plant Physiol. 2017;174:572–82.
Farquhar GD, Richards RA. Isotopic composition of plant carbon correlates with water-use efficiency in wheat genotypes. Aust J Plant Physiol. 1984;11:539–52.
Kunrath TR, Lemaire G, Sadras VO, Gastal F. Water use efficiency in perennial forage species: interactions between nitrogen nutrition and water deficit. Field Crops Res. 2018;222:1–11.
Lüscher A, Daepp M, Blum H, Hartwig UA, Nösberger J. Fertile temperate grassland under elevated atmospheric CO2—role of feed-back mechanisms and availability of growth resources. Eur J Agron. 2004;21:379–98.
Jenkinson DS, Potts JM, Perry JN, Barnett V, Coleman K, Johnston AE. Trends in herbage yields over the last century on the Rothamsted Long-term continuous Hay experiment. J Agric Sci. 1994;122:365–74.
Köhler IH, Macdonald AJ, Schnyder H. Last-century increases in intrinsic water-use efficiency of grassland communities have occurred over a wide range of vegetation composition, nutrient inputs, and soil pH. Plant Physiol. 2016;170:881–90.
Storkey J, Macdonald AJ, Poulton PR, Scott T, Köhler IH, Schnyder H, et al. Grassland biodiversity bounces back from long-term nitrogen addition. Nature. 2015;528:401–4.
Friedli H, Lötscher H, Oeschger H, Siegenthaler U, Stauffer B. Ice core record of the 13C/12C ratio of atmospheric CO2 in the past two centuries. Nature. 1986;324:237–8.
Dlugokencky EJ, Mund JW, Crotwell AM, Crotwell MJ, Thoning KW. Atmospheric carbon dioxide dry air mole fractions from the NOAA ESRL carbon cycle cooperative global air sampling network, 1968-2018. NOAA ESRL. 2019. https://doi.org/10.15138/wkgj-f215.
Crawley MJ, Johnston AE, Silvertown J, Dodd M, de Mazancourt C, Heard MS, et al. Determinants of species richness in the Park Grass Experiment. Am Nat. 2005;165:179–92.
Lemaire G, Jeuffroy M-H, Gastal F. Diagnosis tool for plant and crop N status in vegetative stage: theory and practices for crop N management. Eur J Agron. 2008;28:614–24.
Liebisch F, Bünemann EK, Huguenin-Elie O, Jeangros B, Frossard E, Oberson A. Plant phosphorus nutrition indicators evaluated in agricultural grasslands managed at different intensities. Eur J Agron ronomy. 2013;44:67–77.
Güsewell SN. P ratios in terrestrial plants: variation and functional significance. New Phytol. 2004;164:243-66.
Ma WT, Tcherkez G, Wang XM, Schäufele R, Schnyder H, Yang Y, et al. Accounting for mesophyll conductance substantially improves 13C-based estimates of intrinsic water-use efficiency. New Phytol. 2021;229:1326–38.
Cornwell WK, Wright IJ, Turner J, Maire V, Barbour MM, Cernusak LA, et al. Climate and soils together regulate photosynthetic carbon isotope discrimination within C3 plants worldwide. Glob Ecol Biogeogr. 2018;27:1056–67.
Cowan IR, Farquhar GD. Stomatal function in relation to leaf metabolism and environment. In: Jennings DH, editor. Integration of activity in the higher plant. Cambridge: University Press; 1977. p. 471–505.
Farquhar GD, Barbour MM, Henry BK. Interpretation of oxygen isotope composition of leaf material. In: Griffiths H, editor. Stable isotopes: integration of biological, ecological and geochemical processes. Oxford: BIOS Scientific Publishers; 1998. p. 27–61.
Barbour MM, Fischer RA, Sayre KD, Farquhar GD. Oxygen isotope ratio of leaf and grain material correlates with stomatal conductance and grain yield in irrigated wheat. Funct Plant Biol. 2000;27:625–37.
Bowen GJ, Cai Z, Fiorella RP, Putman AL. Isotopes in the water cycle: regional- to global-scale patterns and applications. Annu Rev Earth Planet Sci. 2019;47:453–79.
Grossiord C, Buckley TN, Cernusak LA, Novick KA, Poulter B, Siegwolf RTW, et al. Plant responses to rising vapor pressure deficit. New Phytol. 2020;226:1550–66.
Robertson GP, Vitousek PM. Nitrogen in agriculture: balancing the cost of an essential resource. Annu Rev Environ Resour. 2009;34:97–125.
Cheng L, Zhang L, Wang Y-P, Canadell JG, Chiew FHS, Beringer J, et al. Recent increases in terrestrial carbon uptake at little cost to the water cycle. Nat Commun. 2017;8:110.
Lawes JB, Gilbert JH. Report of experiments with different manures on permanent meadow land: part II: produce of constituents per acre. J R Agric Soc Engl. 1859;20:228–46.
Silvertown J, Poulton P, Johnston E, Edwards G, Heard M, Biss PM. The Park Grass Experiment 1856–2006: its contribution to ecology. J Ecol. 2006;94:801–14.
Dodd ME, Silvertown J, McConway K, Potts J, Crawley M. Application of the British national vegetation classification to the communities of the Park Grass Experiment through time. Folia Geobot. 1994;29:321–34.
Brendel O, Iannetta PPM, Stewart D. A rapid and simple method to isolate pure alpha-cellulose. Phytochem Anal. 2000;11:7–10.
Gaudinski JB, Dawson TE, Quideau S, Schuur EAG, Roden JS, Trumbore SE, et al. Comparative analysis of cellulose preparation techniques for use with 13C, 14C, and 18O isotopic measurements. Anal Chem. 2005;77:7212–24.
Hanson WC. The photometric determination of phosphorus in fertilizers using the phosphovanado-molybdate complex. J Sci Food Agric. 1950;1:172–3.
Bowley HE, Mathers AW, Young SD, Macdonald AJ, Ander EL, Watts MJ, et al. Historical trends in iodine and selenium in soil and herbage at the Park Grass Experiment, Rothamsted research. UK Soil Use Manag. 2017;33:252–62.
Ehleringer JR, Hall AE, Farquhar GD. Stable isotopes and plant carbon-water relations. San Diego: Academic Press; 1993.
Brooks A, Farquhar GD. Effect of temperature on the CO2/O2 specificity of ribulose-1,5-bisphosphate carboxylase/oxygenase and the rate of respiration in the light. Planta. 1985;165:397–406.
Farquhar GD, Oleary MH, Berry JA. On the relationship between carbon isotope discrimination and the intercellular carbon dioxide concentration in leaves. Aust J Plant Physiol. 1982;9:121–37.
Keeling RF, Graven HD, Welp LR, Resplandy L, Bi J, Piper SC, et al. Atmospheric evidence for a global secular increase in carbon isotopic discrimination of land photosynthesis. Proc Natl Acad Sci. 2017;114:10361–6.
Wittmer MHOM, Auerswald K, Tungalag R, Bai YF, Schäufele R, Bai CH, et al. Carbon isotope discrimination of C3 vegetation in Central Asian grassland as related to long-term and short-term precipitation patterns. Biogeosciences. 2008;5:913–24.
Köhler IH, Poulton PR, Auerswald K, Schnyder H. Intrinsic water-use efficiency of temperate seminatural grassland has increased since 1857: an analysis of carbon isotope discrimination of herbage from the Park Grass Experiment. Glob Chang Biol. 2010;16:1531–41.
White J, Vaughn B, Michel S. https://www.esrl.noaa.gov/gmd/dv/data/. Stable isotopic composition of atmospheric carbon dioxide (13C and 18O) carbon cycle cooperative global air sampling network, 1990–2014. ftp://aftp.cmdl.noaa.gov/data/trace_gases/co2c13/flask/ (2020). Accessed 10 Oct 2020.
Scheidegger Y, Saurer M, Bahn M, Siegwolf R. Linking stable oxygen and carbon isotopes with stomatal conductance and photosynthetic capacity: a conceptual model. Oecologia. 2000;125:350–7.
Sullivan PF, Welker JM. Variation in leaf physiology of Salix arctica within and across ecosystems in the High Arctic: test of a dual isotope (Δ13C and Δ18O) conceptual model. Oecologia. 2006;151:372.
Grams TEE, Kozovits AR, Häberle K-H, Matyssek R, Dawson TE. Combining δ13C and δ18O analyses to unravel competition, CO2 and O3 effects on the physiological performance of different-aged trees. Plant Cell Environ. 2007;30:1023–34.
Moreno-Gutiérrez C, Dawson TE, Nicolás E, Querejeta JI. Isotopes reveal contrasting water use strategies among coexisting plant species in a Mediterranean ecosystem. New Phytol. 2012;196:489–96.
Global Network of Isotopes in Precipitation. The GNIP Database. IAEA/WMO https://nucleus.iaea.org/wiser/index.aspx. Accessed 25 Oct 2020.
Werner M. ECHAM5-wiso simulation data- present-day, mid-Holocene, and last glacial maximum. PANGAEA; 2019. https://doi.org/10.1594/PANGAEA.902347.
DeNiro MJ, Epstein S. Relationship between the oxygen isotope ratios of terrestrial plant cellulose, carbon dioxide, and water. Science. 1979;204:51–3.
Liu HT, Gong XY, Schäufele R, Yang F, Hirl RT, Schmidt A, et al. Nitrogen fertilization and δ18O of CO2 have no effect on 18O-enrichment of leaf water and cellulose in Cleistogenes squarrosa (C4) – is VPD the sole control? Plant Cell Environ. 2016;39:2701–12.
Cernusak LA, Barbour MM, Arndt SK, Cheesman AW, English NB, Feild TS, et al. Stable isotopes in leaf water of terrestrial plants. Plant Cell Environ. 2016;39:1087–102.
Barbour MM. Stable oxygen isotope composition of plant tissue: a review. Funct Plant Biol. 2007;34:83–94.
Hirl RT, Ogée J, Ostler U, Schäufele R, Baca Cabrera JC, Zhu J, et al. Temperature-sensitive biochemical 18O-fractionation and humidity-dependent attenuation factor are needed to predict δ18O of cellulose from leaf water in a grassland ecosystem. New Phytol. 2021;229:3156–71.
Song X, Farquhar GD, Gessler A, Barbour MM. Turnover time of the non-structural carbohydrate pool influences δ18O of leaf cellulose. Plant Cell Environ. 2014;37:2500–7.
Roden JS, Farquhar GD. A controlled test of the dual-isotope approach for the interpretation of stable carbon and oxygen isotope ratio variation in tree rings. Tree Physiol. 2012;32:490–503.
Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods. 2012;9:671–5.
Baca Cabrera JC, Hirl RT, Zhu J, Schäufele R, Schnyder H. Atmospheric CO2 and VPD alter the diel oscillation of leaf elongation in perennial ryegrass: compensation of hydraulic limitation by stored-growth. New Phytol. 2020;227:1776–89.
R Core Team. R: a language and environment for statistical computing. Vienna, Austria. 2020. https://www.R-project.org/. Accessed 8 Jun 2020.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer Verlag; 2016.
We thank Iris H. Köhler for preliminary analyses and discussions, Jianjun Zhu for the support with sampling, Monika Michler for sample preparation, Anja Schmidt for the cellulose extraction, and the Rothamsted data curators for access to data from the Electronic Rothamsted Archive (e-RA).
Consent to participate
This work was supported by the Deutsche Forschungsgemeinschaft (DFG SCHN 557/9-1). The Rothamsted Long-term Experiments National Capability (LTE-NC) is supported by the UK BBSRC (Biotechnology and Biological Sciences Research Council, BBS/E/C/000 J0300) and the Lawes Agricultural Trust. Open Access funding enabled and organized by Projekt DEAL.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. Spring (April–June) meteorological data during the last 100 years. Figure S2. δ18Orain long-term trends. Figure S3. Yields at Park Grass 1917–2018. Table S1. Sample size for all treatments. Table S2. List of publications used for estimation of the average sensitivity of stomatal conductance (gs) to δ18Ocellulose. Table S3. Pre-1960 yield offset-correction.
About this article
Cite this article
Baca Cabrera, J.C., Hirl, R.T., Schäufele, R. et al. Stomatal conductance limited the CO2 response of grassland in the last century. BMC Biol 19, 50 (2021). https://doi.org/10.1186/s12915-021-00988-4