Last-Century Increases in Intrinsic Water-Use Efficiency of Grassland Communities Have Occurred over a Wide Range of Vegetation Composition, Nutrient Inputs, and Soil pH1[OPEN]

Intrinsic water-use efficiency of grasslands increased over a wide range of nutrient inputs, soil pH, and plant community compositions during the last century. Last-century climate change has led to variable increases of the intrinsic water-use efficiency (Wi; the ratio of net CO2 assimilation to stomatal conductance for water vapor) of trees and C3 grassland ecosystems, but the causes of the variability are not well understood. Here, we address putative drivers underlying variable Wi responses in a wide range of grassland communities. Wi was estimated from carbon isotope discrimination in archived herbage samples from 16 contrasting fertilizer treatments in the Park Grass Experiment, Rothamsted, England, for the 1915 to 1929 and 1995 to 2009 periods. Changes in Wi were analyzed in relation to nitrogen input, soil pH, species richness, and functional group composition. Treatments included liming as well as phosphorus and potassium additions with or without ammonium or nitrate fertilizer applications at three levels. Wi increased between 11% and 25% (P < 0.001) in the different treatments between the two periods. None of the fertilizers had a direct effect on the change of Wi (ΔWi). However, soil pH (P < 0.05), species richness (P < 0.01), and percentage grass content (P < 0.01) were significantly related to ΔWi. Grass-dominated, species-poor plots on acidic soils showed the largest ΔWi (+14.7 μmol mol−1). The ΔWi response of these acidic plots was probably related to drought effects resulting from aluminum toxicity on root growth. Our results from the Park Grass Experiment show that Wi in grassland communities consistently increased over a wide range of nutrient inputs, soil pH, and plant community compositions during the last century.

