Rothamsted Repository Download

Background: The anthropogenic increase of atmospheric CO 2 concentration ( c a ) 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 c a -increase and climate change over the last century have already affected canopy stomatal conductance and its links with C and N processes in grassland. Results: Here, we assessed two independent proxies of (growing season-integrating canopy-scale) stomatal conductance changes over the last century: trends of δ 18 O in cellulose ( δ 18 O cellulose ) 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 CO 2 concentration gradient between the atmosphere and the leaf internal gas space ( c a – c i ). The two proxies correlated closely ( R 2 = 0.70), in agreement with the hypothesis. In addition, the sensitivity of δ 18 O cellulose 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 δ 18 O cellulose changes for historical reconstruction of stomatal conductance changes at Park Grass. Trends of δ 18 O cellulose 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. Conclusions: 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 CO 2 responsive forbs and N-fixing legumes. feed-backs to the climate system. Our work also shows that


Background
Atmospheric CO 2 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 CO 2 by photosynthesis and evaporative loss of water vapor in transpiration occur through the same path, the stomata [8]. These small pores in the leaf epidermis adjust their conductance in response to environmental conditions [9,10], such as atmospheric CO 2 concentration (c a ) [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 CO 2 uptake, they simultaneously expose the leaf internal moisture to a comparatively dry atmosphere. This leafto-air vapor pressure gradient drives transpiration [10]. 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 c a , becomes evident when A and E are expressed as the products of stomatal conductance (g s ) and the concentration gradients of the respective gases [12,19]: where c a and v a are the CO 2 and water vapor mole fractions in air, c i and v i that in the substomatal cavity, and 1.6 is the ratio of the diffusivities of CO 2 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 flowdependent nutrient uptake [13,14]. Such a change may alter the nitrogen nutrition status of vegetation [20]. Also, reduced transpiration has been identified as one of the controls limiting N acquisition in elevated CO 2 studies in the field [3]. Understanding relationships between stomatal conductance, photosynthesis, transpiration and N acquisition are fundamental for explaining why for some temperate humid grasslands the CO 2 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 [22], although water use efficiency has increased [23]. It has been hypothesized from elevated CO 2 studies in the field that the limited fertilization effect of elevated c a 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 c a and associated climate change in the last century.
The Park Grass Experiment [24] (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 [22] and intrinsic water use efficiency [23], by investigating changes during the last 100 years by using two periods in the early twentieth (1917)(1918)(1919)(1920)(1921)(1922)(1923)(1924)(1925)(1926)(1927)(1928)(1929)(1930)(1931) and twenty-first century (2004-2018), each of 15 years. Over that time span (1917 to 2018), atmospheric CO 2 concentration (c a ) 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 nitrateor 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).
We determined intrinsic water use efficiency (W i ), the physiological, i.e., non-climatic component of water use efficiency, that accounts separately for atmospheric water demand (v a -v i ) [12,19]: c i /c a was determined from 13 C discrimination (the "Methods" section), in accounting for effects of mesophyll conductance and photorespiration [31], factors that were not previously considered [23]. 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 c a if 1 − c i /c a , the relative CO 2 concentration gradient between the atmosphere and the leaf intercellular space, does not vary. Our results demonstrated that c i /c a 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, c i / c a 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 c i /c a across the century agrees with other 13 C-based observations over geological timescales with greatly varying c a [11] 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 c i /c a is an expected result of coordinated reciprocal adjustments of photosynthesis and stomatal conductance that serve Lime application is represented with the letters "L" (limed) or "U" (unlimed) in the treatment name. Ground chalk was applied as necessary to maintain soil at pH 7 or 6 on sub-plots "a" and "b", respectively; P amount applied to the plots was decreased from 3.5 to 1.7 g P m −2 year −1 since 2017; functional group composition was estimated from botanical separation data between 1915 and 1976, 1991 and 2000, and 2010 and 2012 (e-RA); species richness data report on a 10-year period from 1991 to 2000; total herbage yields and nitrogen nutrition index (NNI) data are given for the last 15 years of the study (2004-2018); phosphorus nutrition index (PNI) and N to P ratios are presented for a 10-year period (2000-2009); soil pH corresponds to samples taken in 1995; NNI was calculated as in Lemaire et al. [28], with parameters for C3 temperate grassland if yield > 100 g m −2 and N crit = 4.8%-N in dry matter if yield < 100 g m −2 . PNI was calculated according to Liebisch et al. [29]. Nutrient limitation was defined based on N to P ratios according to Güsewell [30], with N to P ratios < 10 and > 20 corresponding to N-and P-limitation and N to P ratios between 10 and 20 indicating N and P co-limitation 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 c a .
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 c i /c a . Yield trends over the century diverged for the dicotand grass-rich treatments ( Fig. 3b, P < 0.001). While the    1 (1917-1931). Results are presented for each fertilizer treatment (in the order of the treatment effect on the long-term δ 18 O cellulose -change), the mean of all dicot-rich and grass-rich treatments, and all treatments combined. Bar plots and error bars represent the calculated difference ± SE (SE calculated as SE from period 1 + SE from period 2; n = 11-15 in each period for each treatment depending on data availability; see Additional file 1: Table S1 for details). t tests were used for determining if the long-term changes were statistically significant and significance is presented when P < 0.05. Significance levels: *P < 0.05; **P < 0.01; ***P < 0.001. Significant differences between grass-rich and dicot-rich treatments are designated with different capital letters 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 [22]. Changes of yield (Y) are linked to changes of growing season-integrated water use efficiency via canopy photosynthesis as Y = A (1 -ϕ) (1 -r) [19], with ϕ the proportion of carbon respired and r that allocated to roots and nonharvested or senesced shoot biomass. Growing seasonintegrated canopy photosynthesis of the different treatments may have changed over the century in response to rising c a , 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 (R 2 = 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 δ 18 O of cellulose (δ 18 O cellulose , the "Methods" section) (Fig. 4b, R 2 = 0.52, P < 0.01). Increases of δ 18 O cellulose 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 (δ 18 O rain ) and atmospheric moisture (δ 18 O vapour ) (the "Methods" section). Climate warming should increase δ 18 O rain [36] and has caused an increase of vapor pressure deficit in Europe, which contributes to reduce stomatal conductance [37]. 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 δ 18 O rain (Additional file 1: Fig. S2). And thirdly, the longterm trends of δ 18 O cellulose 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 δ 18 O cellulose between treatments was not attributable to eventual changes of δ 18 O rain or environmental conditions, as all treatments experienced the same site and weather conditions.
On average of all treatments, δ 18 O cellulose 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 δ 18 O cellulose was much stronger in grass-rich (+ 1.1‰; P < 0.001) than dicot-rich (+ 0.3‰; P = 0.06) plots. The trends of δ 18 O cellulose were Relationship between the long-term change in δ 18 O cellulose and the long-term change in stand-scale, growing season-integrated stomatal conductance, estimated as the ratio of C-yield of hay to the CO 2 concentration gradient between the atmosphere and the leaf internal gas space (the "Methods" section) of individual treatments (a), and relationship between relative stomatal conductance change (in % gs change) and δ 18 O-increase (in permil) at the Park Grass Experiment for dicot-rich (n = 6) or grass-rich (n = 6) treatments (mean ± SE) compared to reports of δ 18 O cellulose change versus leafor canopy-scale stomatal conductance change over a range of plant species or genotypes in different environmental conditions (see Additional file 1: Table S2 for details on reported studies) (b). The dashed line and the shaded area in a represent the regression line ± CI95%, while the continuous black line and the shaded area in b represent the mean sensitivity ± CI95% calculated from the reported studies closely proportional (R 2 = 0.70; P < 0.001) to trends in the ratio of hay yield to the CO 2 concentration gradient between the atmosphere and the leaf internal gas space, c a -c i (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 δ 18 O cellulose only occurred when stomatal conductance changed. Additionally, the sensitivity of δ 18 O cellulose 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 δ 18 O cellulose change as a well-grounded measure of longterm 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 δ 18 O cellulose correlated negatively with c i /c a , both at the beginning and end of the century (P < 0.001), with greater ranges of variation in grassrich 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 CO 2 studies, which found consistently decreased rates of N acquisition in grassland, forests and crops under conditions where elevated CO 2 failed to enhance yields [3]. The Park Grass data indicate that such c a -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 δ 18 O cellulose (Fig. 4c, R 2 = 0.56; P < 0.01), in agreement with stomatal conductance-mediated reductions of transpiration and associated mass flowdependent 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 [14], and atmospheric deposition of N and biological N fixation account for a larger proportion of total N acquisition in these plots on average [24]. But, regardless of the mechanism underlying the reduced N acquisition [3] in grass-rich plots, any impairment of N status would reduce photosynthetic capacity and feedback on stomatal conductance to adjust c i /c a 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 CO 2 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 [28] (− 0.08, P < 0.001). Conversely, the nitrogen nutrition index decreased much less in the dicotrich communities (− 0.04, P = 0.01). That difference must have contributed to the convergence of yields of grassand 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 [38]. Much of the Western European grassland receives very high N inputs [39].

Conclusions
The present work provides evidence for an interaction between plant functional groups and the c a -and related climate change-response of temperate-humid permanent grassland in the last century, managed under a haycutting 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 CO 2 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 c a -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 [23], 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 CO 2 enrichment experiments concerning c a -saturation of the CO 2 fertilizer effect on aboveground biomass production [4,6,21], N acquisition [3] 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 [40]. 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 [41].
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 [42], but it changed rapidly following the introduction of the different treatments, creating very contrasting community compositions, which have reached a dynamic equilibrium since 1910 [41]. 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%) dicotrich (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 δ 18 O of α-cellulose (δ 18 O cellulose ), which is associated with interannual variation of meteorological conditions (see the "Stomatal conductance" section and Fig. 5). The subsamples were used to determine the δ 13 C of bulk plant material (δ 13 C p ) and δ 18 O cellulose . Carbon and oxygen isotope composition were expressed as: with R sample the 13 C/ 12 C (or 18 O/ 16 O) ratio of the sample and R standard 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. [43] protocol as modified by Gaudinski et al. [44]. 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 δ 13 C p , the plant material was dried at 40°C for 48 h, ball milled to a fine powder, redried 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 CO 2 gas, previously calibrated against a secondary isotope standard (IAEA-CH6 for δ 13 C, 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 phosphovanadomolybdate colorimetric method following acid digestion [45].

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 (T max ,°C), average air temperature (T mean ,°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 offsetcorrection 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 CO 2 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 (W i ) is a physiological efficiency [47] that represents the ratio of net photosynthesis rate (A) and stomatal conductance to water vapor (g w ): That relationship is controlled by the CO 2 concentration gradient between the atmosphere (c a ) and the leaf internal gas space (c i ), c a -c i , and the ratio of the diffusivities of water vapor and CO 2 (1.6) [19,47]. The CO 2 concentration gradient can also be expressed as the product of c a and the relative CO 2 concentration gradient between the atmosphere and the leaf internal gas space (1 -c i /c a ) (Eq. 5). Thus, solving the equation required two parameters, c a , known from measurements of free air and air bubbles separated from ice cores [25], and c i /c a , that we estimated according to Ma et al. [31] (the bracketed term in Eq. 5). In that term, a is the 13 C discrimination during diffusion of CO 2 in air through the stomatal pore (4.4‰), a m (1.8‰) the fractionation associated with CO 2 dissolution and diffusion in the mesophyll, and b (29‰) and f (11‰) the fractionations due to carboxylation and photorespiration, respectively; Г* is the CO 2 compensation point in the absence of mitochondrial respiration calculated following Brooks and Farquhar [48], and g s /g m the ratio of stomatal and mesophyll conductance [31]. Relying on theory developed by Farquhar et al. [49] and Farquhar and Richards [19], this model accounts for the effects of mesophyll conductance and photorespiration on Δ 13 C-based estimations of intrinsic water use efficiency. Traditionally, estimations of W i from Δ 13 C have used a more abbreviated version of the Farquhar model of photosynthetic carbon isotope discrimination (Δ 13 C) in C 3 plants. That version neglected the photorespiration term and was based on the simplifying assumption that mesophyll conductance is infinite. However, W i estimated from Δ 13 C using the abbreviated model systematically overestimates gas exchange-based estimates of W i , an error that can be corrected by using a constant g s /g m 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 [31]. Also, global scale effects of mesophyll conductance and photorespiration of C 3 vegetation on 13 C discrimination of the terrestrial biosphere have been evident for the last four decades of atmospheric CO 2 increase and concurrent change of the 13 CO 2 / 12 CO 2 ratio [50]. Estimates of intrinsic water use efficiency at Park Grass calculated following Ma et al. [31] were closely correlated (R 2 = 0.97) with those presented by Köhler et al. [23] using the abbreviated W i model for the same treatments and time spans.
Δ 13 C was obtained from the measured δ 13 C p of samples and the δ 13 C of atmospheric CO 2 (δ 13 C a ), as Estimates of c a and δ 13 C a were obtained as in Wittmer et al. [51] and Köhler et al. [52]. The time series generated by Köhler et al. [52] with c a and δ 13 C a 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 c a and δ 13 C a estimation. Calculations of Δ 13 C and W i were performed for each treatment, individually.
As intrinsic water use efficiency was estimated from a representative sample of a hay harvest, it represents a growingseason assimilation-and allocation-weighted measure.

Stomatal conductance
Extrapolating from elevated CO 2 studies [7], we hypothesized that the atmospheric CO 2 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 [23] 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 CO 2 concentration and 13 C 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 [19], canopy-scale, growing season-integrated stomatal conductance to CO 2 (g s , 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 nonharvested 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 (%g s -change) between the beginning (period 1, subscript 1) and end (period 2, subscript 2) of the century for every treatment as: %g s -change ¼ ððg s 2 − g s 1 Þ=ðg s 2 ÞÞ Â 100: ð8Þ Note that this procedure eliminated the constants for ϕ and r, so that the estimated %g s -changes were independent of ϕ and r, except if they had changed during the century.
The effect of stomatal conductance on Δ 18 O cellulose 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 18 O-enriched during daytime and imprints its 18 O-signal on photosynthetic products used in cellulose synthesis [62], and (3) evaporative 18 O-enrichment of leaf water decreases with increasing stomatal conductance [63]. The δ 18 O cellulose 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 [64]. However, the relation between stomatal conductance and Δ 18 O cellulose is influenced by multiple morpho-physiological vegetation characteristics and environmental variables [65], so that stomatal conductance cannot be simply and quantitatively estimated from Δ 18 O cellulose or δ 18 O cellulose at present [66]. As no direct measurements of stomatal conductance were available at Park Grass during the last century, we compared δ 18 O cellulose 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 δ 18 O cellulose . Further, we compared the δ 18 O cellulose change versus %g s -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 [67]. In this analysis, we also included (1) the relationship between Δ 18 O cellulose and canopy-scale stomatal conductance predicted by a process-based 18 O-enabled soil-plantatmosphere model for a multi-seasonal data set of grassland [64] and (2) an unpublished data set from monocultures of Lolium perenne grown at CO 2 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 [68].
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 δ 18 O cellulose ) 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 %g s -decrease (Eq. 8) per unit increase in δ 18 O cellulose . 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 δ 18 O cellulose relationship of that study. On average of all the studies, stomatal conductance decreased by 39% per 1‰ increase of δ 18 O cellulose .

Statistical analysis
We tested the long-term response of intrinsic water use efficiency, yield, δ 18 O cellulose , 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 δ 18 O cellulose ) was tested using ordinary leastsquares 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 [69]. The R package ggplot2 [70] was used for data plotting.