The intrinsic water-use efficiency (W i ) of plants is controlled by photosynthetic carbon assimilation and stomatal conductance via the leaf-level coupling of CO 2 and water fluxes. A general, but variable, increase of W i under rising atmospheric CO 2 has been observed in long-term studies (Peñuelas et al., 2011;Franks et al., 2013;Saurer et al., 2014), but little is known about other environmental or ecosystem factors, which may interact with the effect of increasing CO 2 on W i . An improved understanding of putative interactive mechanisms is important because changes in W i may have significant effects on the global terrestrial carbon and water cycles (Gedney et al., 2006;Betts et al., 2007). This study explores the interactive effects of the increase in atmospheric CO 2 (observed over the last century), nutrient loading, and soil pH together with other related effects on plant species richness and functional group composition on the coupling of plant CO 2 and water fluxes in a seminatural grassland in southeastern England.
W i is a leaf-level efficiency that has also been termed potential water-use efficiency or physiological water-use efficiency, as it excludes the direct influence of vapor pressure deficit (VPD), a parameter determined by environmental conditions, on leaf-level water-use efficiency (Farquhar et al., 1989;Franks et al., 2013). W i reports the relationship between net CO 2 assimilation rate (A n ) and stomatal conductance for water vapor (g H2O ): According to the first law of Fick, A n can be given as the product of the stomatal conductance for CO 2 (g CO2 ) and the concentration gradient between the atmosphere (c a ) and the leaf internal gas space (c i ): A n = g CO2 (c ac i ).
Using g CO2 (c ac i ) instead of A n in Equation 1, replacement of g H2O /g CO2 by the numerical value of g H2O /g CO2 (1.6) and rearrangement yields the following alternative expression of W i : Equation 2 reveals that past changes of W i must have been controlled by two parameters: the change of c a and the concurrent change of 1c i /c a , the relative gradient for CO 2 diffusion into the leaf (Franks et al., 2013). A change in the relative gradient is determined by the changes in A n relative to g H2O , as leaves respond to changing c a and other environmental factors. In particular, Equation 2 shows that any variation in the climate change response of W i is determined by the c i /c a response, if the comparison is made for vegetation at the same location and in the same period of time. Studies with C 3 vegetation, including trees/forests and C 3 grasslands, have revealed a general increase of W i in the last century (Bert et al., 1997;Duquesnay et al., 1998;Feng, 1999;Arneth et al., 2002;Saurer et al., 2004;Barbosa et al., 2010;Köhler et al., 2010;Andreu-Hayles et al., 2011). In many cases, c i /c a , estimated by 13 C discrimination (Farquhar et al., 1989), varied relatively little. Indeed, it has been suggested, based on theoretical grounds and empirical evidence from studies over geological/evolutionary to short time scales, that adaptive feedback responses will tend to maintain c i /c a approximately constant (Ehleringer and Cerling, 1995;Franks et al., 2013), as plants optimize carbon gain with respect to water loss (Cowan and Farquhar, 1977). Yet, c i /c a -dependent variation in the W i response to climate change has also been noted (Peñuelas et al., 2011;Köhler et al., 2012) over the last century, indicating that additional factors, perhaps including other global change drivers, can modify the W i response over this time scale, at least transiently. A meta-analysis by Peñuelas et al. (2011) reports c i /c a -dependent increases of W i for different forests between 6% and 36% from the early 1960s to 2000s. A recent study by Saurer et al. (2014) on European forest trees found increases in W i ranging from 1% to 53% during the last century. The strongest increase of W i was recorded in regions where summer soil-water availability decreased in the last century. For different grassland communities, the c i /c a -dependent increases of W i varied between 13% and 28% at one site (Köhler et al., 2012(Köhler et al., ) from 1915(Köhler et al., to 2009. Evidently, such variation can have important repercussions for the coupling of terrestrial CO 2 and water fluxes. Yet, little is known about the mechanism(s) underlying the variation.
At the Park Grass Experiment (PGE) at Rothamsted, England, Köhler et al. (2012) observed a nitrogen supplydependent enhancement of the W i response on plots receiving nitrate fertilizer and maintained at a nearneutral soil pH by liming. However, the actual relationship between nitrogen supply and W i response did not hold when the unlimed control (soil pH approximately 5.2) was included in the comparison. Remarkably, however, there was a significant positive relationship between the grass content of the community and the W i response of the experimental plots in the investigation. These results suggested that the effect of nutrient supply on the W i response of the grassland communities was indirect, perhaps working via effects on soil pH and/or vegetation composition (plant species richness or functional group composition).
The PGE provides a unique opportunity to study century-scale variation in the c i /c a -dependent variation of W i for a wide range of diverse grassland communities. Much of the extant ecosystem-scale variability of plant species richness and soil pH in temperate grasslands of Europe (Ceulemans et al., 2014) is included in the range of plot-scale plant species richness and soil pH at the PGE (which is reported in this investigation). The different long-term applications of fertilizer and lime over the past century have resulted in substantial changes in soil pH, species richness, and grass content on the experimental plots, but in most cases, within-plot changes over the study period considered here ) were comparatively small (Crawley et al., 2005;Silvertown et al., 2006). All experimental plots are located at the same site and are exposed to the same weather conditions. Consequently, trends in climate as a direct driver for differences in W i between plots can be ruled out.
Here, we explore putative mechanisms underlying eventual c i /c a -dependent variation of W i during the last century at the PGE by, first, quantifying the sustained effect of a wide range of contrasting fertilizer treatments (n = 16) on the change of W i during the last century and, second, analyzing the relationships between the observed W i response of treatments and the respective nutrient status, soil pH, plant species richness, and plant functional group composition of the grassland communities.

Long-Term Changes of W i and c i /c a
For the different treatments, the 15-year average W i (⌀W i ) ranged between 49.6 and 62.5 mmol mol 21 in the 1915 to 1929 period and between 60.2 and 76.8 mmol mol 21 in the 1995 to 2009 period (Table II). The change of W i (DW i ) between the two time periods was positive for all treatments, ranging between +5.9 and +14.7 mmol mol 21 (Fig. 1A). This effect was significant for all treatments (P , 0.001, adjusted for multiple comparisons) and corresponded to a relative enhancement of W i (⌀W i 1995-2009 / ⌀W i 1915-1929 ) of 11% to 25% between the two periods ( Table II).
The treatment-specific mean c i /c a varied between 0.67 and 0.74 in the 1915 to 1929 period, the same as in the 1995 to 2009 period (Table II). Treatment effects on changes of c i /c a between the 1915 to 1929 and 1995 to 2009 periods were not significant, except for the unlimed control (CONTROL.U; for treatment definitions, see "Materials and Methods" and Table I), where the increase was significant (P , 0.001), although very small (+0.028).
Treatment Effects on DW i : Nitrogen Input, Nitrogen Form, and Soil pH The highest increases of W i were observed in the unlimed NH 4 + treatments (N1PK.U, N2PK.U, and N3PK. U: +13.8, +14.2, and +14.7 mmol mol 21 ), and the smallest increase was observed in the limed control (CONTROL. L: +5.9 mmol mol 21 ). Significant differences between treatments existed between the limed control (CONTROL.L) and two of the unlimed NH 4 + plots (N2PK.U and N3PK.U; P , 0.05, adjusted for multiple comparisons).
In general, no significant effect of nitrogen input on DW i was found at the 5% a level ( Fig. 2A). Even within the NH 4 + treatments with the largest differences in nitrogen input (48, 96, and 144 kg ha 21 a 21 ), no significant effect of the amount of applied NH 4 -nitrogen on DW i was evident. Also, there were no significant differences between mean DW i when all treatments were grouped by nitrogen form (ANOVA, P = 0.06). However, when the limed and unlimed NH 4 + treatments were analyzed separately, a significant difference was found between the unlimed NH 4 + treatments and the NO 3nitrogen or no-nitrogen treatments (ANOVA, P , 0.01).
Increasing soil pH led to a significant decrease of DW i ( Fig. 2B; DW i = 16.7 -1.1x; r 2 = 0.24, P , 0.05, degrees of freedom [df] = 14). Furthermore, NH 4 + treatments showed a significant difference between limed and unlimed plots (P , 0.01, adjusted for multiple comparisons). But such a difference was not observed in the other treatments. No significant relationships were found when the unlimed NH 4 + plots were omitted from the analysis.

Species Richness
DW i was negatively related to species richness ( Fig.  2D; DW i = 13.6 -0.1x; r 2 = 0.43, P , 0.01, df = 14), with Table I. Summary of studied treatments at the PGE grouped by fertilizer nitrogen form (no nitrogen, nitrate-nitrogen, and ammonium-nitrogen) Treatment abbreviations and main plot numbers (subplots after 1965 in parentheses), annual nitrogen, phosphorus (P) and potassium (K) fertilizer application (in kg ha 21 ; nitrogen includes estimated biologically fixed nitrogen [N bf ]), functional group composition (in percentage dry weight for grasses [G], nonlegume forbs [F], and legumes [L] in harvested biomass), species richness, and soil pH (in water) are shown. Limed and unlimed subplots are distinguished by L and U in the treatment abbreviation. b Calculated as total number of species counted during a 10-year period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). c Soil pH in water was calculated as the average from available measurements in the sampling period (Warren and Johnston 1964;Williams, 1978;Rothamsted Research, 2006). variation in richness explaining 43% of the total variation of DW i . Significant and quantitatively similar relationships between DW i and species richness were also found when only the unlimed (DW i = 14.4 -0.1x; r 2 = 0.64, P , 0.05, df = 6) or the NH 4 + (DW i = 13.8 -0.1x; r 2 = 0.55, P , 0.05, df = 8 [no-nitrogen treatments included]) treatments were considered. All the observed regressions were strongly dependent on the inclusion of the species-poor unlimed NH 4 + treatments (two to six species per plot) in the analysis. No significant relationship with species richness was evident when these treatments were omitted. The other treatments showed a broad spectrum of species richness (16-46), with no significant trend of DW i versus number of plant species.

Correlations between Variables
Correlations between the directly controlled variables (nitrogen input and soil pH) and the indirectly manipulated variables (percentage abundance of grass and plant species richness) were highly significant, indicating that they were not independent of each other (Table III). Nitrogen input and soil pH were not correlated, as both variables were manipulated separately.

W i Showed a Variable Increase between Different Grassland Communities at a Single Site
In principle, the variability of DW i responses during the last-century CO 2 increase can be attributed either to inherent differential physiological responses (control of A n or g H2O ) of different species and the communities they form (direct effects) or to interacting effects of other environmental factors that modify the aforementioned community-scale physiological responses (indirect effects). The observed range of DW i of the present grassland communities (n = 16) was smaller (11%-25%, average 19%) than that reported by Saurer et al. (2014) for a variety of forest trees in a wide range of sites across Europe (1%-53%, average 28%; n = 35) in a Figure 1. A, Absolute increases in W i (mmol mol 21 ) between the two studied periods (DW i = ⌀W i [1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009] -⌀W i [1915][1916][1917][1918][1919][1920][1921][1922][1923][1924][1925][1926][1927][1928][1929] ) in ascending order B, Annual nitrogen input as the sum of fertilizer nitrogen and biologically fixed nitrogen (N bf ) in kg ha 21 . C, Soil pH in water (circles) and species richness (bars). D, Contribution of grass biomass to total harvested sward biomass in percentage. Error bars represent SE. similar time period (1901-1910 and 1991-2000). The generally lower W i and lower DW i of the grasslands at the PGE compared with forests may be related to the observed lower carbon cost per unit of water used in C 3 grasses compared with gymnosperm and angiosperm trees (Lin et al., 2015). The strongest increases in tree DW i were related to decreasing soil water content at the respective sites. Being located at a single site, the observed treatment-specific variation of DW i at the PGE during the past 100 years is unrelated to direct effects of weather factors or other characteristics of geography. Instead, all the evidence for mechanisms controlling the differences in DW i between treatments at the PGE point to site conditions and plant community characteristics that were modified by the (differential) nutrient and lime inputs. These mechanisms may include the controls of plant community assembly (Silvertown, 1980;Crawley et al., 2005), modifications of the genetic diversity of plant species and niche dimensionality (Harpole and Tilman, 2007), soil microorganism composition (Rousk et al., 2011), plant-soil microbe interactions in processing/sequestration of plant carbon (Fornara et al., 2011), alterations of the relative availabilities of (both fertilized and unfertilized) macronutrients and micronutrients (White et al., 2012), and, possibly, aluminum toxicity (Davies and Snaydon, 1973). Relationships between DW i (mmol mol 21 ) and annual nitrogen input as the sum of fertilizer nitrogen and N bf in kg ha 21 (A), soil pH in water, y = 16.7 + 1.1x, r 2 = 0.24, P , 0.05 (B), contribution of grass biomass to total harvested sward biomass in percentage, y = 3.6 + 0.1x, r 2 = 0.36, P , 0.01 (C), and species richness, y = 13.6 + 0.1x, r 2 = 0.43, P , 0.01 (D). Simple linear regression lines are only shown where significant (P , 0.05). Symbols represent the type of nitrogen addition (triangles = no nitrogen, circles = NO 3 2 , and squares = NH 4 + ) and liming treatment (white symbols = unlimed and black symbols = limed). Bars denote the SE. Horizontal error bars in B and C are not shown where smaller than the symbols. The DW i was significantly different from zero in all treatments (P , 0.001). A significant increase of c i /c a was only observed in the CONTROL.L treatment (P , 0.001).

Treatment
Average SE 1915SE -1929SE 1995SE -2009SE 1915SE -1929SE 1995SE -2009  Although the variability of DW i resulted from the diversity of fertilizer and liming treatments, we did not detect any direct, general relationship between a single nutrient factor, such as annual nitrogen inputs (including biologically fixed nitrogen) or phosphorus + potassium additions, on DW i . For every nutrient (or lime) added or not (singly or in combination with other nutrients), the observed relative range of W i increase varied widely. Moreover, investigations of the plot-specific type of nutrient limitation (phosphorus or nitrogen) for biomass production in the last 50 years revealed no relationship with DW i (R. Hirl, I.H. Köhler, A. Macdonald, and H. Schnyder, unpublished data). Conversely, soil pH, as controlled via targeted amounts of applied lime, was related to DW i , with the highest DW i observed on the most acidic soils. Similar relationships of DW i with species richness and the relative abundance of grasses indicate that this treatment effect may also have worked indirectly via changes in plant community composition.

Does Aluminum Toxicity Lead to Higher DW i via Increased Drought Stress?
The greatest relative increase of W i occurred on the plots with most acidic soils (pH approximately 3.5-4), which are a result of long-term annual ammonium sulfate additions. The low pH on these plots, which has effected a significant mobilization of aluminum (Gould et al., 2014) and other heavy metals, has also led to a dramatic loss of (acid-sensitive) plant species and the selection of a very small number of acid-tolerant and metal-resistant grasses (e.g. Holcus lanatus, Anthoxanthum odoratum, and Agrostis capillaris; Blake and Goulding, 2002;Crawley et al., 2005;Silvertown et al., 2006). These species-poor, grass-rich plots have shown much greater interannual variability of biomass yields than speciesrich plots, perhaps due to the better buffering of species-rich communities against climatic variability (Dodd et al., 1994b). High aluminum concentrations may have exacerbated the effect of climatic variability by causing greater soil moisture deficits on the acidic plots; aluminum toxicity is known to inhibit root growth (Kochian, 1995;Horst et al., 2010), leading to reduced root length density and shallower rooting depth (Jentschke et al., 2001;Tang et al., 2001) in a variety of species. Physiological relationships between aluminum toxicity, impeded root growth, and drought sensitivity of plants are well established (Yang et al., 2013). Also, Dodd et al. (1994b) suggested that the large yield variability of the species-poor, acidic plots on the PGE were related (at least in part) to soil moisture deficits. There is no evidence that the dominant species (H. lanatus, A. odoratum, and A. capillaris) on the acidic plots differ inherently in their reactivity to drought; indeed, their Ellenberg moisture indicator values (Hill et al., 1999) do not differ from the dominant grass species (Alopecurus pratensis, Arrhenatherum elatius, and Dactylis glomerata) on the other plots. However, effects of soil moisture deficits, resulting from aluminum effects on rooting depth, may interact with increasing atmospheric CO 2 at the level of stomata, enhancing the increase in W i via a drought stress-induced decrease in g H2O . Accordingly, soil acidification in combination with aluminum toxicity may have a similar effect on DW i to the climate-related decrease in water availability reported by Saurer et al. (2014).

On the Mechanism and Implication of Treatment-Specific W i Responses
W i is controlled via c i /c a (Eq. 2). This physiological set point seems to have been maintained approximately constant under increasing c a (Ward et al., 2005;Franks et al., 2013), which leads to a concurrent increase of W i . Mean atmospheric CO 2 concentration has increased by 23.5% between the first (1915-1929) and second (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) period of this study. In comparison, W i increased by 11% to 25% between the two periods. This means that c i /c a was approximately maintained constant in the treatments with the highest increase of W i (the species-poor, grass-rich communities on the acidic plots), whereas c i /c a increased to some extent on the plots with the least increase of W i (the species-and forb-rich communities on plots with neutral pH). This (small but significant) increase in c i /c a more than halved the CO 2 -driven enhancement of W i on the limed control treatment.
Despite a large variation in nutrient supply, soil pH, plant species richness, and functional group composition, the majority of the treatments showed intermediate increases in W i (Fig. 1A). Only the most extreme plots (acidic soil and nutrient poor) differed significantly from the others. This general response suggests that a wide range of European seminatural C 3 grasslands may show increasing W i under rising atmospheric CO 2 by maintaining c i /c a approximately constant or increasing it just slightly. Depending on changes in atmospheric VPD, the increases in W i could translate directly into a corresponding response of leaf-level transpiration efficiency (W t ), as W t = W i 3 VPD 21 (Farquhar and Richards, 1984;Farquhar et al., 1989). During spring growth (March-June), VPD at the PGE did not change between the two studied periods. If this is also true for other European regions, then, in combination with the observation of constant c i /c a , this implies that ecosystem evapotranspiration in European seminatural Table III. Correlations (Pearson coefficients) for average percentage abundance of grasses in total harvested herbage biomass, species richness, soil pH, and nitrogen input (calculated as fertilizer nitrogen plus N bf ) Asterisks denote the significance level (Pearson's product-moment correlation: ***, P , 0.001; **, P , 0.01; and n.s., not significant).

CONCLUSION
This study demonstrates significant correlations between the last-century increase of W i (DW i ) and the proportion of grasses in the plant community, species richness of communities, and soil pH, which were all interdependent site-community characteristics resulting from long-term differential fertilizer and liming treatments. Remarkably, however, the correlations depended very strongly on the inclusion of the ammonium sulfate fertilizer treatments, which led to the highest DW i and was associated with low soil pH, a high grass proportion, and low species richness (associated with the selection of aluminum-tolerant grass species). Therefore, drought sensitivity, resulting from the inhibition of root growth by aluminum at low soil pH, seemed to be the most important interactive factor causing an above-average last-century increase of W i at the PGE.
All grassland communities exhibited a consistent increase in W i , and this was true for communities with greatly differing plant species richness and functional group composition and nutrient and lime inputs during the last century. The analysis covered virtually the entire range of species richness and soil pH at the PGE, which is representative of a significant proportion of the extensively managed (two cuts per year) European temperate seminatural grasslands. Based on this observation, we suggest that the increase in W i on the PGE is representative of a large part of the extensively managed temperate grassland in Europe and perhaps other parts of the world. To determine the extent to which this also applies to more intensively managed grassland would require further study.

The Park Grass Experiment (PGE)
The PGE was established in 1856 at Rothamsted Experimental Station (40 km north of London in Hertfordshire, England, 0°219W, 51°499N, 128 m above sea level). Different fertilizer and manure treatments were applied to 20 experimental plots on 2.8 ha of old grassland. In 1903, most of the plots were divided in two to introduce a test of lime on one half; 4 tons ha 21 CaCO 3 was applied every fourth year. In 1965, most of the plots were further divided into four subplots (a, b, c, and d). The a, b, and c subplots receive lime every 3 years, if necessary, to maintain a target soil pH of 7, 6, and 5, respectively. The d subplots are unlimed (for more details, see Köhler et al., 2012). The experiment is located on a moderately well-drained silty clay loam overlying clay-with-flints. The soil pH was slightly acidic when the experiment began (5.4-5.6), and the nutrient status was poor (Silvertown et al., 2006).
The original vegetation of the PGE has been classified by Dodd et al. (1994a) as a dicotyledon-rich Cynosurus cristatus-Centaurea nigra grassland, one of the mesotrophic grassland communities in the British National Vegetation Classification system. Since the start of the experiment, the herbage has been cut around mid-June for hay production; samples were taken from the material dried in situ until 1960. In subsequent years, strips were cut with a forage harvester, and fresh herbage samples were taken for drying and archiving. The herbage remaining on the plot was made into hay, as before. Originally, the following regrowth was grazed by sheep, which were penned on the individual plots, but since 1875, grazing was abandoned and a second harvest cut was removed green. Dried samples from all plots and both cuts have been stored in the Rothamsted Sample Archive since the beginning of the experiment (Silvertown et al., 2006). This work presents data from the first cut of two study periods: 1915 to 1929 and 1995 to 2009. Rainfall and temperature have been measured at Rothamsted since 1853 and 1878, respectively. In the months preceding the first cut (March-June), mean daily temperature 6 SD was 9.2°C 6 0.5°C versus 10.5°C 6 0.6°C and VPD 6 SD was 0.32 6 0.04 kPa versus 0.34 6 0.04 kPa in the two study periods (1915-1929 versus 1995-2009). The average sum of annual rainfall 6 SD in 1915 to 1929 versus 1995 to 2009 was very similar (805 6 140 mm versus 796 6 146 mm), and this was also the case for seasonal rainfall from March to June (208 6 49 mm versus 219 6 62 mm). For more climatic information, see Köhler et al. (2010).

Fertilizer and Liming Treatments
A previous study (Köhler et al., 2012) investigated a smaller range of treatments (n = 5), including nitrate fertilizer inputs, on plots where a nearneutral soil pH was maintained by liming. This work extends the range of treatments (n = 16) and now reports on both limed and unlimed treatments combined with or without additions of phosphorus (35 kg ha 21 a 21 from triple superphosphate) plus potassium (225 kg ha 21 a 21 from potassium sulfate) and one of two forms of nitrogen fertilizers given at elemental rates of 0, 48, or 96 kg ha 21 a 21 (as sodium nitrate or ammonium sulfate) or 144 kg ha 21 a 21 (as ammonium sulfate). The treatments are designated thus: CONTROL denotes unfertilized plots; PK indicates plots receiving phosphorus plus potassium fertilizer; N* refers to nitrate and N to ammonia fertilizer; and 1, 2, or 3 give the dosage (48, 96, or 144 kg ha 21 a 21 ) of nitrogen fertilizer. In the following, limed and unlimed subplots are distinguished by L and U. Thus, for instance, N*1PK. U stands for the treatment receiving 48 kg ha 21 a 21 nitrogen in the form of sodium nitrate and phosphorus plus potassium, but no lime. For details of treatments and the abbreviations used, see Table I. Most treatments have not been replicated or randomized. However, the meadow was reasonably uniform before the experiment began, and the size of the plots (greater than 100 m 2 ) compensates to some extent for the lack of replication (Crawley et al., 2005). The experiment included two replicate plots for the CONTROL.U and PK.U treatments for the whole studied period and for the N*1.U and N*1PK.U treatments for the period after 1965.
Fertilizer nitrogen was applied in spring as one dressing except for the ammonium and nitrate treatments with the largest nitrogen inputs (N3PK.U, N3PK.L, N*2PK.U, and N*2PK.L), where the dressing was split. The nitrogen applied to the N3PK.U and N3PK.L plots was split in two, 96 kg ha 21 was applied on the first date and 48 kg ha 21 on the second, about 4 weeks after the first. Similarly, the nitrogen applied to the N*2PK.U and N*2PK.L plots was also split, 48 kg ha 21 being applied on both dates. Since 1980, these split dressings have stopped, and all of the plots given nitrogen received the total amount in one dressing. The other fertilizers were applied in winter. Nitrogen input (Table I; Fig. 1B) was defined as fertilizer-nitrogen supply plus an estimate of biologically fixed nitrogen (N bf ). The average N bf input for the spring growth was estimated using the model of Høgh-Jensen et al. (2004) with the parameterization for cuts of 1-to 2 year-old mixed grass and red clover (Trifolium pratense) systems as described by Köhler et al. (2012). Estimates of N bf were included for treatments with an average nitrogen fixation of greater than 10 kg ha 21 a 21 , which included treatments PK.L, PK.U, N*1.PK.L, and N*1.PK.U (Table I; Fig. 1B). Results from soil pH measurements (measured in water) were available for the years 1923 and 1959 (Warren and Johnston, 1964), 1974to 1977(Williams, 1978), and 1991to 2000(Rothamsted Research, 2006. Average soil pH ( Fig. 1C; Table I) in the sampling period differed greatly between treatments but was relatively constant within a treatment, except for the limed ammonium sulfate plots, where pH increased over time.

Functional Group Composition and Species Richness
Functional group composition and species richness differed greatly between treatments (Table I; Fig. 1, C and D), but according to Silvertown et al. (2006), functional group composition within treatments has reached a dynamic equilibrium and was relatively constant since 1915. Between years, the contribution of different functional groups has varied widely within many plots, but no general long-term trend was observed (Supplemental Fig. S1). Therefore, vegetation composition at the functional group level (grasses, nonlegume forbs, and legumes) was assumed to be constant within treatments, and means for each group were calculated from partial botanical separation data available for individual years from 1915 to 1976 (data set PARKPARTCOMP, eRA, 2013) and for each year from 1991 to 2000 (data set PARKCOMPIC, eRA, 2013).
We compared species richness (the number of species present) between the 1915 to 1929 period and the 1995 to 2009 period using available species-level botanical separation data (data sets PARKCOMP and PARKCOMPIC, eRA, 2013). For the first period, however, botanical separation data at the species level were only available for a single year on some treatments. For estimation of the average differences in species richness between the treatments, therefore, we used data from a comprehensive study (Crawley et al., 2005;data set PARKCOMPIC, eRA, 2013), which reflects the total number of species counted during a 10-year period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000). Slight changes over time between the two studied periods were observed on the PK.U, N1.U, and N2PK.U treatments but were small in comparison with between-treatment variation (Supplemental Fig. S2).
Estimation of W i from c a and Carbon Isotope Composition W i was obtained from Equation 2. For c a , we used published data of c a obtained from measurements of CO 2 in free air (Keeling et al., 2009) and in gas bubbles from ice cores for the time before CO 2 measurements in free air began (Friedli et al., 1986;Francey et al., 1999). c i /c a was estimated from the known relationship with carbon isotope discrimination, 13 D (Farquhar et al., 1982;Farquhar and Richards, 1984) as obtained from the linear model of 13 D: where a is the fractionation of 13 CO 2 (relative to 12 CO 2 ) during diffusion through the stomatal pores (4.4‰) and b is the net fractionation by carboxylases (27‰; Eq. 2). 13 D was obtained from the carbon isotope composition of the archived plant material (d 13 C p ) and that of atmospheric CO 2 (d 13 C a ; Friedli et al., 1986;Francey et al., 1999;White and Vaughn, 2011) as The estimates for atmospheric d 13 C and CO 2 concentration were obtained for each year in the two sampling periods as described by Köhler et al. (2010). From 1915 to 2009, atmospheric CO 2 increased by 96 mmol mol 21 . The average difference between the two 15-year-long periods was 72 mmol mol 21 . The mean atmospheric CO 2 concentration 6 SD during the first period was 306 6 1 mmol mol 21 , compared with 378 6 8 mmol mol 21 in the second period.

Carbon Isotope Analysis
Sample preparation and carbon isotope analysis were performed with the same procedures/protocols and instrumentation as described by Köhler et al. (2012). Representative subsamples of plant material were collected in the Rothamsted Sample Archive, dried at 40°C for 48 h, milled to homogeneity, and dried again at 60°C for 24 h. Aliquots of 0.7 6 0.05 mg were weighed into tin cups (IVA Analysentechnik) and combusted in an elemental analyzer (NA 1110; Carlo Erba) interfaced (Conflo III; Finnigan MAT) with an isotope ratio mass spectrometer (Delta Plus; Finnigan MAT). Carbon isotope data were obtained as d 13 C, with d 13 C = [(R sample /R standard ) -1], where R is the 13 C/ 12 C ratio in the sample or standard (Vienna Pee Dee Belemnite, V-PDB). Each sample was measured against a laboratory working standard CO 2 gas, which was previously calibrated against a secondary isotope standard (IAEA-CH6, International Atomic Energy Agency, Sucrose Reference Material; accuracy of calibration 6 0.06‰ SD). After every 10th sample, a solid internal laboratory standard with similar carbon-nitrogen ratio to the respective sample material (fine ground wheat [Triticum aestivum] flour) was run as a control. The solid internal laboratory standards were previously calibrated against the IAEA-CH6 standard. The precision for sample repeats was better than 0.1‰.

Statistical Analysis
The climate change response of W i , DW i , of the different plots was obtained as the difference between the average of W i in the 1915 to 1929 period (⌀W i 1915-1929 ) and the 1995 to 2009 period (⌀W i 1995-2009 ), 1995-2009 2 ⌀W i 1915-1929 : The length of the two periods was chosen so that it was long enough to absorb the interannual variability of climatic conditions and short enough to still retain a significant difference in CO 2 concentration between the two periods. As a proxy for the interannual climatic variability, we used the plant available soil water for a grass reference crop (PAW ref ), which has been shown to explain interannual and intraannual variation in 13 D better than single climatic factors like VPD or precipitation (Schnyder et al., 2006;Köhler et al., 2010). PAW was calculated as PAW i ¼ PAW i 2 1 þ P i 2 AET i , where PAW i-1 is the modeled plant available water of the previous day, P i is the precipitation on day i, and AET i is the modeled actual evapotranspiration on day i. AET i equaled PET i as long as PAW i /PAW capacity . 0.3. Otherwise, AET i was calculated as AET ¼ PET 0:3 $ PAWi PAWcapacity (Schnyder et al., 2006). PET was estimated with the FAO Penman-Monteith equation for a standard grass reference crop (Allen et al., 1998). The maximum plant available soil water (PAW capacity ) of the soil has not been measured directly but was inferred to be 135 mm (for the top 70 cm of the soil) from measurements on similar soils at Rothamsted (Avery and Catt, 1995). For further details on the calculation of PAW ref , see Köhler et al. (2012).
We calculated the mean and SD of PAW ref for periods of increasing duration (2-30 years) after 1915 and before 2009 and compared the values between the periods. This procedure led to the choice of a period length of 15 years, for which the mean and SD of both periods were similar (PAW ref1915-1929 , 55.6 6 30.7 mm; PAW ref1995-2009 , 58.7 6 28.3 mm).
For each treatment, we tested if DW i was significantly different from zero using linear regression, adjusting for multiple comparisons. The relationships between DW i and direct and indirect treatment parameters (nitrogen input, soil pH, grass content in percentage, and species richness) were analyzed using simple linear regression. When analysis was conducted on subsets grouped by nitrogen form (NH 4 or NO 3 treatments), the no-nitrogen treatments were included in the regressions as controls. All statistical analysis was done with R (R Core Team, 2014). The multcomp package (Hothorn et al., 2008) in R was used for multiple comparisons of model parameters. Adjusted P values (Holm method) are reported for multiple comparisons.

Supplemental Data
The following supplemental materials are available.
Supplemental Figure S1. Percentage contribution of functional groups to total harvested dry matter on the studied treatments Supplemental Figure S2. Species counts in individual years on the studied treatments