Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Long-term trends in yield variance of temperate managed grassland

Long-term trends in yield variance of temperate managed grassland The management of climate-resilient grassland systems is important for stable livestock fodder production. In the face of climate change, maintaining productivity while minimizing yield variance of grassland systems is increasingly challenging. To achieve climate-resilient and stable productivity of grasslands, a better understanding of the climatic drivers of long- term trends in yield variance and its dependence on agronomic inputs is required. Based on the Park Grass Experiment at Rothamsted (UK), we report for the first time the long-term trends in yield variance of grassland (1965–2018) in plots given different fertilizer and lime applications, with contrasting productivity and plant species diversity. We implemented a statisti - cal model that allowed yield variance to be determined independently of yield level. Environmental abiotic covariates were included in a novel criss-cross regression approach to determine climatic drivers of yield variance and its dependence on agronomic management. Our findings highlight that sufficient liming and moderate fertilization can reduce yield variance while maintaining productivity and limiting loss of plant species diversity. Plots receiving the highest rate of nitrogen ferti- lizer or farmyard manure had the highest yield but were also more responsive to environmental variability and had less plant species diversity. We identified the days of water stress from March to October and temperature from July to August as the two main climatic drivers, explaining approximately one-third of the observed yield variance. These drivers helped explain consistent unimodal trends in yield variance—with a peak in approximately 1995, after which variance declined. Here, for the first time, we provide a novel statistical framework and a unique long-term dataset for understanding the trends in yield variance of managed grassland. The application of the criss-cross regression approach in other long-term agro-ecological trials could help identify climatic drivers of production risk and to derive agronomic strategies for improving the climate resilience of cropping systems. Keywords Agronomic management · Biomass production · Climate resilience · Fertilizer input · Food security · Liming · Plant species diversity · Soil pH · Temperature · Water stress * Janna Macholdt 1 Introduction janna.macholdt@landw.uni-halle.de The management of climate-resilient grassland systems is Professorship of Agronomy, Institute of Agriculture and Nutritional Sciences, Martin-Luther-University Halle- important for stable livestock fodder production over time Wittenberg, Betty-Heimann-Strasse 5, 06120 Halle (Saale), (Schmidhuber and Tubiello 2007; Arata et al. 2020; Trnka Germany et al. 2021; Bengtsson et al. 2019; Bommarco et al. 2013; Biostatistics Unit, Institute of Crop Science, University Reckling et al. 2021). However, in the face of climate change of Hohenheim, Fruwirthstrasse 23, 70599 Stuttgart, Germany and the associated increases in abiotic stresses, maintain- Protecting Crops and Environment, Rothamsted Research, ing productivity while minimizing temporal yield variance Harpenden AL5 2JQ, Hertfordshire, UK (or rather improving yield stability) of grassland systems Computational and Analytical Sciences Department, will become increasingly challenging (Olesen and Bindi Rothamsted Research, Harpenden  AL5 2JQ, Hertfordshire, 2002; Ray et al. 2019). Observed climatic changes, such as UK increasing temperatures or weather anomalies, have nega- Section of Environmental Chemistry and Physics, tively affected global grassland productivity (Kipling et al. Department of Plant and Environmental Sciences, University 2016; Höglind et al. 2013; Addy et al. 2022; Brookshire of Copenhagen, Thorvaldsensvej 40, 1871 Copenhagen, Denmark Vol.:(0123456789) 1 3 37 Page 2 of 19 J. Macholdt et al. and Weaver 2015; Wilcox et al. 2017; Hall and Scurlock et al. 1997) because estimates of yield variance based on 1991). For northern Europe as the experimental region of a shorter period might be imprecise and long-term trends this study, Trnka et al. (2021) reported that the annual aver- cannot be detected (Hadasch et al. 2020; Macholdt et al. age temperature as well as the frequency of combined heat 2021; Reckling et al. 2021). Data over decadal time scales and drought stress for grassland has increased in recent also allow the adaptation of grasslands to climatic trends decades. In addition, they showed that by 2050, the area (such as increasing temperatures) to be quantified in addi - of grassland exposed to combined heat and drought may tion to the buffering of short-term climatic variability. double compared with that in 2021 and that northern regions Long-term experiments (LTEs) and their associated data- will experience higher temperatures and drought more fre- sets provide a unique opportunity to examine the effects quently (Trnka et al. 2021). As well as long-term temporal of climate change (Berti et al. 2016). trends, these changes in climatic conditions are expected to Appropriate statistical methods are necessary to handle increase the interannual yield variance of grassland systems the often complex design of LTEs and any experimental and decrease plant species diversity, which could put future modifications that have occurred over time (Reckling et al. food security at risk (Graux et al. 2013; Piseddu et al. 2021). 2021; Payne 2015; Macdonald et al. 2018), such as those The intensification of grassland production through the used for genotype × environment × management interac- use of inorganic fertilizers and simplified swards of fast- tion analyses in other disciplines, including plant breed- growing grass cultivars has delivered increased primary bio- ing (Hadasch et al. 2020) and agronomy (Macholdt et al. mass and live weight gain of livestock (Carswell et al. 2019). 2021). Although there are some recent studies on the yield However, the use of inorganic fertilizer is a major contribu- variance of field crops, detailed knowledge about the long- tor to greenhouse gas emissions (both in manufacturing and term effects of climatic changes and agronomic manage- application) and results in swards with low plant species ment on temporal trends in yield variance for grasslands is diversity that may also compromise resilience (Robertson limited (Dodd et al. 1997). Previous studies reporting yield and Vitousek 2009; Storkey et al. 2015). Recent findings variance of grasslands often cover only short periods with of short-term grassland studies have confirmed that addi- less than 10 years (Tilman et al. 2006; Prieto et al. 2015; tional nutrient inputs often increase not only yield but also Hautier et  al. 2014, 2020; Haughey et  al. 2018; Zhang interannual yield variance and reduce plant species diversity et al. 2016. 2019; Sanderson 2010; Dodd et al. 1997). To (Hautier et al. 2020, 2014; Zhang et al. 2016, 2019; Crawley date, there have been only few analyses of permanent, et al. 2005). Moreover, studies reported that relative climatic managed grassland systems that include temporal trends in adaptability decreased with increasing land use intensity yield variance covering a long-term period (Craven et al. (Deguines et al. 2014). In particular, mineral N (nitrogen) 2018; Isbell et al. 2017; Dodd et al. 1997). These studies fertilizer input reduces the diversity of terrestrial vegeta- often neglect a further requirement, which is that yield tion by favouring fast-growing grass species adapted to high variance should be determined independently of yield nutrient availability (Midolo et al. 2018). There is evidence level, otherwise yield variance can be incorrectly inter- that greater plant species diversity in managed grasslands preted if the time span is too short and there is a systematic may enhance their resilience to climate change and result dependency of variation on the mean (Preissel et al. 2015). in more stable yields in response to disturbance (Tracy and In this study, we used a long-term dataset (1965–2018) Sanderson 2004; Tilman et al. 2006; Haughey et al. 2018; from The Park Grass Experiment at Rothamsted (UK) with Sanderson 2010; Baca Cabrera et al. 2021). However, the consistent long-term management and determined yield relationship between stable productivity (low interannual variance independently of yield level, for most accurate yield variance) and plant species diversity can also strongly statistical estimates. The novelty of our analysis is that we vary depending on agronomic management and soil condi- included environmental abiotic covariates, such as atmos- tions (Bullock et al. 2001; Hector et al. 1999; Tracy and pheric chemistry (wet and dry N deposition, SO, CO ) 2 2 Sanderson 2004; Crawley et al. 2005; Storkey et al. 2015). A and other climatic parameters (air temperature, precipi- very important agronomic management factor for stabilizing tation, soil moisture deficit, etc.), in a novel criss-cross yields over time and for maintaining plant species diversity regression approach (extended Finlay–Wilkinson regres- is liming, particularly in soils prone to acidification (Fornara sion) to determine the main climatic drivers of long-term et al. 2011; Storkey et al. 2015). yield variance and to evaluate the effects of liming and To achieve resilient, sustainable, and stable productiv- fertilizer applications on the relative impact of environ- ity of grasslands, a better understanding of the climatic mental variability on yield sensitivity (or responsiveness) drivers of long-term trends in temporal yield variance and across a range of plots with contrasting productivity and its dependence on agronomic inputs and biological diver- plant species richness. sity is required. Ideally, such assessments should be done We specifically addressed the following three research on long-term datasets (>20 years) (Piepho 1998; Dodd questions: 1 3 Long-term trends in yield variance of temperate managed grassland Page 3 of 19 37 I. Are there temporal trends in interannual yield vari- provided in Fig. A1-b Supplementary material and available ance over the study period, and do these vary in rela- online (https://doi. or g/10. 23637/ ERADOC- 1- 143 ; p. 43). In tion to fertilizer and lime applications? 1856, starting values for PGE soil parameters in the topsoil II. What are the most important environmental abiotic (0–23 cm) were estimated as pH 5.7, 11.6% sand, 66.3% silt, −1 3 drivers explaining yield variance, and does related 22.1% clay, soil weight 2430 t  ha , bulk density 1.1 g  cm , −1 yield sensitivity to these climatic drivers depend on and total soil nitrogen content of 5830 kg N  ha (Lawes agronomic management? and Gilbert 1859). III. Is there a correlation between plant species diversity, The mean annual air temperature and rainfall at the mean yield, and yield variance? Rothamsted site (1981–2010) were 9.8  °C and 733  mm, respectively (Perryman et  al. 2019). In recent decades, the annual mean air temperature (1989–2018) was 1.1 °C 2 Material and methods warmer compared to the previous period (1878–1988: 9.04 °C) (Macdonald et al. 2018). The increasing trend in 2.1 T he Park Grass Experiment (PGE) temperature anomaly at Rothamsted (1880–2018) is pro- vided in Fig.  A2 Supplementary material and available This study was based on the Park Grass Experiment (PGE) online (https:// doi. org/ 10. 23637/ rms- RMAAt empan omaly- at Rothamsted (Fig. 1), initiated by Lawes and Gilbert in 1). The effective plant available water capacity of 135 mm 1856 to examine the effects of different mineral fertilizers in the effective root zone of the Park Grass site implies that and organic manures on the productivity of permanent pas- plant growth is often limited by lack of water in summer ture cut for hay (Silvertown et al. 2006). A detailed descrip- (Avery and Catt 1995). tion of the experiment is available in the Rothamsted Guide The current analysis is focused on the period from 1965 to Classical Experiments (Macdonald et al. 2018) and the to 2018 using plots with constant treatments and a consist- e-RA website (http://www .er a.r othams ted.ac. uk/ e xperiment/ ent harvesting methodology to ensure comparability and rpg5). The Park Grass soil is a moderately well-drained silty accurate estimates of temporal trends in yield variance. The clay loam overlying clay-with-flints (Avery and Catt 1995), design of the PGE is provided in Fig. 1 (more details are a Chromic or Vertic Luvisol according to the FAO clas- provided in Fig. A1-a Supplementary material), consisting sification. The experimental site shows a relatively uniform of 24 main plots with contrasting fertilizer treatments, each soil, based on comprehensive soil analyses made by Avery divided into four sub-plots ‘a, b, c, and d’ for liming treat- and Catt (1995)—the ‘Soils at Rothamsted Colour Map’ is ments. Sub-plots a, b, and c receive lime, if needed, every 3 Fig. 1 Park Grass Experi- ment aerial view (left) and plot layout (right). Loca- tion: Harpenden, UK, Herts, AL5 2JQ (51°48′12.33″N; 0°22′21.66″W; 130 m a.s.l.). Detailed information about plot layout and treatments are shown in Tables A1 and A3 Sup- plementary material. Source: electronic Rothamsted Archive (http:// www. era. rotha msted. ac. uk/ Park# images/). 1 3 37 Page 4 of 19 J. Macholdt et al. years to maintain soil pH 7, 6, and 5, respectively; sub-plot d (Crawley et al. 2005; Storkey et al. 2015; Ray et al. 2015). receives no chalk. The yield data selected for this study were Although plant communities respond to these drivers, the taken from the a, b, c, and d sub-plots of seven contrasting relative differences between the plots in species richness main plots (red marked in Fig. A1 Supplementary material) and diversity are conserved in time (dynamic equilibrium); chosen to represent a range of productivity and species rich- this meant that the relationship between plant diversity and ness (plot3: Nil/unfertilized, plot7/2: ‘P K Na Mg’, plot6: yield variance could be included in our analysis (Silvertown −1 ‘N1 P K Na Mg’ with 48 kg N  ha ammonium sulfate, et al. 2006). −1 plot 9/2: ‘N2 P K Na Mg’ with 96 kg N  ha ammonium −1 sulfate, plot11/1: ‘N3 P K Na Mg Si’ with 144 kg N  ha 2.2 Environmental abiotic covariates ammonium sulfate, plot13/2: manure applied in a 4-year −1 cycle, plot17: N*1 with 48 kg N  ha sodium nitrate), except Based on statistical analyses (see Section  2.3), the effect on plot 6, where only the a and b sub-plots were available. of the following environmental abiotic covariates on tem- The plots 7/2, 6, 9/2, 11/1, and 17 represent different lev - poral yield variance was tested. Different climatic covari - els and forms of inorganic fertilization, whereas the plot ables (air temperature, humidity, rainfall, hours of sun- 13/2 represents organic fertilization (poultry manure since shine, soil moisture, radiation, wind) were measured daily 2003, before farmyard manure). Additional treatment details at Rothamsted Research (available at: http:// www. era. rotha are provided in Table A3 Supplementary material. Yield msted. ac. uk/# measu remen ts), except for measurements of data were obtained from the electronic Rothamsted Archive atmospheric carbon dioxide concentration (Fig. A4 Sup- ‘e-RA’ and have been made publicly available on the e-RA plementary material) that was recorded monthly. In addi- website ( https:// doi. org/ 10. 23637/ r pg5- yield s1965- 2018- tion, air chemistry parameters were measured, N deposition, 01) (Perryman and Ostler 2021) and included plot-specific and SO emissions (1965–2018, UK), which are shown in twice-yearly yields (total aboveground biomass; sum of 1st Fig. A5 Supplementary material. Besides these measured and 2nd cuts; 100% dry matter). The first cut was made into covariates, the accumulated numbers of water stress days for hay and removed in mid-June, and the second cut was taken the vegetation periods 1st March–15th of June (time period with a forage harvester while still green (end-October, bio- up to typical date for 1st cut/harvest), 16th June–31st Octo- mass removed). ber (time period up to typical date for 2nd cut/harvest), and The plant communities on the PGE are naturally assem- for the entire period (Fig. 2) were calculated yearly in the bled from a regional species pool that is classified as dicot- following manner. The amount of maximum plant available yledon-rich Cynosurus cristatus–Centaurea nigra grassland, water (MPAW [mm]), which acted as a ‘bucket of water’, one of the mesotrophic grassland communities in the British was specified for the soil. At the start of the calculation (1 National Vegetation Classification system (Dodd et al. 1994). Jan 1965), it is assumed that the actual plant available water The botanical composition of all sub-plots was studied annu- (APAW(t)) equaled the MPAW (‘the bucket was full’). The ally from 1991 to 2000 and from 2010 to 2012 by recording daily water surplus was calculated as rainfall minus potential the dry mass of each plant species in early June (number evaporation over grass (WS(t)), with potential evaporation of species, Shannon’s diversity index) (Storkey et al. 2015; over grass as a derived meteorological variable (see formula Crawley et  al. 2005). The species of grasses, forbs, and details: http:// www . er a. r o t ha ms ted. ac. uk/ info/ me t/ der iv legumes comprising at least 5% of the aboveground bio-ed_ varia bles# EVAPG). A daily balance was calculated as mass found in these surveys and used for this study have DB(t) = APAW(t − 1) + WS(t). For positive WS(t)-values, been made publicly available on the e-RA website (https:// APAW increased to the maximum of MPAW. Above this doi. org/ 10. 23637/ rpg5- speci es_ 1991- 2000- 01) (Perryman value, water was assumed to run off or drain away. For nega - et al. 2021). An overview about the changes in the number tive values, APAW was reduced, and when reaching 0, no of plant species over time is provided in Fig. A11 Supple- more water extraction was possible. mentary material; a more detailed description of the PGE biodiversity data has been reported in the Rothamsted Guide If DB(t) > MPAW; APAW(t) = MPAW to Classical Experiments on pp. 25–27 (Macdonald et al. If MPAW > DB(t) > 0; APAW(t) = DB(t) 2018) and is available online: http:// www. era. rotha msted. If DB(t) < 0; APAW(t) = 0 ac.uk/ home/ W eb_L TE_Guide book_ 2018_ 2019- r eprint. pdf . Because of the resource required for vegetation assessments, Days where APAW was 0 were counted as stress days. data on plant diversity are only available for a sub-set of the The numbers of water stress days accumulated for the 1st years and plots used in this analysis. However, the data cover March–15th of June (harvest time of the 1st cut), 16th a significant proportion of the time period analysed in our June–31st of October (harvest time of the 2nd cut), and for study that includes large interannual variation in weather the entire period each year were determined. No water stress and yields and a period of change in atmospheric chemistry occurred before March or after October. Several values of 1 3 Long-term trends in yield variance of temperate managed grassland Page 5 of 19 37 Fig. 2 Accumulated number of water stress days for the vegeta- tion period from March to Octo- ber and temporal development of the mean air temperature (°C) for the months July–August at Rothamsted (1965–2018). Water stress was defined as a limited plant available soil water content (formula described in “Material and methods” sec- tion). Further information about the temperature anomaly is provided in Fig. A2 Supplemen- tary material. MPAW were tested, but an MPAW of 135 mm provided the spaced knots (set to 10, approximately 5-year intervals). For best correlation for both the 2nd cut and the yearly yield, in each of these 10 sub-periods, a separate residual variance line with the soil description (Kohler et al. 2016). was fitted to assess changes in temporal yield variance, with lower values indicating less unexplained variability between 2.3 Statistical analysis years (=more stable yields). Each ‘fertilizer × liming’ treat- ment combination was assumed to have a specific set of To account for the experimental design, a mixed model variance components ([t/ha] ) for the different periods. We was used based on REML, as recommended by Raman denote the period-specific yield variances as ‘environmental et al. (Raman et al. 2011) and Onofri et al. (Onofri et al. variance’, a term coined by Römer (1917). Here, we replace 2016). Each plot had a different ‘fertilizer × liming’ treat- the arithmetic treatment mean with the spline estimate of ment combination, with no replication or randomization the temporal trend. The model, therefore, allowed yield vari- (Fig. A1 Supplementary material), so that plot errors and ance to be determined independently of yield level, which ‘fertilizer × liming’ interactions (=residual) could not be highlights that the statistical approach used here differs from separated. The size of the main plots (306–1912  m ) might earlier approaches based on the PGE, which used the clas- partly compensate for the lack of replication, particularly sic coefficient of variation for a measure for yield variabil- because the experimental site was reasonably uniform when ity (Dodd et al. 1997). This is an important improvement the experiment started in 1856 (Crawley et al. 2005; Lawes because yield variance can be incorrectly interpreted if there and Gilbert 1859). is a systematic dependency of the measure of variation on To answer the first research question comparing tem- the mean. No such dependencies were found in this analysis. poral trends in mean yield and yield variance across fer- The model was fitted separately for each ‘fertilizer × liming’ tilizer × liming treatments, we fitted for each plot (‘ferti- treatment. In this analysis, any systematic year effects (trend, lizer × liming’ treatment combination) a smoothing spline etc.) are captured by the fitted splines, whereas the random for the mean yield trend via a random-effects specification residual captures any unexplained year effects. as a mixed model (Verbyla et al. 1999). An advantage of the To answer the second research question of how much of smoothing spline approach over other regression approaches the variation in yield can be explained by climatic drivers is that no specific assumption is needed as regards the func- and the interaction with agronomic management, treatment- tional form of the trend. Preliminary inspection of the data specific multiple regression analyses were performed. The revealed that such flexibility was needed. The entire model strength of the relationship between the response vari- syntax (using R version 4.0.0) is provided in Table A6 Sup- able ‘yield’ and several explanatory environmental abiotic plementary material. The trends were fitted using ASReml- covariates (as described in former Chapter 2.2), as well as R using the specification random = ~spl(t ,k = 10), where t the importance of each of the predictors to the relationship, is continuous time in years and k is the number of evenly was assessed. We further tested for significant differences 1 3 37 Page 6 of 19 J. Macholdt et al. in regression coefficients of the main effects between treat- covariates. With regard to the nonnormal distributions of ments by using a general linear model approach, but no sig- the underlying data (outliers), correlation coefficients were nificant differences between treatments were found. Further - calculated from the ranks of the data, not from their actual more, a novel criss-cross regression approach with extended values. Kendall’s tau (τ) was used to ensure that the results Eberhart–Russell regression analyses (Eberhart and Russell were accurate because the same ranks were repeated too 1966) was used, in which the environmental mean was mod- many times in the partial datasets (e.g. botanical surveys). eled using the previously as significant selected environ - The strength of the correlation increased both from 0 to +1 mental abiotic covariates ((1) accumulated days of water and 0 to −1, where −1/+1 indicated the strongest correla- stress from March to October; (2) mean air temperature tion and 0 indicated no correlation. The sign of τ showed the from May to June; (3) mean air temperature from July to direction of the correlation; if negative, the variables were August based on the multiple regression analyses; shown inversely related. To gain further insight into the temporal in Table  1). This new method allowed us to assess the dynamics of the plant communities, a multivariate analy- treatment-specific yield sensitivity (or responsiveness) to sis was done using the data on relative biomass recorded variability in these climatic conditions. The Eberhart–Rus- on all study plots for the 10-year period (1991–2000) for sell model was further fitted using the criss-cross regres- which data were available. Although this only covered part sion approach proposed for Finlay–Wilkinson regression by of the total period covered by the yield variance analysis, Digby (2009) and extended for the mixed model version it included a large range of yields. Two partial canonical by Nabugoomu et al. (1999). This iterative scheme allowed correspondence analyses (pCCA) were done after remov- for (i) regression of the environmental index on covariates, ing species that were recorded in less than 5% of samples to (ii) heterogeneity of variance for the independent devia- avoid bias owing to rare species. First, the effect of year was tions from the regression, and (iii) serial correlation of the analysed, including plot as a covariate. Second, the variance residuals. The intercept of Finlay–Wilkinson regressions in community composition explained by plot was quanti- provided information about the general yield level of a treat- fied, including year as a factorial covariate. In both cases, ment compared to others. The slopes of regression lines can the proportion of functional group (grass/forb/legume) was be interpreted as ‘yield sensitivity’ to climatic perturbation included in the ordination plots as a supplementary variable. (slope > 1: higher sensitivity; slope < 1: less sensitivity/ better resilience). Treatments with a slope of approximately 1 showed an average yield reaction across the treatments, 3 Results and discussion similar to the average response indicated by the environ- mental mean (reference, black regression line with a slope 3.1 Temporal trends in yield variance of 1; Fig. 5a–d). We would like to stress that this criss-cross and identifying main climatic drivers regression approach for the extended Finlay–Wilkinson regression was specifically developed for the PGE and, to the Based on the unique long-term dataset and the new designed best of our knowledge, constitutes a novel method. Briefly, statistical approach, this study provides, for the first-time, the Finlay–Wilkinson regression model can be written as insights into trends in temporal yield variance of grassland =  +  w , where  is the expected performance of the in response to climatic drivers and its dependence on agro- ij i i j ij ith treatment in the jth environment,  and  are intercept nomic management; yield variance was determined indepen- i i and slope (‘yield sensitivity’) for the ith genotype, and w is dently of yield level, which is important to avoid any misin- the environmental mean of the jth environment. The envi- terpretation of resilience. Our work adds a new perspective ronmental mean, in turn, is modelled by a linear regression to earlier productivity analyses of grassland experiments, as w =  +  x +  x + .... +  x , where x is the value which were based mostly on shorter time periods (Prieto j 0 1 1j 2 2j p pj hj of the hth covariate in the jth environment and  ,… ,  are et al. 2015; Tilman et al. 2001; Graux et al. 2013; Trnka 0 p regression parameters. A detailed description of the model et al. 2021; Storkey et al. 2015; Haughey et al. 2018; Sander- and its estimation is provided in Table A7 Supplementary son 2010), and reveals the importance of long-term experi- material. ments for detecting possible trends over time. It is comple- Regarding the third research question, we aimed to iden- mentary to a recent analysis of trends in productivity on the tify the potential correlation between plant species diversity PGE that showed a consistent decline in yields in response to (number of species and Shannon’s diversity index) and the climate change across a longer time period (1902–2016) for mean yield or temporal yield variance depending on the four treatments and using yield data from just the first hay agronomic treatment (fertilizer × liming). For this correla- cut (Addy et al. 2022). By focusing on yield variance (using tion analysis, we used the mean yield and temporal yield data from both cuts) across a wider range of treatments, our variance data based on the statistical analyses as described study provides additional insights into the potential for man- above to account for the experimental design and other agement to impact adaptability and resilience of grasslands 1 3 Long-term trends in yield variance of temperate managed grassland Page 7 of 19 37 1 3 Table 1 Treatment-specific multiple regression analyses for the dependent variable ‘yield’ and selected climatic drivers. Analyses based on data (1965–2018) of total yield (sum of 1st + 2nd cut). Abbreviations: VIF = variance influence factor (measures the strength of the correlation between the independent variables in regression analysis); regression coefficient B = slope ‘b’ of regression line (how much Y changes for each change in X); standardized coefficient Beta = analogous to the interpretation of ‘b’, except that Beta expresses change in standard scores. Climatic drivers with significant main effects were selected based on a preceding general linear model analysis: accumulated number of water stress days from March to October (limited plant available soil water content with expected drought stress) and mean air temperature (=Temp.) in summer months (Fig. 2). No significant differences in regression coefficients B between treatments were found. Treatment explanations are provided in Table A3 Supplementary material (FYM/PM = farmyard/poultry manure; N* = sodium nitrate). Treatment Fertilization Liming General model statistic Second selected covariable First selected covariable Durbin Adjusted VIF Covariable Regression Standardized p value Covariable Regression Standardized p value Watson R Square Collinearity coefficient coefficient coefficient coefficient Beta test B Beta B 3 Nil pH 7 2.02 0.29 1.00 Water stress -0.02 -0.55 0.00 none selected n/a n/a n/a pH 6 1.68 0.23 1.00 Water stress -0.02 -0.49 0.00 none selected n/a n/a n/a pH 5 1.62 0.29 1.00 Water stress -0.02 -0.54 0.00 none selected n/a n/a n/a no chalk 1.58 0.25 1.00 Water stress -0.02 -0.50 0.00 none selected n/a n/a n/a 7/2 P K Na Mg pH 7 1.97 0.49 1.36 Water stress -0.03 -0.47 0.00 Temp.July-Aug. -0.66 -0.35 0 00 pH 6 2.05 0.37 1.36 Water stress -0.02 -0.37 0.01 Temp.July-Aug. -0.63 -0.35 0 01 pH 5 1.76 0.34 1.01 Water stress -0.04 -0.52 0.00 Temp.Mar.-Apr. 0.96 0.35 0 00 no chalk 1.77 0.25 1.00 Water stress -0.03 -0.52 0.00 none selected n/a n/a n/a 6 N1 P K Na Mg pH 7 2.05 0.48 1.29 Water stress -0.02 -0.42 0.00 Temp.July-Aug. -0.76 -0.40 0 00 pH 6 2.06 0.40 1.29 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.70 -0.37 0.01 pH 5 n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a no chalk n/a n/a n/a  n/a n/a n/a n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 1.71 0.36 1.36 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.79 -0.31 0.02 pH 6 1.51 0.36 1.36 Water stress -0.02 -0.39 0.01 Temp.July-Aug. -0.64 -0.38 0.01 pH 5 1.70 0.35 1.36 Water stress -0.02 -0.36 0.01 Temp.July-Aug. -0.74 -0.35 0.01 no chalk 1.57 0.38 1.36 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.85 -0.32 0.02 11/1 N3 P K Na Mg pH 7 1.64 0.29 136 Water stress -0.02 -0.32 0.02 Temp.July-Aug. -0.73 -0.33 0.02 pH 6 1.62 0.27 136 Water stress -0.02 -0.29 0.04 Temp.July-Aug. -0.67 -0.34 0.02 pH 5 1.71 0.26 136 Temp.July-Aug. -1.30 -0.52 0.00 Waterstress -0.02 -0.25 0.05 no chalk 1.82 0.35 136 Temp.July-Aug. -1.47 -0.60 0.00 Waterstress -0.01 -0.16 0.05 13/2 FYM/PM pH 7 1.72 0.34 136 Temp.July-Aug. -0.91 -0.40 0.00 Waterstress -0.02 -0.29 0.03 pH 6 2.16 0.38 136 Temp.July-Aug. -0.91 -0.44 0.00 Waterstress -0.02 -0.29 0.03 pH 5 1.72 0.33 136 Water stress -0.03 -0.39 0.01 Temp.July-Aug. -0.67 -0.30 0.03 no chalk 1.64 0.21 1.00 Water stress -0.03 -0.47 0.00 none selected n/a n/a n/a 37 Page 8 of 19 J. Macholdt et al. to abiotic stress. Our treatment-specific model allowed tem- poral trends, using a spline component, to be dissected from residual variance, fitted as separate parameters for ten con- secutive sub-periods. As a striking finding, we identified a similar pattern of temporal trends in yield variance (overview results Fig. 3) (Römer 1917), particularly for plots with low rates of N −1 (48 kg N  ha , N1) or no N fertilizer ‘Nil’ (plot specific results Fig. 4). This pattern can be described as follows: in the 1960s and 1970s, the plots showed relatively low yield variance. In subsequent years, we observed greater yield variance with a peak at approximately 1995 before becom- ing more stable in the past two decades (Fig. 3). This pattern of temporal trends in yield variance appeared to be less pro- nounced in treatments with a high soil pH; in contrast, peaks were more pronounced in treatments with low pH (Fig. A8 Supplementary material). The exceptions were treatments with reduced or no liming, but with a high mineral N supply (Fig. 4e) or farmyard manure (FYM) (Fig. 4g), where these patterns were not as evident, and periods of high yield vari- ability were observed throughout the study period. Key message for ‘Research question I’: There were con- sistent unimodal trends observed in yield variance in plots with low to moderate or no nitrogen fertilizer additions, with a peak in 1990s, after which variability declined. Yield was most variable in plots with higher nutrient inputs and lower soil pH. Based on treatment-specific multiple regression analy - ses including environmental abiotic covariates, we explored the possible causes of the observed temporal trends in yield variance. We identified the accumulated days of water stress from March to October and the mean air temperature from July to August as the two main climatic drivers, explain- ing around one-third of the observed yield variance across treatments, or even up to 49%/48% of yield variance in treat- ment P K Na Mg and N1 P K Na Mg, respectively (Table 1 and Fig. 2). Overall, the impact of temperature driving yield variance was lower than the stronger impact of water stress (Table 1: see columns ‘standardized coefficient Beta’ with lower absolute values indicating less impact). The identi- fied main climatic drivers tally with previous analyses of the PGE made by Dodd et al. (1994, 1997) and recent ones by Addy et al. (2021, 2022); any differences can be explained by the fact that we included the second cut in the yield data. The seasonality of temperate grassland production is primar- ily affected by soil moisture and temperature, which con- strain the length and determine the intensity of the growing season (Trnka et al. 2011). Generally, significant tempera- ture changes, particularly hot temperatures in summer with a limiting water balance, have a negative effect on grassland productivity (Kipling et al. 2016; Höglind et al. 2013) and are expected to increase the interannual and seasonal pro- duction variability of grassland systems (Graux et al. 2013; 1 3 Table 1 (continued) Treatment Fertilization Liming General model statistic Second selected covariable First selected covariable Durbin Adjusted VIF Covariable Regression Standardized p value Covariable Regression Standardized p value Watson R Square Collinearity coefficient coefficient coefficient coefficient Beta test B Beta B 17 N*1 pH 7 1.95 0.31 1.19 Water stress -0.02 -0A2 0.00 Temp.May-June -0.25 -0.27 0.04 pH 6 1.91 0.23 136 Water stress -0.01 -0.29 0.04 Temp.May-June -0.24 -0.29 0.05 pH 5 1.80 0.32 136 Temp.May-June -0.25 -0.58 0.00 Waterstress -0.02 -0.34 0.01 no chalk 1.73 0.28 1.00 Temp.July-Aug. -0.58 -0.46 0.00 Waterstress -0.01 -0.22 0.05 Long-term trends in yield variance of temperate managed grassland Page 9 of 19 37 Fig. 3 Summary plot for the overview temporal trend in mean yield (blue dotted line; grey crosses real harvest data) and yield variance (red bars) including standard errors (grey error bars) based on the total mean over all liming × fertiliza- tion treatments (1965–2018). Underlying plot specific results are provided in Fig. 4; overview of liming treatments (Fig. A8 Supplementary material), and overview of fertilization treat- ments (Fig. A9 Supplementary material). Detailed informa- tion about treatments is shown in Table A3 Supplementary material. Chang et al. 2017). However, the water stress index we cal- this study). A recent study by Addy et al. (2021) identified culated is partly based on a derived meteorological (potential clusters of years with similar weather patterns between 1900 evaporation from grass) and an estimated parameter (soil and 2020 at Rothamsted, which might help to clarify the water storage), which may differ from actual measurements. unexplained rest of unimodal yield variance. They found a This may have contributed to the lower correlations seen climate cluster characterized by cool and dry springs from in the Nil (no fertilizer) and N*1 plots together with the approximately the 1960 to 1970s (in our study: stable yields fact that other nutrients may have limited growth. Further- were observed at the beginning), followed by a period with more, water stress often occurs in July–August, so the fact a variety of clusters and widely varying weather patterns that both factors come out as important may indicate that (in our study: more variable yields occurred from approxi- plants may be hit harder by high temperatures when transpi- mately 1980 to 2000), and a transition since 2000 with an ration and thereby cooling are limited. The effect of water increased tendency toward higher temperatures in springs stress was consistently present in all treatments and can be and drier periods in June (in our study: associated with more assumed to be the main driving factor in the PGE. Regard- stable yields at the end; Fig.  3). The study suggests that ing the effect of agronomic practices on resistance to water the positive effects of sufficient water availability can off- stress (Table 1 and Fig. 2), sufficient liming is important for set the negative effects of warmer temperatures on pasture sustaining grassland productivity due to positive effects on performance (Addy et  al. 2021). An additional explana- soil pH, root growth, soil structure, nutrient availability, soil tion could be that this peak in yield variance between 1980 carbon, and soil biota (Holland et al. 2018; Fornara et al. and 2000 coincides with a period of decreasing nitrogen 2011). deposition and SO emissions (see Fig. A5 Supplementary Key message for ‘Research question II—part A’: The material), which is reflected in shifts in plant community accumulated days of water stress from March to October composition, especially an increase in the relative propor- and mean air temperature from July to August were the most tion of legumes. Regarding the SO emissions, there was a −1 −1 important climatic drivers, explaining approximately one- decline from approximately 65 kg  ha in 1980 to 5 kg  ha third of the observed interannual yield variance. in 2006 (Anon 2006). We might expect this to affect species The remaining unexplained yield variance in the PGE diversity only in S-limited plots, but S is applied together might be driven by different environmental factors, which with K, Na, and Mg in most plots (such as in the ‘P K Na could not be determined further (i.e. due to lack of data in Mg’ treatment) in the PGE, so any effect of changes in air 1 3 37 Page 10 of 19 J. Macholdt et al. Trend in mean yield Temporal yield variance [(t/ha)²] [t/ha] 15 12 A) Nil Trend in mean yield Yield variance 10 8 5 5 4 0 0 15 12 B) P K Na Mg - pH 7 10 8 5 4 0 0 15 12 C) N1 P K Na Mg 5 4 0 0 15 12 D) N2 P K Na Mg 10 8 5 4 0 0 15 12 E) N3 P K Na Mg 10 8 5 4 0 0 15 12 F) N*1 10 8 5 4 0 0 15 12 G) FYM/PM 10 8 5 4 0 0 1965 1985 2005 1965 1985 2005 1965 1985 2005 1965 1985 2005 pH = 7 pH = 6 pH = 5 no chalk 1 3 Long-term trends in yield variance of temperate managed grassland Page 11 of 19 37 ◂Fig. 4 Temporal trends in mean yield (blue splines with approxi- and transpiration, as affected by increasing temperature mated confidence intervals) and temporal yield variance (red bars (Fig. A2 Supplementary material) and atmospheric C O con- incl. standard errors) depending on the treatment (fertilizer  ×  lim- centrations (C ) (Fig. A4 Supplementary material). It was ing) for the experimental period of 1965–2018. Analysis based on noted that grasses appeared to be more sensitive to increas- year  ×  plot specific yields; these raw yield data are shown as black dots in each graphic. A Treatment no. 3: Nil (no fertilizer input)—pH ing C than forbs, resulting in lower water use efficiency, 7/6/5/no chalk. B Treatment no. 7/2: P K Na Mg—pH 7/6/5/no chalk. decreased N uptake, and declining biomass production in C Treatment no. 6: N1 P K Na Mg—pH 7/6 (restricted data availabil- grass-rich communities (Baca Cabrera et al. 2021). This indi- ity: only 1972–2018 and pH 7/6). D Treatment no. 9/2: N2 P K Na cated that plant communities comprising only a few species Mg—pH 7/6/5/no chalk. E Treatment no. 11/1: N3 P K Na Mg—pH 7/6/5/no chalk. F Treatment no. 17: N*1 (N* = sodium nitrate)—pH may be disproportionately affected if the species present are 7/6/5/no chalk. G Treatment no. 13/2: FYM/PM (farmyard/poultry poorly adapted to changing climatic conditions and are also manure)—pH 7/6/5/no chalk. Yield variance denoted Römer’s envi- less adaptable to long-term climatic trends (Eisenhauer et al. ronmental variance, with lower values indicating more stable yields 2019). Our results on trends in yield variance and mean yield and higher values indicating more variable yields. Detailed informa- tion about treatments is shown in Table A3 Supplementary material. support this hypothesis (Table 2 and Fig. 4). Summary plots are provided as an overall overview (Fig. 3), overview of liming treatments (Fig. A8 Supplementary material), and overview 3.2 Determine yield sensitivity to climatic changes of fertilization treatments (Fig. A9 Supplementary material). using a novel criss‑cross regression analysis chemistry is more likely a result of decreasing N deposition We present a novel criss-cross regression analysis (extended in these plots. The highest peak in yield variance in the mid- Finlay–Wilkinson regression), in which we modeled the dle period was observed in plots lacking N but supplied with environmental mean using the three main identified cli - P, K, Na, and Mg (Fig. A9). These plots have the highest matic factors (Table 3 and Fig. 5). This approach allowed proportion of legumes (Table 2) that are sensitive to changes us to determine the treatment-specific yield sensitivity to in N deposition. This suggests that where yield is dependent variation in climatic conditions. A great advantage of this on biological N fixation (see key message IV), resilience approach is its parsimony, providing a single slope for each will also be determined by the specific response of legumes treatment that assesses sensitivity to all included covariates to environment change—such as atmospheric N deposition simultaneously. The interaction between fitted environmen - (Storkey et al. 2015). tal mean and treatments was highly significant (F  = 3.17; In addition to yield variance, the temporal trends in mean P < 0.0001; Table A7 Supplementary material), showing yield were in some cases relatively stable, as for the unfer- that the slopes are different between treatments. The results tilized ‘Nil’ treatment with a relatively constant yield of shown in Table 3 can be interpreted as follows: a lower slope −1 approximately 3 t  ha 286 (Fig.  4a). In other cases, the indicates less sensitivity (or responsiveness) of yield to vari- mean yield decreased over time (negative trends), supporting ation in climatic conditions; a lower variance of deviation the conclusions of Addy et al. (2022), who analysed PGE indicates a more stable yield. A higher slope indicates a data over a longer time period (1902–2016) and modelled higher sensitivity of yield to changing climatic conditions; spring hay yields under four fertilizer regimes in response a higher variance of deviation indicates more yield fluc- to seasonal temperature and rainfall. The PGE modelling tuations. Regarding the predicted environmental mean, the study showed that warmer and drier years in the twentieth signs of all three regression coefficients (theta1–theta3) and twenty-first centuries resulted in yield reductions and were negative, meaning the predicted environmental mean are forecasted to decline further up to 50% under future increases with decreasing values for x1, x2, and x3. Thus, (2020–2080) climate scenarios (Addy et al. 2022). less water stress (x1) and lower temperatures from May to This was the case particularly for treatments with greater June (x2) and from July to August (x3) resulted in higher inputs of N and those provided with organic manures, includ- grassland yields during the observed experimental period. ing the ‘N3 P K Na Mg—pH 7 or 6’ (Fig. 4e) and the ‘FYM/ Key message for ‘Research question II—part B’: Liming PM—pH 7 or 6’ (Fig. 4g) treatments in which mean yields to attain a soil pH of 6–7 and moderate N supply (e.g. treat- −1 −1 decreased from >9 t  ha to nearly 5 t  ha from 1965 to ment ‘N2 P K Na Mg—pH 6/7) were identified as the most 2018. In the PGE, the higher input plots (see ‘N3 P K Na promising agronomic practices for sustaining yield under Mg—pH 5 or no chalk’) were dominated by only a few grass varying climatic conditions and for reducing yield sensitivity species (Table 2), which might have resulted in greater yield to abiotic stresses. variance when the climatic conditions were unfavourable for In the intensively fertilized treatment without liming ‘N3 these species and reduced adaptability. Baca Cabrera et al. P K Na Mg—no chalk’, the regression lines had a slope >1 (2021) reported that the declining yields observed in grass- (b = 1.32), which indicates a high yield sensitivity to abiotic rich plant communities of the PGE over the last century were stress—where water is limiting or there is heat stress in the associated with decreases in N uptake, stomatal conductance, summer, proportionally more of the potential yield is lost 1 3 37 Page 12 of 19 J. Macholdt et al. Table 2 Overview results for treatment-specific mean yield, tempo- 2010–2012; species number: 1974, 1991–2000, 2010–2012; propor- ral yield variance, and plant species diversity.  n/a: no data available. tion legumes: 1991–2000). Treatment explanations are provided in Mean values based on available years of surveys (soil pH: 1998– Table  A3 Supplementary material (FYM/PM  =  farmyard/poultry 2014; yield: 1965–2018; Shannon’s diversity index: 1991–2000, manure; N* = sodium nitrate). Treatment Fertilization Liming Soil pH Mean yield [t/ha] Yield variance [(t/ha) ] Shannon's Species Proportion diversity number legumes index [%] Parameter Standard error Parameter Standard error estimate estimate error 3 Nil pH 7 7.2 3.14 0.16 1.11 0.21 2.6 31 5.8 pH 6 6.3 3.52 0.18 1.51 0.35 2.7 30 3.6 pH 5 5.1 2.12 0.23 0.93 0.26 2.1 27 1.1 no chalk 5.2 2.67 0.29 1.40 0.12 2.2 26 0.3 7/2 P K Na Mg pH 7 7.0 7.02 0.26 2.10 0.17 2.4 23 28.6 pH 6 6.2 7.46 0.16 1.97 0.15 2.3 23 16.5 pH 5 5.2 5.94 0.35 2.59 0.19 2.3 22 25.7 no chalk 5.0 4.15 0.32 2.63 0.16 2.1 23 17.4 6 N1 P K Na Mg pH 7 7.0 7.16 0.24 1.94 0.41 2.5 22 17.1 pH 6 5.9 7.06 0.24 2.01 0.42 2.0 20 11.3 pH 5 n/a n/a n/a n/a n/a n/a n/a n/a no chalk n/a n/a n/a n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 7.1 7.32 0.26 1.62 0.21 2.1 17 2.7 pH 6 6.2 7.37 0.26 1.82 0.19 1.8 17 1.6 pH 5 5.0 6.45 0.22 2.69 0.22 1.5 15 10.8 no chalk 3.6 5.39 0.35 2.38 0.29 0.6 3 0.0 11/1 N3 P K Na Mg pH 7 7.0 8.45 0.38 2.18 0.31 1.5 11 0.0 pH 6 6.1 7.54 0.30 1.99 0.19 1.6 12 0.0 pH 5 5.1 7.16 0.36 2.75 0.27 0.8 10 0.0 no chalk 3.5 6.54 0.27 3.38 0.41 0.0 2 0.0 13/2 FYM/PM pH 7 6.9 7.07 0.21 1.87 0.18 2.3 20 0.6 pH 6 6.0 7.75 0.19 2.27 0.23 2.3 20 4.6 pH 5 5.3 6.86 0.35 2.79 0.39 2.3 20 0.2 no chalk 5.1 6.17 0.40 2.85 0.27 2.0 19 0.5 17 N*1 pH 7 7.1 3.95 0.22 1.07 0.25 2.2 23 1.6 pH 6 6.3 4.03 0.18 1.25 0.22 2.1 24 0.1 pH 5 5.8 4.11 0.14 1.24 0.17 2.0 22 0.0 no chalk 5.8 3.92 0.12 0.82 0.11 2.2 24 0.0 (Table 3 and Fig. 5d). In comparison, for treatment ‘N2 P K In addition to slopes, the intercept of regression lines is the Na Mg—pH 6/7’, lower slopes (b = 0.95/0.99) of regression second important criterion of Finlay–Wilkinson regression lines were found, suggesting less yield sensitivity (or respon- analysis. The intercept provides information about the overall siveness) to climatic perturbation (Table 3 and Fig. 5a, b). yield level of a treatment (Table 3 and Fig. 5), with a higher This confirms findings of short-term grassland studies show - intercept indicating higher yield performance (e.g. treatment ing that higher/intensive nutrient inputs often increase yield ‘N3 P K Na Mg—pH 6/7’) and a smaller (or more negative) but can destabilize productivity (Hautier et al. 2014, 2020; value referring to treatments with low yield level (e.g. unfer- Zhang et al. 2016, 2019; Crawley et al. 2005). Moreover, the tilized treatment ‘Nil—no chalk’). Overall, the combined results showed that plots with greater fertilizer inputs were evaluation of slope, intercept, and variance of deviation can less resistant to climatic perturbation (higher yield variance), point to treatments that show a favourable combination of which is in agreement with the findings of Deguines et al. resistance to climatic perturbation (slope < 1) together with (2014), who reported that relative adaptability decreased a high and stable yield level (high intercept, low variance of with increasing land use intensity. deviation). In this study, liming (pH 7/6) and moderate N 1 3 Long-term trends in yield variance of temperate managed grassland Page 13 of 19 37 Table 3 Treatment-specific criss-cross regression analyses (see predicted mean yields for each year (=environmental mean, x-axis) Fig.  5; detailed explanation Table  A7 Supplementary material).  The and regressed on the fertilization  ×  liming treatment yields (y-axis); calculation of the environmental mean considered the three selected getting regression lines with treatment-specific intercepts and slopes. climatic drivers (x1 ‘accumulated number of water stress days’; x2 Treatment explanations are provided in Table  A3  Supplementary ‘mean air temperature from May to June’; x3 ‘mean air temperature material (FYM/PM = farmyard/poultry manure; N* = sodium nitrate). from July to August’ (Fig.  2 and Table  1), and were used to obtain Treatment Fertilization Liming Criss-cross regression analyses Slope Intercept Variance of deviation Parameter Standard error Parameter Standard error Parameter estimate estimate estimate 3 Nil pH 7 0.70 0.14 -1.21 0.96 0.04 pH 6 0.76 0.14 -1.18 0.95 0.03 pH 5 0.69 0.15 -1.93 1.01 0.12 no chalk 0.88 0.16 -2.74 1.06 0.20 7/2 P K Na Mg pH 7 1.18 0.17 -0.21 1.11 0.29 pH 6 1.99 0.17 1.14 1.13 0.33 pH 5 1.17 0.18 -1.16 1.16 0.39 no chalk 1.02 0.17 -1.60 1.14 0.36 6 N1 P K Na Mg pH 7 1.18 0.17 -0.11 1.12 0.30 pH 6 1.07 0.18 0.38 1.14 0.33 pH 5 n/a n/a n/a n/a n/a no chalk n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 0.93 0.17 1.86 1.14 0.34 pH 6 0.98 0.18 1.55 1.17 0.42 pH 5 1.21 0.17 -1.17 1.11 0.29 no chalk 1.19 0.19 -1.90 1.27 0.62 11/1 N3 P K Na Mg pH 7 1.02 0.21 2.21 1.35 0.83 pH 6 0.88 0.17 2.41 1.14 0.36 pH 5 1.09 0.20 0.76 1.29 0.69 no chalk 1.36 0.25 -1.79 1.60 1.50 13/2 FYM/PM pH 7 1.11 0.21 -0.01 1.37 0.87 pH 6 1.19 0.22 0.31 1.43 1.02 pH 5 1.20 0.20 -0.33 1.32 0.76 no chalk 1.09 0.21 -0.32 1.36 0.83 17 N*1 pH 7 0.86 0.15 -1.41 0.98 0.08 pH 6 0.77 0.15 -0.67 0.97 0.05 pH 5 0.87 0.17 -1.20 1.09 0.25 no chalk 0.62 0.16 0.19 1.07 0.22 Predicited x1 (Water stress days) -0.02 0.00 14.66 1.45 n/a environmental x2 (Temperature May-June) -0.18 0.13 14.66 1.45 n/a mean x3 (Temperature July-August) -0.64 0.11 14.66 1.45 n/a supply (e.g. treatment ‘N2 P K Na Mg—pH 6/7’) were identi- suggesting a similar environmental adaptability and yield fied as the most important agronomic practices for sustaining performance. The largest yield susceptibility to lower pH yield under changing climatic conditions (Table 3 and Fig. 5). values was observed in the ‘P K Na Mg’ treatment (possibly A striking finding was that the position and slope of explained by the specific response for legumes, see above), regression lines, except ‘Nil’ (unfertilized) and ‘N*1’ (with and the least reactivity was shown by the ‘N3 P K Na Mg’ −1 48 kg N  ha sodium nitrate), showed the least differentiation treatment, followed by the ‘N2 P K Na Mg’ and FYM/PM’ between treatments in the liming variant ‘pH 6’ (Fig. 5b), treatments (Table 3 and Fig. 5). 1 3 37 Page 14 of 19 J. Macholdt et al. Fig. 5 Graphical visualization Treatment specific yield of the treatment-specific criss- [t/ha] cross regression analyses (see Table 2; detailed explanation (a) pH 7 (b) pH 6 provided in Table A7 Supple- mentary material) depending on liming. a Treatments with Environmental pH 7. b Treatments with pH mean 6. c Treatments with pH 5. d N3PKNaMg Treatments with no chalk. The N2PKNaMg calculation of the environmen- N1PKNaMg tal mean considered the three PKNaMg selected climatic drivers: the 4 FYM/PM ‘accumulated number of water stress days’ and ‘mean air N*1 temperature from May–June Nil and July–August’ (Fig. 2 and Table 1) were used to obtain predicted mean yields for each year (=environmental mean, x-axis) and regressed on the fertilization × liming treat- (c) pH 5 (d) no chalk ment yields (y-axis); getting regression lines with treatment- specific intercepts and slopes. Treatment explanations are provided in Table A3 Sup- plementary material (FYM/ PM = farmyard/poultry manure; N* = sodium nitrate). 0246 810 0246 810 Environmental mean yield [t/ha] Key message for ‘Research question III—part A’: Higher 3.3 Correlation between plant species diversity plant species diversity was correlated not only with more and yield variance stable grassland yields, but also with lower yield levels. In addition, the pCCA of the effect of year on plant com- All correlations were assessed via the non-parametric Kend- munity composition including plot as a covariate (Fig. A12 all’s tau, which provides a robust measure of association. In Supplementary material) explained 16.4% of total variance the PGE, an increased plant species diversity, expressed as (P < 0.001). Years tended to cluster in discrete periods char- Shannon’s diversity index, number of species, and proportion acterized by the dominance of different species supporting of legumes, was found to be associated with lower, but more the idea of environmental perturbation promoting species stable, biomass yields and vice versa (Table 2 and Table A10 coexistence. For example, Crepis capillaris appears to have Supplementary material), which is in line with findings of been favoured by the environmental conditions in 1991 and recent grassland studies (Haughey et al. 2018; Sanderson 1992 and Leontodon autumnalis in 1999 and 2000. The first 2010). For mean yield and temporal yield variance, the cor- axis discriminated between years dominated by grasses and relations (Kendall’s tau) with plant species diversity indices those dominated by forbs. These results imply that the pro- across treatments were significantly negative ( P < 0.05): for ductivity of diverse plots with a more balanced community yield variance and Shannon’s diversity index (τ = −0.30), in terms of the ratio of grasses to forbs will be more resilient species number (τ  =  −0.30), and proportion of legumes to variability in the environment. The pCCA of the effect of (τ = −0.20); for mean yield and Shannon’s diversity index ‘fertilizer × liming treatment’ (sub-plot) on plant community (τ = −0.15), species number (τ = −0.50), and proportion of composition including year as a covariate explained 64.2% legumes (τ = −0.18) (Table A10 Supplementary material). 1 3 Long-term trends in yield variance of temperate managed grassland Page 15 of 19 37 of total variance (P < 0.001) and discriminated between from aluminum toxicity effects on root growth (Kohler et al. unlimed plots with higher N-input (e.g. 9/2d: ‘N2 P K Na 2016), which may have magnified temporal yield variance. Mg—no chalk’; 11/1d: ‘N3 P K Na Mg—no chalk’) that Hence, liming is a very important management factor for are dominated by grasses and the unfertilized plots (3b/a: stabilizing yields and maintaining legumes in soils prone to ‘Nil—pH 6/7’) that had a higher proportion of forbs (Fig. acidification (Fornara et al. 2011; Storkey et al. 2015), which A13 Supplementary material). in this study has also been proven relevant for supporting sta- Although our results are correlative and plant biodiversity ble grassland productivity over time (Table 1 and Fig. 4) and data were only available for a sub-set of years, they are sup- under varying climatic conditions, particularly limiting the ported by evidence showing that greater plant species rich- water balance and higher temperatures from July to August ness and phylogenetic diversity in managed grasslands may (Table 3 and Fig. 5). enhance their resilience to climate change via enhanced asyn- chrony in the performance of co-occurring species and result in more stable biomass production in response to disturbance 4 Conclusion (Isbell et al. 2017; Eisenhauer et al. 2019; Hautier et al. 2020). This assumes that a large species pool is likely, by chance, Overall, our analysis led to the conclusion that liming, fol- to possess one or two stress-tolerant species that are able to lowed by moderate nutrient supply, promoted plant species resist abiotic stress (e.g. drought) (Kahmen et al. 2005). Such diversity, yield stability, and environmental adaptability stress-tolerant species could compensate for less tolerant spe- and enhanced the long-term sustainability of grassland pro- cies and thus help stabilize productivity in grasslands, and this duction (in terms of stable productivity, biodiversity, and could become more important with respect to climate change reduced synthetic fertilizer inputs). The three research ques- (Loreau and de Mazancourt 2013; Trnka et al. 2021; Haughey tions can be answered as the following: et al. 2018). In the PGE, treatments with high proportions of legumes (‘P K Na Mg’ and ‘N1 P K Na Mg’) had mean yields I. Yes, there were temporal trends in interannual yield close to those in plots with higher rates of inorganic N fertilizer variance over the study period and these varied in while also maintaining higher species richness. However, the relation to fertilizer and lime applications. In particu- yield variance in those plots was similar to that in plots with lar, there were consistent unimodal trends observed moderate fertilizer inputs and sufficient liming (treatment ‘N2 in yield variance in plots with low to moderate or no P K Na Mg – pH 7’; Table 2). These results could be explained nitrogen fertilizer additions, with a peak in the 1990s, by the diversity of fast vs slow functional traits (Storkey and after which variability declined. Yield was most vari- Macdonald 2022; Reich and Cornelissen 2014). Grassland able in plots with higher nutrient inputs and lower communities dominated by slow species were found to stabi- soil pH. lize biomass productivity by increasing mean yield relative to II. The accumulated days of water stress from March temporal yield variance (Craven et al. 2018). to October and mean air temperature from July to Key message for ‘Research question III–part B’: Yield August were identified as the most important climatic variance increased, and plant species diversity decreased with drivers, explaining approximately one-third of the greater fertilizer inputs and reduced liming. Treatments with observed interannual yield variance. Yes, yield sen- high proportions of legumes had mean yields close to those sitivity to these climatic drivers depended on agro- in plots with higher rates of inorganic N fertilizer while also nomic management. Liming and moderate N supply maintaining higher species richness. reduced yield sensitivity to abiotic stresses. The relationship between stable productivity (low yield III. Yes, there was a correlation between plant species variance) and plant species diversity strongly varies depend- diversity, mean yield, and yield variance. Higher ing on agronomic management and soil conditions (Bullock plant species diversity was correlated not only with et al. 2001; Hector et al. 1999; Tracy and Sanderson 2004; more stable grassland yields but also with lower yield Crawley et al. 2005; Storkey et al. 2015). For example, spe- levels. Yield variance increased, and plant species cies diversity decreased sharply with reduced soil pH, par- diversity decreased with greater fertilizer inputs and ticularly under enhanced N supply in the form of ammonium reduced liming. sulfate (Crawley et al. 2005). Acidification due to long-term application of ammonium sulfate on the PGE has selected As a limitation of this study, it should be noted that there for a few grass species tolerant of low pH and has signifi- is a lack of replication in the PGE and that plot errors and cantly reduced the proportion of legumes on plot 11/1 (see ‘fertilizer × liming’ interactions (=residual) could not be Table 2, ‘N3 PK Na Mg—no chalk’, pH 3.5). This acidic separated. Furthermore, biodiversity data were only avail- plot was probably subject to physiological stresses imposed able for a sub-set of years and plots. For this reason, our by low pH (Dodd et al. 1994) and drought effects resulting findings should be interpreted carefully and validated by 1 3 37 Page 16 of 19 J. Macholdt et al. Data Availability The yield dataset of the PGE and the mean annual air further detailed analyses including prospective yield data, temperatures from 1965 to 2018 at Rothamsted were provided by the biodiversity surveys, and soil–climate measurements. electronic Rothamsted Archive and are available from the e-RA web- As possible features of future studies, meta-analyses site (see https://doi. or g/10. 23637/ r pg5-yield s1965- 2018- 01 and https:// of various grassland LTEs under different climate and site doi.or g/10. 23637/ r ms-RMAA tem p-02 , respectively). The precipitation chemistry data from 1992 to 2015 were provided by the UK Environ- conditions may provide further valuable information about mental Change Network (ECN) and are publicly available (see https:// temporal trends in yield variance depending on agronomic doi. or g/ 10. 5285/ 18b7c 387- 037d- 4949- 98bc- e8db5 ef426 4c). The management. In addition to retrospective analyses, grassland annual emissions of sulfur dioxide were provided by the UK National LTEs are also a valuable source for application of agroecosys- Atmospheric Emissions Inventory and are publicly available (https:// naei.beis. go v.uk/ dat a/dat a-selec t or?vie w=air pollut ants ). The species of tem models to simulate grassland responses under contrasting grasses, forbs, and legumes comprising at least 5% of the aboveground soil conditions and under future climate scenarios (Qi et al. biomass found in the surveys of the PGE are available from the e-RA 2018), which should be addressed more in upcoming studies. website (see https://doi. or g/10. 23637/ r pg5-speci es_ 1991- 2000- 01 and Overall, the analysis of long-term grassland experiments, http:// www. era. rotha msted. ac. uk/ datas et/ rpg5/ 01- OAPGs pecies). like this study based on the PGE, could help to improve the Code availability The R-scripts for the statistical analyses generated climate resilience and sustainability of grassland systems by during the current study are available in Tables A6–A7 Supplementary identifying climatic drivers and optimizing the agronomic material or from the corresponding author on reasonable request. management accordingly. In particular, the application of the new designed criss-cross regression approach, in which the Declarations environmental mean was modeled using the selected environ- Ethics approval and consent to participate Not applicable. mental abiotic covariates, allows the assessment of the yield sensitivity (or responsiveness) to changes in climatic condi- Consent for publication Not applicable. tions. The application of this criss-cross regression approach in other agro-ecological trials could help to identify climatic Conflict of interest The authors declare no competing interests. drivers of production risk and to derive agronomic manage- Open Access This article is licensed under a Creative Commons Attri- ment strategies for improving the climate resilience of crop- bution 4.0 International License, which permits use, sharing, adapta- ping systems. This will become increasingly important for tion, distribution and reproduction in any medium or format, as long stable agricultural production in the face of climate change as you give appropriate credit to the original author(s) and the source, and the associated growing risk for abiotic stresses. provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are Supplementary Information The online version contains supplementary included in the article's Creative Commons licence, unless indicated material available at https:// doi. org/ 10. 1007/ s13593- 023- 00885-w. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not Acknowledgements We thank the Rothamsted farm staff, who man- permitted by statutory regulation or exceeds the permitted use, you will aged the PGE and the large teams of people who were involved in the need to obtain permission directly from the copyright holder. To view a fieldwork, vegetation sampling, and sorting between 1965 and 2018. copy of this licence, visit http://cr eativ ecommons. or g/licen ses/ b y/4.0/ . We thank the UK Biotechnology and Biological Sciences Research Council (BBS/E/C/000J0300) and the Lawes Agricultural Trust for supporting the PGE. We thank the Lawes Agricultural Trust and Rothamsted Research for data from the e-RA database. We thank the References UK Environmental Change Network (ECN) for the Precipitation Chem- istry Data. Further thanks are extended to the anonymous reviewers and Addy JWG, Ellis RH, Macdonald AJ, Semenov MA, Mead A (2021) editors for their feedback on earlier versions of this paper. Changes in agricultural climate in South-Eastern England from 1892 to 2016 and differences in cereal and permanent grassland Authors' contributions Janna Macholdt initiated and conceived the yield. Agr Forest Meteorol 308–309:108560. https:// doi. org/ 10. study. All the authors contributed to the further research process. Sarah 1016/j. agrfo rmet. 2021. 108560 Perryman and Tony Scott compiled the datasets. Janna Macholdt, Stef- Addy JWG, Ellis RH, MacLaren C, Macdonald AJ, Semenov MA, fen Hadasch, Hans-Peter Piepho, and Merete Styczen designed and Mead A (2022) A heteroskedastic model of Park Grass spring carried out the statistical analysis. The first draft of the manuscript was hay yields in response to weather suggests continuing yield written by Janna Macholdt, and all the authors commented on previous decline with climate change in future decades. J R Soc Interface versions of the manuscript. All the authors have read and approved the 19(193):20220361. https:// doi. org/ 10. 1098/ rsif. 2022. 0361 final manuscript. Anon (2006) Guide to the classical and other long-term experiments, datasets and sample archive. Lawes Agricultural Trust Ltd, Funding Open Access funding enabled and organized by Projekt DEAL. Harpenden UK. https:// doi. org/ 10. 23637/ rotha msted- long- term- The PGE is supported by the UK Biotechnology and Biological Sciences exper iments- guide- 2006 Research Council (BBS/E/C/000J0300) and the Lawes Agricultural Arata L, Fabrizi E, Sckokai P (2020) A worldwide analysis of trend in Trust. The authors Janna Macholdt and Hans-Peter Piepho acknowledge crop yields and yield variability: evidence from FAO data. Econ support by DFG (Deutsche Forschungsgemeinschaft) grants MA Model 90:190–208. https://doi. or g/10. 1016/j. econm od. 2020. 05. 006 7094/1-1 and PI 377/20-2, respectively. The author Jonathan Storkey Avery BW, Catt JA (1995) The soil at Rothamsted. 44 e-RA. Lawes acknowledges the support by the Biological Sciences Research Council Agricultural Trust Co. Ltd, Harpenden UK. https:// doi. org/ 10. (BBSRC) funded Soils to Nutrition project (BBS/E/C/000I0320). 23637/ ERADOC- 1- 143 1 3 Long-term trends in yield variance of temperate managed grassland Page 17 of 19 37 Baca Cabrera JC, Hirl RT, Schaufele R, Macdonald A, Schnyder H O, Freschet GT, Giling DP, Hattenschwiler S, Hillebrand H, (2021) Stomatal conductance limited the CO response of grass- Hines J, Isbell F, Koller-France E, Konig-Ries B, de Kroon H, land in the last century. BMC Biol 19(1):50. https:// doi. org/ 10. Meyer ST, Milcu A, Muller J, Nock CA, Petermann JS, Roscher 1186/ s12915- 021- 00988-4 C, Scherber C, Scherer-Lorenzen M, Schmid B, Schnitzer SA, Bengtsson J, Bullock JM, Egoh B, Everson C, Everson T, O'Connor Schuldt A, Tscharntke T, Turke M, van Dam NM, van der Plas T, O'Farrell PJ, Smith HG, Lindborg R (2019) Grasslands-more F, Vogel A, Wagg C, Wardle DA, Weigelt A, Weisser WW, important for ecosystem services than you might think. Ecosphere Wirth C, Jochum M (2019) A multitrophic perspective on biodi- 10(2). https:// doi. org/ 10. 1002/ ecs2. 2582 versity-ecosystem functioning research. Adv Ecol Res 61:1–54. Berti A, Marta AD, Mazzoncini M, Tei F (2016) An overview on long-term https:// doi. org/ 10. 1016/ bs. aecr. 2019. 06. 001 agro-ecosystem experiments: present situation and future potential. Fornara DA, Steinbeiss S, McNamara NP, Gleixner G, Oakley S, Eur J Agron 77:236–241. https:// doi. org/ 10. 1016/j. eja. 2016. 01. 004 Poulton PR, Macdonald AJ, Bardgett RD (2011) Increases in Bommarco R, Kleijn D, Potts SG (2013) Ecological intensification: soil organic carbon sequestration can reduce the global warm- harnessing ecosystem services for food security. Trends Ecol Evol ing potential of long-term liming to permanent grassland. Glob 28(4):230–238. https:// doi. org/ 10. 1016/j. tree. 2012. 10. 012 Chang Biol 17(5):1925–1934. https:// doi. org/ 10. 1111/j. 1365- Brookshire EN, Weaver T (2015) Long-term decline in grassland 2486. 2010. 02328.x productivity driven by increasing dryness. Nat Commun 6:7148. Graux A-I, Bellocchi G, Lardy R, Soussana J-F (2013) Ensemble https:// doi. org/ 10. 1038/ ncomm s8148 modelling of climate change risks and opportunities for man- Bullock JM, Pywell RF, Burke MJW, Walker KJ (2001) Restoration of aged grasslands in France. Agric for Meteorol 170:114–131. biodiversity enhances agricultural production. Ecol Lett 4(3):185–https:// doi. org/ 10. 1016/j. agrfo rmet. 2012. 06. 010 189. https:// doi. org/ 10. 1046/j. 1461- 0248. 2001. 00215.x Hadasch S, Laidig F, Macholdt J, Bönecke E, Piepho HP (2020) Carswell AM, Gongadze K, Misselbrook TH, Wu L (2019) Impact of Trends in mean performance and stability of winter wheat and transition from permanent pasture to new swards on the nitrogen winter rye yields in a long-term series of variety trials. Field use efficiency, nitrogen and carbon budgets of beef and sheep Crops Res 252. https:// doi. org/ 10. 1016/j. fcr. 2020. 107792 production. Agric Ecosyst Environ 283:106572. https:// doi. org/ Hall DO, Scurlock JMO (1991) Climate change and productivity 10. 1016/j. agee. 2019. 106572 of natural grassland. Ann Bot 67:49–55. https:// www. jstor. org/ Chang J, Ciais P, Viovy N, Soussana JF, Klumpp K, Sultan B (2017) stable/ 42758 390 Future productivity and phenology changes in European grass- Haughey E, Suter M, Hofer D, Hoekstra NJ, McElwain JC, Luscher A, lands for different warming levels: implications for grassland Finn JA (2018) Higher species richness enhances yield stability management and carbon balance. Carbon Balance Manag in intensively managed grasslands with experimental disturbance. 12(1):11. https:// doi. org/ 10. 1186/ s13021- 017- 0079-8 Sci Rep 8(1):15047. https://doi. or g/10. 1038/ s41598- 018- 33262-9 Craven D, Eisenhauer N, Pearse WD, Hautier Y, Isbell F, Roscher Hautier Y, Seabloom EW, Borer ET, Adler PB, Harpole WS, Hillebrand C, Bahn M, Beierkuhnlein C, Bonisch G, Buchmann N, Byun H, Lind EM, MacDougall AS, Stevens CJ, Bakker JD, Buckley C, Catford JA, Cerabolini BEL, Cornelissen JHC, Craine JM, YM, Chu C, Collins SL, Daleo P, Damschen EI, Davies KF, Fay De Luca E, Ebeling A, Griffin JN, Hector A, Hines J, Jentsch A, PA, Firn J, Gruner DS, Jin VL, Klein JA, Knops JM, La Pierre Kattge J, Kreyling J, Lanta V, Lemoine N, Meyer ST, Minden V, KJ, Li W, McCulley RL, Melbourne BA, Moore JL, O’Halloran Onipchenko V, Polley HW, Reich PB, van Ruijven J, Schamp B, LR, Prober SM, Risch AC, Sankaran M, Schuetz M, Hector A Smith MD, Soudzilovskaia NA, Tilman D, Weigelt A, Wilsey (2014) Eutrophication weakens stabilizing effects of diversity in B, Manning P (2018) Multiple facets of biodiversity drive the natural grasslands. Nature 508(7497):521–525. https://doi. or g/10. diversity-stability relationship. Nat Ecol Evol 2(10):1579–1587. 1038/ natur e13014 https:// doi. org/ 10. 1038/ s41559- 018- 0647-7 Hautier Y, Zhang P, Loreau M, Wilcox KR, Seabloom EW, Borer Crawley MJ, Johnston AE, Silvertown J, Dodd M, de Mazancourt ET, Byrnes JEK, Koerner SE, Komatsu KJ, Lefcheck JS, Hector C, Heard MS, Henman DF, Edwards GR (2005) Determinants A, Adler PB, Alberti J, Arnillas CA, Bakker JD, Brudvig LA, of species richness in the Park Grass Experiment. Am Nat Bugalho MN, Cadotte M, Caldeira MC, Carroll O, Crawley M, 165(2):179–192. https:// doi. org/ 10. 1086/ 427270 Collins SL, Daleo P, Dee LE, Eisenhauer N, Eskelinen A, Fay Deguines N, Jono C, Baude M, Henry M, Julliard R, Fontaine C PA, Gilbert B, Hansar A, Isbell F, Knops JMH, MacDougall AS, (2014) Large-scale trade-off between agricultural intensification McCulley RL, Moore JL, Morgan JW, Mori AS, Peri PL, Pos ET, and crop pollination services. Front Ecol Environ 12(4):212– Power SA, Price JN, Reich PB, Risch AC, Roscher C, Sankaran 217. https:// doi. org/ 10. 1890/ 130054 M, Schutz M, Smith M, Stevens C, Tognetti PM, Virtanen R, Digby PGN (2009) Modified joint regression analysis for incom - Wardle GM, Wilfahrt PA, Wang S (2020) General destabilizing plete variety × environment data. The J Agric Sci 93(1):81–86. effects of eutrophication on grassland productivity at multiple https:// doi. org/ 10. 1017/ S0021 85960 00861 59 spatial scales. Nat Commun 11(1):5375. https:// doi. org/ 10. 1038/ Dodd ME, Silvertown J, McConway K, Potts J, Crawley M (1994) s41467- 020- 19252-4 Application of the British national vegetation classification to Hector A, Schmid B, Beierkuhnlein C, Caldeira MC, Diemer M, Dimi- the communities of the park grass experiment through time. trakopoulos PG, Finn JA, Freitas H, Giller PS, Good J, Harris R, Folia Geobot Phytotx 29(3):321–334. https:// doi. org/ 10. 1007/ Hogberg P, Huss-Danell K, Joshi J, Jumpponen A, Korner C, Lead- bf028 82911 ley PW, Loreau M, Minns A, Mulder CP, O’Donovan G, Otway SJ, Dodd M, Silvertown J, McConway K, Potts J, Crawley M (1997) Sta- Pereira JS, Prinz A, Read DJ, Scherer-Lorenzen M, Schulze ED, bility in the plant communities of the Park Grass Experiment: Siamantziouras ASD, Spehn EM, Terry AC, Troumbis AY, Wood- the relationships between species richness, soil pH and biomass ward FI, Yachi S, Lawton JH (1999) Plant diversity and productiv- variability. Philos Trans R Soc Lond B Biol Sci 346(1316):185– ity experiments in European grasslands. Science 286(5442):1123– 193. https:// doi. org/ 10. 1098/ rstb. 1994. 0140 1127. https:// doi. org/ 10. 1126/ scien ce. 286. 5442. 1123 Eberhart SA, Russell WA (1966) Stability Parameters for Comparing Höglind M, Thorsen SM, Semenov MA (2013) Assessing uncertain- Varieties. Crop Sci 6(1):36–40. https:// doi. org/ 10. 2135/ crops ties in impact of climate change on grass production in North- ci1966. 00111 83X00 06000 10011x ern Europe using ensembles of global climate models. Agric for Eisenhauer N, Schielzeth H, Barnes AD, Barry K, Bonn A, Brose Meteorol 170:103–113. https://doi. or g/10. 1016/j. ag rfor met.2012. 02. 010 U, Bruelheide H, Buchmann N, Buscot F, Ebeling A, Ferlian 1 3 37 Page 18 of 19 J. Macholdt et al. Holland JE, Bennett AE, Newton AC, White PJ, McKenzie BM, George Perryman S, Scott T, Hall C (2019) Rothamsted 30-year meteorological TS, Pakeman RJ, Bailey JS, Fornara DA, Hayes RC (2018) Liming means 1981–2010. Electronic Rothamsted Archive, Rothamsted impacts on soils, crops and biodiversity in the UK: a review. Sci Research. https:// doi. org/ 10. 23637/ OARES 30YrM eans8 180 Total Environ 610–611:316–332. https:// doi. org/ 10. 1016/j. scito Perryman S, Crawley M, Ostler R, Storkey J (2021) Dataset: Park Grass tenv. 2017. 08. 020 species, fertilizer and lime treatments 1991–2000. Electronic Isbell F, Adler PR, Eisenhauer N, Fornara D, Kimmel K, Kremen C, Rothamsted Archive, Rothamsted Research. https:// doi. org/ 10. Letourneau DK, Liebman M, Polley HW, Quijas S, Scherer-Lor-23637/ rpg5- speci es_ 1991- 2000- 01 enzen M, Bardgett R (2017) Benefits of increasing plant diversity Piepho HP (1998) Methods for comparing the yield stability of crop- in sustainable agroecosystems. J Ecol 105(4):871–879. https://doi. ping systems—a review. J Agron Crop Sci 180(4):193–213. org/ 10. 1111/ 1365- 2745. 12789https:// doi. org/ 10. 1111/j. 1439- 037X. 1998. tb005 26.x Kahmen A, Perner J, Buchmann N (2005) Diversity-dependent produc- Piseddu F, Bellocchi G, Picon-Cochard C (2021) Mowing and warm- tivity in semi-natural grasslands following climate perturbations. ing effects on grassland species richness and harvested biomass: Funct Ecol 19(4):594–601. https:// doi.or g/ 10. 1111/j.13 65-24 35. meta-analyses. Agron Sustain Dev 41(6). https://d oi.o rg/1 0.10 07/ 2005. 01001.xs13593- 021- 00722-y Kipling RP, Virkajarvi P, Breitsameter L, Curnel Y, De Swaef T, Preissel S, Reckling M, Schläfke N, Zander P (2015) Magnitude and Gustavsson AM, Hennart S, Hoglind M, Jarvenranta K, Minet farm-economic value of grain legume pre-crop benefits in Europe: J, Nendel C, Persson T, Picon-Cochard C, Rolinski S, Sandars a review. Field Crops Res 175:64–79. https:// doi. org/ 10. 1016/j. DL, Scollan ND, Sebek L, Seddaiu G, Topp CFE, Twardy S, Van fcr. 2015. 01. 012 Middelkoop J, Wu L, Bellocchi G (2016) Key challenges and pri- Prieto I, Violle C, Barre P, Durand JL, Ghesquiere M, Litrico I (2015) orities for modelling European grasslands under climate change. Complementary effects of species and genetic diversity on pro - Sci Total Environ 566–567:851–864. https:// doi. org/ 10. 1016/j. ductivity and stability of sown grasslands. Nat Plants 1(4):15033. scito tenv. 2016. 05. 144https:// doi. org/ 10. 1038/ NPLAN TS. 2015. 33 Kohler IH, Macdonald AJ, Schnyder H (2016) Last-century increases Qi A, Holland RA, Taylor G, Richter GM (2018) Grassland futures in in intrinsic water-use efficiency of grassland communities fave Great Britain—productivity assessment and scenarios for land use occurred over a wide range of vegetation composition, nutrient change opportunities. Sci Total Environ 634:1108–1118. https:// inputs, and soil pH. Plant Physiol 170(2):881–890. https:// doi. doi. org/ 10. 1016/j. scito tenv. 2018. 03. 395 org/ 10. 1104/ pp. 15. 01472 Raman A, Ladha JK, Kumar V, Sharma S, Piepho HP (2011) Stability Lawes JB, Gilbert JH (1859) Report of experiments with different analysis of farmer participatory trials for conservation agriculture manures on permanent meadow land: part i: produce of hay per using mixed models. Field Crops Res 121(3):450–459. https://doi. acre. J Roy Agr Soc Engl 19:552–573org/ 10. 1016/j. fcr. 2011. 02. 001 Loreau M, de Mazancourt C (2013) Biodiversity and ecosystem stabil- Ray DK, Gerber JS, MacDonald GK, West PC (2015) Climate varia- ity: a synthesis of underlying mechanisms. Ecol Lett 16:106–115. tion explains a third of global crop yield variability. Nat Commun https:// doi. org/ 10. 1111/ ele. 12073 6:5989. https:// doi. org/ 10. 1038/ ncomm s6989 Macdonald A, Poulton P, Clark I, Scott T, Glendining M, Perryman Ray DK, West PC, Clark M, Gerber JS, Prishchepov AV, Chatterjee S, Storkey J, Bell J, Shield I, McMillan V, Hawkins J (2018) S (2019) Climate change has likely already affected global food Rothamsted long-term experiments—guide to the classical production. Plos One 14(5):e0217148. https:// doi. org/ 10. 1371/ and other long-term experiments, datasets and sample archive. journ al. pone. 02171 48 Rothamsted Research (UK). https:// www. rotha msted. ac. uk/ sites/ Reckling M, Ahrends H, Chen T-W, Eugster W, Hadasch S, Knapp defau lt/ files/ RRes_ LTE% 20Gui debook_ 2018_% 20web% 20AW. S, Laidig F, Linstädter A, Macholdt J, Piepho H-P, Schiffers K, pdf. Accessed 24 Feb 2023 Döring TF (2021) Methods of yield stability analysis in long-term Macholdt J, Hadasch S, Piepho HP, Reckling M, Taghizadeh-Toosi A, field experiments. A review. Agron Sustain Dev 41(2). https://doi. Christensen BT (2021) Yield variability trends of winter wheat org/ 10. 1007/ s13593- 021- 00681-4 and spring barley grown during 1932–2019 in the Askov Long- Reich PB, Cornelissen H (2014) The world-wide ‘fast-slow’ plant term Experiment. Field Crops Res 264:108083. https:// doi. org/ economics spectrum: a traits manifesto. J Ecol 102(2):275–301. 10. 1016/j. fcr. 2021. 108083https:// doi. org/ 10. 1111/ 1365- 2745. 12211 Midolo G, Alkemade R, Schipper AM, Benítez-López A, Perring Robertson GP, Vitousek PM (2009) Nitrogen in agriculture: balancing MP, De Vries W, Xu X (2018) Impacts of nitrogen addition on the cost of an essential resource. Annu Rev Environ Resour 34:97– plant species richness and abundance: a global meta-analysis. 125. https:// doi. org/ 10. 1146/ annur ev. envir on. 032108. 105046 Glob Ecol Biogeogr 28(3):398–413. https:// doi. or g/ 10. 1111/ Römer T (1917) Are high-yielding varieties more yield stable? Com- geb. 12856 mun DLG 31:87–89 Nabugoomu F, Kempton RA, Talbot M (1999) Analysis of series of Sanderson MA (2010) Stability of production and plant species diver- trials where varieties differ in sensitivity to locations. J Agric Biol sity in managed grasslands: a retrospective study. Basic Appl Ecol Environ Stat 4(3). https:// doi. org/ 10. 2307/ 14003 88 11(3):216–224. https:// doi. org/ 10. 1016/j. baae. 2009. 08. 002 Olesen JE, Bindi M (2002) Consequences of climate change for Euro- Schmidhuber J, Tubiello FN (2007) Global food security under climate pean agricultural productivity, land use and policy. Eur J Agron change. Proc Natl Acad Sci U S A 104(50):19703–19708. https:// 16(4):239–262. https:// doi. org/ 10. 1016/ S1161- 0301(02) 00004-7doi. org/ 10. 1073/ pnas. 07019 76104 Onofri A, Seddaiu G, Piepho H-P (2016) Long-Term Experiments Silvertown J, Poulton P, Johnston E, Edwards G, Heard M, Biss PM with cropping systems: case studies on data analysis. Eur J Agron (2006) The Park Grass Experiment 1856–2006: its contribution 77:223–235. https:// doi. org/ 10. 1016/j. eja. 2016. 02. 005 to ecology. J Ecol 94(4):801–814. https://d oi.o rg/1 0.1 111/j.1 365- Payne RW (2015) The design and analysis of long-term rotation experi-2745. 2006. 01145.x ments. Agron J 107(2):772–785. https:// doi. org/ 10. 2134/ agron Storkey J, Macdonald AJ (2022) The role of long-term experiments in j2012. 0411 validating trait-based approaches to achieving multifunctionality Perryman S, Ostler R (2021) Dataset: Park Grass hay yields, ferti- in grasslands. Front Agric Sci Eng 9(2):187–196. https:// doi. org/ lizer and lime treatments 1965–2018. Electronic Rothamsted 10. 15302/j- fase- 20214 38 Storkey J, Macdonald AJ, Poulton PR, Scott T, Kohler IH, Schnyder H, Archive, Rothamsted Research https:// doi. org/ 10. 23637/ r pg5- Goulding KW, Crawley MJ (2015) Grassland biodiversity bounces yield s1965- 2018- 01 1 3 Long-term trends in yield variance of temperate managed grassland Page 19 of 19 37 back from long-term nitrogen addition. Nature 528(7582):401– Verbyla A, Cullis B, Kenward M, Welham S (1999) The analysis of 404. https:// doi. org/ 10. 1038/ natur e16444 designed experiments and longitudinal data by using smoothing Tilman D, Reich PB, Knops J, Wedin D, Mielke T, Lehman C (2001) splines. J R Stat Soc Series C Stat 48(3):269–311. https:// www. Diversity and productivity in a long-term grassland experiment. Sci-jstor. org/ stable/ 26808 26 ence 294(5543):843–845. https:// doi. org/ 10. 1126/ scien ce. 10603 91 Wilcox KR, Shi Z, Gherardi LA, Lemoine NP, Koerner SE, Hoover Tilman D, Reich PB, Knops JMH (2006) Biodiversity and ecosys- DL, Bork E, Byrne KM, Cahill J Jr, Collins SL, Evans S, Gilgen tem stability in a decade-long grassland experiment. Nature AK, Holub P, Jiang L, Knapp AK, LeCain D, Liang J, Garcia- 441(7093):629–632. https:// doi. org/ 10. 1038/ natur e04742 Palacios P, Penuelas J, Pockman WT, Smith MD, Sun S, White Tracy BF, Sanderson MA (2004) Productivity and stability relation- SR, Yahdjian L, Zhu K, Luo Y (2017) Asymmetric responses ships in mowed pasture communities of varying species compo- of primary productivity to precipitation extremes: a synthesis of sition. Crop Sci 44(6):2180–2186. https:// doi. org/ 10. 2135/ crops grassland precipitation manipulation experiments. Glob Chang ci2004. 2180 Biol 23(10):4376–4385. https:// doi. org/ 10. 1111/ gcb. 13706 Trnka M, Olesen JE, Kersebaum KC, SkjelvÅG AO, Eitzinger J, Zhang Y, Loreau M, Lu X, He N, Zhang G, Han X (2016) Nitrogen Seguin B, Peltonen-Sainio P, RÖTter R, Iglesias ANA, Orlandini enrichment weakens ecosystem stability through decreased spe- S, DubrovskÝ M, Hlavinka P, Balek J, Eckersten H, Cloppet E, cies asynchrony and population stability in a temperate grassland. Calanca P, Gobin A, VuČEtiĆ V, Nejedlik P, Kumar S, Lalic B, Glob Chang Biol 22(4):1445–1455. https:// doi. org/ 10. 1111/ gcb. Mestre A, Rossi F, Kozyra J, Alexandrov V, SemerÁDovÁ D, 13140 ŽAlud Z (2011) Agroclimatic conditions in Europe under climate Zhang Y, Feng J, Loreau M, He N, Han X, Jiang L (2019) Nitrogen change. Glob Chang Biol 17(7):2298–2318. https:// doi. org/ 10. addition does not reduce the role of spatial asynchrony in stabilis- 1111/j. 1365- 2486. 2011. 02396.x ing grassland communities. Ecol Lett 22(4):563–571. https://doi. Trnka M, Balek J, Semenov MA, Semeradova D, Belinova M, Hlavinka org/ 10. 1111/ ele. 13212 P, Olesen JE, Eitzinger J, Schaumberger A, Zahradnicek P, Kopecky D, Zalud Z (2021) Future agroclimatic conditions and implications Publisher's note Springer Nature remains neutral with regard to for European grasslands. Biol Plant 64:865–880. https://doi. or g/10. jurisdictional claims in published maps and institutional affiliations. 32615/ bp. 2021. 005 1 3 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Agronomy for Sustainable Development Springer Journals

Long-term trends in yield variance of temperate managed grassland

Download PDF
19 pages

Loading next page...
 
/lp/springer-journals/long-term-trends-in-yield-variance-of-temperate-managed-grassland-rOQUNXD0db

References (86)

Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2023
ISSN
1774-0746
eISSN
1773-0155
DOI
10.1007/s13593-023-00885-w
Publisher site
See Article on Publisher Site

Abstract

The management of climate-resilient grassland systems is important for stable livestock fodder production. In the face of climate change, maintaining productivity while minimizing yield variance of grassland systems is increasingly challenging. To achieve climate-resilient and stable productivity of grasslands, a better understanding of the climatic drivers of long- term trends in yield variance and its dependence on agronomic inputs is required. Based on the Park Grass Experiment at Rothamsted (UK), we report for the first time the long-term trends in yield variance of grassland (1965–2018) in plots given different fertilizer and lime applications, with contrasting productivity and plant species diversity. We implemented a statisti - cal model that allowed yield variance to be determined independently of yield level. Environmental abiotic covariates were included in a novel criss-cross regression approach to determine climatic drivers of yield variance and its dependence on agronomic management. Our findings highlight that sufficient liming and moderate fertilization can reduce yield variance while maintaining productivity and limiting loss of plant species diversity. Plots receiving the highest rate of nitrogen ferti- lizer or farmyard manure had the highest yield but were also more responsive to environmental variability and had less plant species diversity. We identified the days of water stress from March to October and temperature from July to August as the two main climatic drivers, explaining approximately one-third of the observed yield variance. These drivers helped explain consistent unimodal trends in yield variance—with a peak in approximately 1995, after which variance declined. Here, for the first time, we provide a novel statistical framework and a unique long-term dataset for understanding the trends in yield variance of managed grassland. The application of the criss-cross regression approach in other long-term agro-ecological trials could help identify climatic drivers of production risk and to derive agronomic strategies for improving the climate resilience of cropping systems. Keywords Agronomic management · Biomass production · Climate resilience · Fertilizer input · Food security · Liming · Plant species diversity · Soil pH · Temperature · Water stress * Janna Macholdt 1 Introduction janna.macholdt@landw.uni-halle.de The management of climate-resilient grassland systems is Professorship of Agronomy, Institute of Agriculture and Nutritional Sciences, Martin-Luther-University Halle- important for stable livestock fodder production over time Wittenberg, Betty-Heimann-Strasse 5, 06120 Halle (Saale), (Schmidhuber and Tubiello 2007; Arata et al. 2020; Trnka Germany et al. 2021; Bengtsson et al. 2019; Bommarco et al. 2013; Biostatistics Unit, Institute of Crop Science, University Reckling et al. 2021). However, in the face of climate change of Hohenheim, Fruwirthstrasse 23, 70599 Stuttgart, Germany and the associated increases in abiotic stresses, maintain- Protecting Crops and Environment, Rothamsted Research, ing productivity while minimizing temporal yield variance Harpenden AL5 2JQ, Hertfordshire, UK (or rather improving yield stability) of grassland systems Computational and Analytical Sciences Department, will become increasingly challenging (Olesen and Bindi Rothamsted Research, Harpenden  AL5 2JQ, Hertfordshire, 2002; Ray et al. 2019). Observed climatic changes, such as UK increasing temperatures or weather anomalies, have nega- Section of Environmental Chemistry and Physics, tively affected global grassland productivity (Kipling et al. Department of Plant and Environmental Sciences, University 2016; Höglind et al. 2013; Addy et al. 2022; Brookshire of Copenhagen, Thorvaldsensvej 40, 1871 Copenhagen, Denmark Vol.:(0123456789) 1 3 37 Page 2 of 19 J. Macholdt et al. and Weaver 2015; Wilcox et al. 2017; Hall and Scurlock et al. 1997) because estimates of yield variance based on 1991). For northern Europe as the experimental region of a shorter period might be imprecise and long-term trends this study, Trnka et al. (2021) reported that the annual aver- cannot be detected (Hadasch et al. 2020; Macholdt et al. age temperature as well as the frequency of combined heat 2021; Reckling et al. 2021). Data over decadal time scales and drought stress for grassland has increased in recent also allow the adaptation of grasslands to climatic trends decades. In addition, they showed that by 2050, the area (such as increasing temperatures) to be quantified in addi - of grassland exposed to combined heat and drought may tion to the buffering of short-term climatic variability. double compared with that in 2021 and that northern regions Long-term experiments (LTEs) and their associated data- will experience higher temperatures and drought more fre- sets provide a unique opportunity to examine the effects quently (Trnka et al. 2021). As well as long-term temporal of climate change (Berti et al. 2016). trends, these changes in climatic conditions are expected to Appropriate statistical methods are necessary to handle increase the interannual yield variance of grassland systems the often complex design of LTEs and any experimental and decrease plant species diversity, which could put future modifications that have occurred over time (Reckling et al. food security at risk (Graux et al. 2013; Piseddu et al. 2021). 2021; Payne 2015; Macdonald et al. 2018), such as those The intensification of grassland production through the used for genotype × environment × management interac- use of inorganic fertilizers and simplified swards of fast- tion analyses in other disciplines, including plant breed- growing grass cultivars has delivered increased primary bio- ing (Hadasch et al. 2020) and agronomy (Macholdt et al. mass and live weight gain of livestock (Carswell et al. 2019). 2021). Although there are some recent studies on the yield However, the use of inorganic fertilizer is a major contribu- variance of field crops, detailed knowledge about the long- tor to greenhouse gas emissions (both in manufacturing and term effects of climatic changes and agronomic manage- application) and results in swards with low plant species ment on temporal trends in yield variance for grasslands is diversity that may also compromise resilience (Robertson limited (Dodd et al. 1997). Previous studies reporting yield and Vitousek 2009; Storkey et al. 2015). Recent findings variance of grasslands often cover only short periods with of short-term grassland studies have confirmed that addi- less than 10 years (Tilman et al. 2006; Prieto et al. 2015; tional nutrient inputs often increase not only yield but also Hautier et  al. 2014, 2020; Haughey et  al. 2018; Zhang interannual yield variance and reduce plant species diversity et al. 2016. 2019; Sanderson 2010; Dodd et al. 1997). To (Hautier et al. 2020, 2014; Zhang et al. 2016, 2019; Crawley date, there have been only few analyses of permanent, et al. 2005). Moreover, studies reported that relative climatic managed grassland systems that include temporal trends in adaptability decreased with increasing land use intensity yield variance covering a long-term period (Craven et al. (Deguines et al. 2014). In particular, mineral N (nitrogen) 2018; Isbell et al. 2017; Dodd et al. 1997). These studies fertilizer input reduces the diversity of terrestrial vegeta- often neglect a further requirement, which is that yield tion by favouring fast-growing grass species adapted to high variance should be determined independently of yield nutrient availability (Midolo et al. 2018). There is evidence level, otherwise yield variance can be incorrectly inter- that greater plant species diversity in managed grasslands preted if the time span is too short and there is a systematic may enhance their resilience to climate change and result dependency of variation on the mean (Preissel et al. 2015). in more stable yields in response to disturbance (Tracy and In this study, we used a long-term dataset (1965–2018) Sanderson 2004; Tilman et al. 2006; Haughey et al. 2018; from The Park Grass Experiment at Rothamsted (UK) with Sanderson 2010; Baca Cabrera et al. 2021). However, the consistent long-term management and determined yield relationship between stable productivity (low interannual variance independently of yield level, for most accurate yield variance) and plant species diversity can also strongly statistical estimates. The novelty of our analysis is that we vary depending on agronomic management and soil condi- included environmental abiotic covariates, such as atmos- tions (Bullock et al. 2001; Hector et al. 1999; Tracy and pheric chemistry (wet and dry N deposition, SO, CO ) 2 2 Sanderson 2004; Crawley et al. 2005; Storkey et al. 2015). A and other climatic parameters (air temperature, precipi- very important agronomic management factor for stabilizing tation, soil moisture deficit, etc.), in a novel criss-cross yields over time and for maintaining plant species diversity regression approach (extended Finlay–Wilkinson regres- is liming, particularly in soils prone to acidification (Fornara sion) to determine the main climatic drivers of long-term et al. 2011; Storkey et al. 2015). yield variance and to evaluate the effects of liming and To achieve resilient, sustainable, and stable productiv- fertilizer applications on the relative impact of environ- ity of grasslands, a better understanding of the climatic mental variability on yield sensitivity (or responsiveness) drivers of long-term trends in temporal yield variance and across a range of plots with contrasting productivity and its dependence on agronomic inputs and biological diver- plant species richness. sity is required. Ideally, such assessments should be done We specifically addressed the following three research on long-term datasets (>20 years) (Piepho 1998; Dodd questions: 1 3 Long-term trends in yield variance of temperate managed grassland Page 3 of 19 37 I. Are there temporal trends in interannual yield vari- provided in Fig. A1-b Supplementary material and available ance over the study period, and do these vary in rela- online (https://doi. or g/10. 23637/ ERADOC- 1- 143 ; p. 43). In tion to fertilizer and lime applications? 1856, starting values for PGE soil parameters in the topsoil II. What are the most important environmental abiotic (0–23 cm) were estimated as pH 5.7, 11.6% sand, 66.3% silt, −1 3 drivers explaining yield variance, and does related 22.1% clay, soil weight 2430 t  ha , bulk density 1.1 g  cm , −1 yield sensitivity to these climatic drivers depend on and total soil nitrogen content of 5830 kg N  ha (Lawes agronomic management? and Gilbert 1859). III. Is there a correlation between plant species diversity, The mean annual air temperature and rainfall at the mean yield, and yield variance? Rothamsted site (1981–2010) were 9.8  °C and 733  mm, respectively (Perryman et  al. 2019). In recent decades, the annual mean air temperature (1989–2018) was 1.1 °C 2 Material and methods warmer compared to the previous period (1878–1988: 9.04 °C) (Macdonald et al. 2018). The increasing trend in 2.1 T he Park Grass Experiment (PGE) temperature anomaly at Rothamsted (1880–2018) is pro- vided in Fig.  A2 Supplementary material and available This study was based on the Park Grass Experiment (PGE) online (https:// doi. org/ 10. 23637/ rms- RMAAt empan omaly- at Rothamsted (Fig. 1), initiated by Lawes and Gilbert in 1). The effective plant available water capacity of 135 mm 1856 to examine the effects of different mineral fertilizers in the effective root zone of the Park Grass site implies that and organic manures on the productivity of permanent pas- plant growth is often limited by lack of water in summer ture cut for hay (Silvertown et al. 2006). A detailed descrip- (Avery and Catt 1995). tion of the experiment is available in the Rothamsted Guide The current analysis is focused on the period from 1965 to Classical Experiments (Macdonald et al. 2018) and the to 2018 using plots with constant treatments and a consist- e-RA website (http://www .er a.r othams ted.ac. uk/ e xperiment/ ent harvesting methodology to ensure comparability and rpg5). The Park Grass soil is a moderately well-drained silty accurate estimates of temporal trends in yield variance. The clay loam overlying clay-with-flints (Avery and Catt 1995), design of the PGE is provided in Fig. 1 (more details are a Chromic or Vertic Luvisol according to the FAO clas- provided in Fig. A1-a Supplementary material), consisting sification. The experimental site shows a relatively uniform of 24 main plots with contrasting fertilizer treatments, each soil, based on comprehensive soil analyses made by Avery divided into four sub-plots ‘a, b, c, and d’ for liming treat- and Catt (1995)—the ‘Soils at Rothamsted Colour Map’ is ments. Sub-plots a, b, and c receive lime, if needed, every 3 Fig. 1 Park Grass Experi- ment aerial view (left) and plot layout (right). Loca- tion: Harpenden, UK, Herts, AL5 2JQ (51°48′12.33″N; 0°22′21.66″W; 130 m a.s.l.). Detailed information about plot layout and treatments are shown in Tables A1 and A3 Sup- plementary material. Source: electronic Rothamsted Archive (http:// www. era. rotha msted. ac. uk/ Park# images/). 1 3 37 Page 4 of 19 J. Macholdt et al. years to maintain soil pH 7, 6, and 5, respectively; sub-plot d (Crawley et al. 2005; Storkey et al. 2015; Ray et al. 2015). receives no chalk. The yield data selected for this study were Although plant communities respond to these drivers, the taken from the a, b, c, and d sub-plots of seven contrasting relative differences between the plots in species richness main plots (red marked in Fig. A1 Supplementary material) and diversity are conserved in time (dynamic equilibrium); chosen to represent a range of productivity and species rich- this meant that the relationship between plant diversity and ness (plot3: Nil/unfertilized, plot7/2: ‘P K Na Mg’, plot6: yield variance could be included in our analysis (Silvertown −1 ‘N1 P K Na Mg’ with 48 kg N  ha ammonium sulfate, et al. 2006). −1 plot 9/2: ‘N2 P K Na Mg’ with 96 kg N  ha ammonium −1 sulfate, plot11/1: ‘N3 P K Na Mg Si’ with 144 kg N  ha 2.2 Environmental abiotic covariates ammonium sulfate, plot13/2: manure applied in a 4-year −1 cycle, plot17: N*1 with 48 kg N  ha sodium nitrate), except Based on statistical analyses (see Section  2.3), the effect on plot 6, where only the a and b sub-plots were available. of the following environmental abiotic covariates on tem- The plots 7/2, 6, 9/2, 11/1, and 17 represent different lev - poral yield variance was tested. Different climatic covari - els and forms of inorganic fertilization, whereas the plot ables (air temperature, humidity, rainfall, hours of sun- 13/2 represents organic fertilization (poultry manure since shine, soil moisture, radiation, wind) were measured daily 2003, before farmyard manure). Additional treatment details at Rothamsted Research (available at: http:// www. era. rotha are provided in Table A3 Supplementary material. Yield msted. ac. uk/# measu remen ts), except for measurements of data were obtained from the electronic Rothamsted Archive atmospheric carbon dioxide concentration (Fig. A4 Sup- ‘e-RA’ and have been made publicly available on the e-RA plementary material) that was recorded monthly. In addi- website ( https:// doi. org/ 10. 23637/ r pg5- yield s1965- 2018- tion, air chemistry parameters were measured, N deposition, 01) (Perryman and Ostler 2021) and included plot-specific and SO emissions (1965–2018, UK), which are shown in twice-yearly yields (total aboveground biomass; sum of 1st Fig. A5 Supplementary material. Besides these measured and 2nd cuts; 100% dry matter). The first cut was made into covariates, the accumulated numbers of water stress days for hay and removed in mid-June, and the second cut was taken the vegetation periods 1st March–15th of June (time period with a forage harvester while still green (end-October, bio- up to typical date for 1st cut/harvest), 16th June–31st Octo- mass removed). ber (time period up to typical date for 2nd cut/harvest), and The plant communities on the PGE are naturally assem- for the entire period (Fig. 2) were calculated yearly in the bled from a regional species pool that is classified as dicot- following manner. The amount of maximum plant available yledon-rich Cynosurus cristatus–Centaurea nigra grassland, water (MPAW [mm]), which acted as a ‘bucket of water’, one of the mesotrophic grassland communities in the British was specified for the soil. At the start of the calculation (1 National Vegetation Classification system (Dodd et al. 1994). Jan 1965), it is assumed that the actual plant available water The botanical composition of all sub-plots was studied annu- (APAW(t)) equaled the MPAW (‘the bucket was full’). The ally from 1991 to 2000 and from 2010 to 2012 by recording daily water surplus was calculated as rainfall minus potential the dry mass of each plant species in early June (number evaporation over grass (WS(t)), with potential evaporation of species, Shannon’s diversity index) (Storkey et al. 2015; over grass as a derived meteorological variable (see formula Crawley et  al. 2005). The species of grasses, forbs, and details: http:// www . er a. r o t ha ms ted. ac. uk/ info/ me t/ der iv legumes comprising at least 5% of the aboveground bio-ed_ varia bles# EVAPG). A daily balance was calculated as mass found in these surveys and used for this study have DB(t) = APAW(t − 1) + WS(t). For positive WS(t)-values, been made publicly available on the e-RA website (https:// APAW increased to the maximum of MPAW. Above this doi. org/ 10. 23637/ rpg5- speci es_ 1991- 2000- 01) (Perryman value, water was assumed to run off or drain away. For nega - et al. 2021). An overview about the changes in the number tive values, APAW was reduced, and when reaching 0, no of plant species over time is provided in Fig. A11 Supple- more water extraction was possible. mentary material; a more detailed description of the PGE biodiversity data has been reported in the Rothamsted Guide If DB(t) > MPAW; APAW(t) = MPAW to Classical Experiments on pp. 25–27 (Macdonald et al. If MPAW > DB(t) > 0; APAW(t) = DB(t) 2018) and is available online: http:// www. era. rotha msted. If DB(t) < 0; APAW(t) = 0 ac.uk/ home/ W eb_L TE_Guide book_ 2018_ 2019- r eprint. pdf . Because of the resource required for vegetation assessments, Days where APAW was 0 were counted as stress days. data on plant diversity are only available for a sub-set of the The numbers of water stress days accumulated for the 1st years and plots used in this analysis. However, the data cover March–15th of June (harvest time of the 1st cut), 16th a significant proportion of the time period analysed in our June–31st of October (harvest time of the 2nd cut), and for study that includes large interannual variation in weather the entire period each year were determined. No water stress and yields and a period of change in atmospheric chemistry occurred before March or after October. Several values of 1 3 Long-term trends in yield variance of temperate managed grassland Page 5 of 19 37 Fig. 2 Accumulated number of water stress days for the vegeta- tion period from March to Octo- ber and temporal development of the mean air temperature (°C) for the months July–August at Rothamsted (1965–2018). Water stress was defined as a limited plant available soil water content (formula described in “Material and methods” sec- tion). Further information about the temperature anomaly is provided in Fig. A2 Supplemen- tary material. MPAW were tested, but an MPAW of 135 mm provided the spaced knots (set to 10, approximately 5-year intervals). For best correlation for both the 2nd cut and the yearly yield, in each of these 10 sub-periods, a separate residual variance line with the soil description (Kohler et al. 2016). was fitted to assess changes in temporal yield variance, with lower values indicating less unexplained variability between 2.3 Statistical analysis years (=more stable yields). Each ‘fertilizer × liming’ treat- ment combination was assumed to have a specific set of To account for the experimental design, a mixed model variance components ([t/ha] ) for the different periods. We was used based on REML, as recommended by Raman denote the period-specific yield variances as ‘environmental et al. (Raman et al. 2011) and Onofri et al. (Onofri et al. variance’, a term coined by Römer (1917). Here, we replace 2016). Each plot had a different ‘fertilizer × liming’ treat- the arithmetic treatment mean with the spline estimate of ment combination, with no replication or randomization the temporal trend. The model, therefore, allowed yield vari- (Fig. A1 Supplementary material), so that plot errors and ance to be determined independently of yield level, which ‘fertilizer × liming’ interactions (=residual) could not be highlights that the statistical approach used here differs from separated. The size of the main plots (306–1912  m ) might earlier approaches based on the PGE, which used the clas- partly compensate for the lack of replication, particularly sic coefficient of variation for a measure for yield variabil- because the experimental site was reasonably uniform when ity (Dodd et al. 1997). This is an important improvement the experiment started in 1856 (Crawley et al. 2005; Lawes because yield variance can be incorrectly interpreted if there and Gilbert 1859). is a systematic dependency of the measure of variation on To answer the first research question comparing tem- the mean. No such dependencies were found in this analysis. poral trends in mean yield and yield variance across fer- The model was fitted separately for each ‘fertilizer × liming’ tilizer × liming treatments, we fitted for each plot (‘ferti- treatment. In this analysis, any systematic year effects (trend, lizer × liming’ treatment combination) a smoothing spline etc.) are captured by the fitted splines, whereas the random for the mean yield trend via a random-effects specification residual captures any unexplained year effects. as a mixed model (Verbyla et al. 1999). An advantage of the To answer the second research question of how much of smoothing spline approach over other regression approaches the variation in yield can be explained by climatic drivers is that no specific assumption is needed as regards the func- and the interaction with agronomic management, treatment- tional form of the trend. Preliminary inspection of the data specific multiple regression analyses were performed. The revealed that such flexibility was needed. The entire model strength of the relationship between the response vari- syntax (using R version 4.0.0) is provided in Table A6 Sup- able ‘yield’ and several explanatory environmental abiotic plementary material. The trends were fitted using ASReml- covariates (as described in former Chapter 2.2), as well as R using the specification random = ~spl(t ,k = 10), where t the importance of each of the predictors to the relationship, is continuous time in years and k is the number of evenly was assessed. We further tested for significant differences 1 3 37 Page 6 of 19 J. Macholdt et al. in regression coefficients of the main effects between treat- covariates. With regard to the nonnormal distributions of ments by using a general linear model approach, but no sig- the underlying data (outliers), correlation coefficients were nificant differences between treatments were found. Further - calculated from the ranks of the data, not from their actual more, a novel criss-cross regression approach with extended values. Kendall’s tau (τ) was used to ensure that the results Eberhart–Russell regression analyses (Eberhart and Russell were accurate because the same ranks were repeated too 1966) was used, in which the environmental mean was mod- many times in the partial datasets (e.g. botanical surveys). eled using the previously as significant selected environ - The strength of the correlation increased both from 0 to +1 mental abiotic covariates ((1) accumulated days of water and 0 to −1, where −1/+1 indicated the strongest correla- stress from March to October; (2) mean air temperature tion and 0 indicated no correlation. The sign of τ showed the from May to June; (3) mean air temperature from July to direction of the correlation; if negative, the variables were August based on the multiple regression analyses; shown inversely related. To gain further insight into the temporal in Table  1). This new method allowed us to assess the dynamics of the plant communities, a multivariate analy- treatment-specific yield sensitivity (or responsiveness) to sis was done using the data on relative biomass recorded variability in these climatic conditions. The Eberhart–Rus- on all study plots for the 10-year period (1991–2000) for sell model was further fitted using the criss-cross regres- which data were available. Although this only covered part sion approach proposed for Finlay–Wilkinson regression by of the total period covered by the yield variance analysis, Digby (2009) and extended for the mixed model version it included a large range of yields. Two partial canonical by Nabugoomu et al. (1999). This iterative scheme allowed correspondence analyses (pCCA) were done after remov- for (i) regression of the environmental index on covariates, ing species that were recorded in less than 5% of samples to (ii) heterogeneity of variance for the independent devia- avoid bias owing to rare species. First, the effect of year was tions from the regression, and (iii) serial correlation of the analysed, including plot as a covariate. Second, the variance residuals. The intercept of Finlay–Wilkinson regressions in community composition explained by plot was quanti- provided information about the general yield level of a treat- fied, including year as a factorial covariate. In both cases, ment compared to others. The slopes of regression lines can the proportion of functional group (grass/forb/legume) was be interpreted as ‘yield sensitivity’ to climatic perturbation included in the ordination plots as a supplementary variable. (slope > 1: higher sensitivity; slope < 1: less sensitivity/ better resilience). Treatments with a slope of approximately 1 showed an average yield reaction across the treatments, 3 Results and discussion similar to the average response indicated by the environ- mental mean (reference, black regression line with a slope 3.1 Temporal trends in yield variance of 1; Fig. 5a–d). We would like to stress that this criss-cross and identifying main climatic drivers regression approach for the extended Finlay–Wilkinson regression was specifically developed for the PGE and, to the Based on the unique long-term dataset and the new designed best of our knowledge, constitutes a novel method. Briefly, statistical approach, this study provides, for the first-time, the Finlay–Wilkinson regression model can be written as insights into trends in temporal yield variance of grassland =  +  w , where  is the expected performance of the in response to climatic drivers and its dependence on agro- ij i i j ij ith treatment in the jth environment,  and  are intercept nomic management; yield variance was determined indepen- i i and slope (‘yield sensitivity’) for the ith genotype, and w is dently of yield level, which is important to avoid any misin- the environmental mean of the jth environment. The envi- terpretation of resilience. Our work adds a new perspective ronmental mean, in turn, is modelled by a linear regression to earlier productivity analyses of grassland experiments, as w =  +  x +  x + .... +  x , where x is the value which were based mostly on shorter time periods (Prieto j 0 1 1j 2 2j p pj hj of the hth covariate in the jth environment and  ,… ,  are et al. 2015; Tilman et al. 2001; Graux et al. 2013; Trnka 0 p regression parameters. A detailed description of the model et al. 2021; Storkey et al. 2015; Haughey et al. 2018; Sander- and its estimation is provided in Table A7 Supplementary son 2010), and reveals the importance of long-term experi- material. ments for detecting possible trends over time. It is comple- Regarding the third research question, we aimed to iden- mentary to a recent analysis of trends in productivity on the tify the potential correlation between plant species diversity PGE that showed a consistent decline in yields in response to (number of species and Shannon’s diversity index) and the climate change across a longer time period (1902–2016) for mean yield or temporal yield variance depending on the four treatments and using yield data from just the first hay agronomic treatment (fertilizer × liming). For this correla- cut (Addy et al. 2022). By focusing on yield variance (using tion analysis, we used the mean yield and temporal yield data from both cuts) across a wider range of treatments, our variance data based on the statistical analyses as described study provides additional insights into the potential for man- above to account for the experimental design and other agement to impact adaptability and resilience of grasslands 1 3 Long-term trends in yield variance of temperate managed grassland Page 7 of 19 37 1 3 Table 1 Treatment-specific multiple regression analyses for the dependent variable ‘yield’ and selected climatic drivers. Analyses based on data (1965–2018) of total yield (sum of 1st + 2nd cut). Abbreviations: VIF = variance influence factor (measures the strength of the correlation between the independent variables in regression analysis); regression coefficient B = slope ‘b’ of regression line (how much Y changes for each change in X); standardized coefficient Beta = analogous to the interpretation of ‘b’, except that Beta expresses change in standard scores. Climatic drivers with significant main effects were selected based on a preceding general linear model analysis: accumulated number of water stress days from March to October (limited plant available soil water content with expected drought stress) and mean air temperature (=Temp.) in summer months (Fig. 2). No significant differences in regression coefficients B between treatments were found. Treatment explanations are provided in Table A3 Supplementary material (FYM/PM = farmyard/poultry manure; N* = sodium nitrate). Treatment Fertilization Liming General model statistic Second selected covariable First selected covariable Durbin Adjusted VIF Covariable Regression Standardized p value Covariable Regression Standardized p value Watson R Square Collinearity coefficient coefficient coefficient coefficient Beta test B Beta B 3 Nil pH 7 2.02 0.29 1.00 Water stress -0.02 -0.55 0.00 none selected n/a n/a n/a pH 6 1.68 0.23 1.00 Water stress -0.02 -0.49 0.00 none selected n/a n/a n/a pH 5 1.62 0.29 1.00 Water stress -0.02 -0.54 0.00 none selected n/a n/a n/a no chalk 1.58 0.25 1.00 Water stress -0.02 -0.50 0.00 none selected n/a n/a n/a 7/2 P K Na Mg pH 7 1.97 0.49 1.36 Water stress -0.03 -0.47 0.00 Temp.July-Aug. -0.66 -0.35 0 00 pH 6 2.05 0.37 1.36 Water stress -0.02 -0.37 0.01 Temp.July-Aug. -0.63 -0.35 0 01 pH 5 1.76 0.34 1.01 Water stress -0.04 -0.52 0.00 Temp.Mar.-Apr. 0.96 0.35 0 00 no chalk 1.77 0.25 1.00 Water stress -0.03 -0.52 0.00 none selected n/a n/a n/a 6 N1 P K Na Mg pH 7 2.05 0.48 1.29 Water stress -0.02 -0.42 0.00 Temp.July-Aug. -0.76 -0.40 0 00 pH 6 2.06 0.40 1.29 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.70 -0.37 0.01 pH 5 n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a no chalk n/a n/a n/a  n/a n/a n/a n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 1.71 0.36 1.36 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.79 -0.31 0.02 pH 6 1.51 0.36 1.36 Water stress -0.02 -0.39 0.01 Temp.July-Aug. -0.64 -0.38 0.01 pH 5 1.70 0.35 1.36 Water stress -0.02 -0.36 0.01 Temp.July-Aug. -0.74 -0.35 0.01 no chalk 1.57 0.38 1.36 Water stress -0.02 -0.40 0.00 Temp.July-Aug. -0.85 -0.32 0.02 11/1 N3 P K Na Mg pH 7 1.64 0.29 136 Water stress -0.02 -0.32 0.02 Temp.July-Aug. -0.73 -0.33 0.02 pH 6 1.62 0.27 136 Water stress -0.02 -0.29 0.04 Temp.July-Aug. -0.67 -0.34 0.02 pH 5 1.71 0.26 136 Temp.July-Aug. -1.30 -0.52 0.00 Waterstress -0.02 -0.25 0.05 no chalk 1.82 0.35 136 Temp.July-Aug. -1.47 -0.60 0.00 Waterstress -0.01 -0.16 0.05 13/2 FYM/PM pH 7 1.72 0.34 136 Temp.July-Aug. -0.91 -0.40 0.00 Waterstress -0.02 -0.29 0.03 pH 6 2.16 0.38 136 Temp.July-Aug. -0.91 -0.44 0.00 Waterstress -0.02 -0.29 0.03 pH 5 1.72 0.33 136 Water stress -0.03 -0.39 0.01 Temp.July-Aug. -0.67 -0.30 0.03 no chalk 1.64 0.21 1.00 Water stress -0.03 -0.47 0.00 none selected n/a n/a n/a 37 Page 8 of 19 J. Macholdt et al. to abiotic stress. Our treatment-specific model allowed tem- poral trends, using a spline component, to be dissected from residual variance, fitted as separate parameters for ten con- secutive sub-periods. As a striking finding, we identified a similar pattern of temporal trends in yield variance (overview results Fig. 3) (Römer 1917), particularly for plots with low rates of N −1 (48 kg N  ha , N1) or no N fertilizer ‘Nil’ (plot specific results Fig. 4). This pattern can be described as follows: in the 1960s and 1970s, the plots showed relatively low yield variance. In subsequent years, we observed greater yield variance with a peak at approximately 1995 before becom- ing more stable in the past two decades (Fig. 3). This pattern of temporal trends in yield variance appeared to be less pro- nounced in treatments with a high soil pH; in contrast, peaks were more pronounced in treatments with low pH (Fig. A8 Supplementary material). The exceptions were treatments with reduced or no liming, but with a high mineral N supply (Fig. 4e) or farmyard manure (FYM) (Fig. 4g), where these patterns were not as evident, and periods of high yield vari- ability were observed throughout the study period. Key message for ‘Research question I’: There were con- sistent unimodal trends observed in yield variance in plots with low to moderate or no nitrogen fertilizer additions, with a peak in 1990s, after which variability declined. Yield was most variable in plots with higher nutrient inputs and lower soil pH. Based on treatment-specific multiple regression analy - ses including environmental abiotic covariates, we explored the possible causes of the observed temporal trends in yield variance. We identified the accumulated days of water stress from March to October and the mean air temperature from July to August as the two main climatic drivers, explain- ing around one-third of the observed yield variance across treatments, or even up to 49%/48% of yield variance in treat- ment P K Na Mg and N1 P K Na Mg, respectively (Table 1 and Fig. 2). Overall, the impact of temperature driving yield variance was lower than the stronger impact of water stress (Table 1: see columns ‘standardized coefficient Beta’ with lower absolute values indicating less impact). The identi- fied main climatic drivers tally with previous analyses of the PGE made by Dodd et al. (1994, 1997) and recent ones by Addy et al. (2021, 2022); any differences can be explained by the fact that we included the second cut in the yield data. The seasonality of temperate grassland production is primar- ily affected by soil moisture and temperature, which con- strain the length and determine the intensity of the growing season (Trnka et al. 2011). Generally, significant tempera- ture changes, particularly hot temperatures in summer with a limiting water balance, have a negative effect on grassland productivity (Kipling et al. 2016; Höglind et al. 2013) and are expected to increase the interannual and seasonal pro- duction variability of grassland systems (Graux et al. 2013; 1 3 Table 1 (continued) Treatment Fertilization Liming General model statistic Second selected covariable First selected covariable Durbin Adjusted VIF Covariable Regression Standardized p value Covariable Regression Standardized p value Watson R Square Collinearity coefficient coefficient coefficient coefficient Beta test B Beta B 17 N*1 pH 7 1.95 0.31 1.19 Water stress -0.02 -0A2 0.00 Temp.May-June -0.25 -0.27 0.04 pH 6 1.91 0.23 136 Water stress -0.01 -0.29 0.04 Temp.May-June -0.24 -0.29 0.05 pH 5 1.80 0.32 136 Temp.May-June -0.25 -0.58 0.00 Waterstress -0.02 -0.34 0.01 no chalk 1.73 0.28 1.00 Temp.July-Aug. -0.58 -0.46 0.00 Waterstress -0.01 -0.22 0.05 Long-term trends in yield variance of temperate managed grassland Page 9 of 19 37 Fig. 3 Summary plot for the overview temporal trend in mean yield (blue dotted line; grey crosses real harvest data) and yield variance (red bars) including standard errors (grey error bars) based on the total mean over all liming × fertiliza- tion treatments (1965–2018). Underlying plot specific results are provided in Fig. 4; overview of liming treatments (Fig. A8 Supplementary material), and overview of fertilization treat- ments (Fig. A9 Supplementary material). Detailed informa- tion about treatments is shown in Table A3 Supplementary material. Chang et al. 2017). However, the water stress index we cal- this study). A recent study by Addy et al. (2021) identified culated is partly based on a derived meteorological (potential clusters of years with similar weather patterns between 1900 evaporation from grass) and an estimated parameter (soil and 2020 at Rothamsted, which might help to clarify the water storage), which may differ from actual measurements. unexplained rest of unimodal yield variance. They found a This may have contributed to the lower correlations seen climate cluster characterized by cool and dry springs from in the Nil (no fertilizer) and N*1 plots together with the approximately the 1960 to 1970s (in our study: stable yields fact that other nutrients may have limited growth. Further- were observed at the beginning), followed by a period with more, water stress often occurs in July–August, so the fact a variety of clusters and widely varying weather patterns that both factors come out as important may indicate that (in our study: more variable yields occurred from approxi- plants may be hit harder by high temperatures when transpi- mately 1980 to 2000), and a transition since 2000 with an ration and thereby cooling are limited. The effect of water increased tendency toward higher temperatures in springs stress was consistently present in all treatments and can be and drier periods in June (in our study: associated with more assumed to be the main driving factor in the PGE. Regard- stable yields at the end; Fig.  3). The study suggests that ing the effect of agronomic practices on resistance to water the positive effects of sufficient water availability can off- stress (Table 1 and Fig. 2), sufficient liming is important for set the negative effects of warmer temperatures on pasture sustaining grassland productivity due to positive effects on performance (Addy et  al. 2021). An additional explana- soil pH, root growth, soil structure, nutrient availability, soil tion could be that this peak in yield variance between 1980 carbon, and soil biota (Holland et al. 2018; Fornara et al. and 2000 coincides with a period of decreasing nitrogen 2011). deposition and SO emissions (see Fig. A5 Supplementary Key message for ‘Research question II—part A’: The material), which is reflected in shifts in plant community accumulated days of water stress from March to October composition, especially an increase in the relative propor- and mean air temperature from July to August were the most tion of legumes. Regarding the SO emissions, there was a −1 −1 important climatic drivers, explaining approximately one- decline from approximately 65 kg  ha in 1980 to 5 kg  ha third of the observed interannual yield variance. in 2006 (Anon 2006). We might expect this to affect species The remaining unexplained yield variance in the PGE diversity only in S-limited plots, but S is applied together might be driven by different environmental factors, which with K, Na, and Mg in most plots (such as in the ‘P K Na could not be determined further (i.e. due to lack of data in Mg’ treatment) in the PGE, so any effect of changes in air 1 3 37 Page 10 of 19 J. Macholdt et al. Trend in mean yield Temporal yield variance [(t/ha)²] [t/ha] 15 12 A) Nil Trend in mean yield Yield variance 10 8 5 5 4 0 0 15 12 B) P K Na Mg - pH 7 10 8 5 4 0 0 15 12 C) N1 P K Na Mg 5 4 0 0 15 12 D) N2 P K Na Mg 10 8 5 4 0 0 15 12 E) N3 P K Na Mg 10 8 5 4 0 0 15 12 F) N*1 10 8 5 4 0 0 15 12 G) FYM/PM 10 8 5 4 0 0 1965 1985 2005 1965 1985 2005 1965 1985 2005 1965 1985 2005 pH = 7 pH = 6 pH = 5 no chalk 1 3 Long-term trends in yield variance of temperate managed grassland Page 11 of 19 37 ◂Fig. 4 Temporal trends in mean yield (blue splines with approxi- and transpiration, as affected by increasing temperature mated confidence intervals) and temporal yield variance (red bars (Fig. A2 Supplementary material) and atmospheric C O con- incl. standard errors) depending on the treatment (fertilizer  ×  lim- centrations (C ) (Fig. A4 Supplementary material). It was ing) for the experimental period of 1965–2018. Analysis based on noted that grasses appeared to be more sensitive to increas- year  ×  plot specific yields; these raw yield data are shown as black dots in each graphic. A Treatment no. 3: Nil (no fertilizer input)—pH ing C than forbs, resulting in lower water use efficiency, 7/6/5/no chalk. B Treatment no. 7/2: P K Na Mg—pH 7/6/5/no chalk. decreased N uptake, and declining biomass production in C Treatment no. 6: N1 P K Na Mg—pH 7/6 (restricted data availabil- grass-rich communities (Baca Cabrera et al. 2021). This indi- ity: only 1972–2018 and pH 7/6). D Treatment no. 9/2: N2 P K Na cated that plant communities comprising only a few species Mg—pH 7/6/5/no chalk. E Treatment no. 11/1: N3 P K Na Mg—pH 7/6/5/no chalk. F Treatment no. 17: N*1 (N* = sodium nitrate)—pH may be disproportionately affected if the species present are 7/6/5/no chalk. G Treatment no. 13/2: FYM/PM (farmyard/poultry poorly adapted to changing climatic conditions and are also manure)—pH 7/6/5/no chalk. Yield variance denoted Römer’s envi- less adaptable to long-term climatic trends (Eisenhauer et al. ronmental variance, with lower values indicating more stable yields 2019). Our results on trends in yield variance and mean yield and higher values indicating more variable yields. Detailed informa- tion about treatments is shown in Table A3 Supplementary material. support this hypothesis (Table 2 and Fig. 4). Summary plots are provided as an overall overview (Fig. 3), overview of liming treatments (Fig. A8 Supplementary material), and overview 3.2 Determine yield sensitivity to climatic changes of fertilization treatments (Fig. A9 Supplementary material). using a novel criss‑cross regression analysis chemistry is more likely a result of decreasing N deposition We present a novel criss-cross regression analysis (extended in these plots. The highest peak in yield variance in the mid- Finlay–Wilkinson regression), in which we modeled the dle period was observed in plots lacking N but supplied with environmental mean using the three main identified cli - P, K, Na, and Mg (Fig. A9). These plots have the highest matic factors (Table 3 and Fig. 5). This approach allowed proportion of legumes (Table 2) that are sensitive to changes us to determine the treatment-specific yield sensitivity to in N deposition. This suggests that where yield is dependent variation in climatic conditions. A great advantage of this on biological N fixation (see key message IV), resilience approach is its parsimony, providing a single slope for each will also be determined by the specific response of legumes treatment that assesses sensitivity to all included covariates to environment change—such as atmospheric N deposition simultaneously. The interaction between fitted environmen - (Storkey et al. 2015). tal mean and treatments was highly significant (F  = 3.17; In addition to yield variance, the temporal trends in mean P < 0.0001; Table A7 Supplementary material), showing yield were in some cases relatively stable, as for the unfer- that the slopes are different between treatments. The results tilized ‘Nil’ treatment with a relatively constant yield of shown in Table 3 can be interpreted as follows: a lower slope −1 approximately 3 t  ha 286 (Fig.  4a). In other cases, the indicates less sensitivity (or responsiveness) of yield to vari- mean yield decreased over time (negative trends), supporting ation in climatic conditions; a lower variance of deviation the conclusions of Addy et al. (2022), who analysed PGE indicates a more stable yield. A higher slope indicates a data over a longer time period (1902–2016) and modelled higher sensitivity of yield to changing climatic conditions; spring hay yields under four fertilizer regimes in response a higher variance of deviation indicates more yield fluc- to seasonal temperature and rainfall. The PGE modelling tuations. Regarding the predicted environmental mean, the study showed that warmer and drier years in the twentieth signs of all three regression coefficients (theta1–theta3) and twenty-first centuries resulted in yield reductions and were negative, meaning the predicted environmental mean are forecasted to decline further up to 50% under future increases with decreasing values for x1, x2, and x3. Thus, (2020–2080) climate scenarios (Addy et al. 2022). less water stress (x1) and lower temperatures from May to This was the case particularly for treatments with greater June (x2) and from July to August (x3) resulted in higher inputs of N and those provided with organic manures, includ- grassland yields during the observed experimental period. ing the ‘N3 P K Na Mg—pH 7 or 6’ (Fig. 4e) and the ‘FYM/ Key message for ‘Research question II—part B’: Liming PM—pH 7 or 6’ (Fig. 4g) treatments in which mean yields to attain a soil pH of 6–7 and moderate N supply (e.g. treat- −1 −1 decreased from >9 t  ha to nearly 5 t  ha from 1965 to ment ‘N2 P K Na Mg—pH 6/7) were identified as the most 2018. In the PGE, the higher input plots (see ‘N3 P K Na promising agronomic practices for sustaining yield under Mg—pH 5 or no chalk’) were dominated by only a few grass varying climatic conditions and for reducing yield sensitivity species (Table 2), which might have resulted in greater yield to abiotic stresses. variance when the climatic conditions were unfavourable for In the intensively fertilized treatment without liming ‘N3 these species and reduced adaptability. Baca Cabrera et al. P K Na Mg—no chalk’, the regression lines had a slope >1 (2021) reported that the declining yields observed in grass- (b = 1.32), which indicates a high yield sensitivity to abiotic rich plant communities of the PGE over the last century were stress—where water is limiting or there is heat stress in the associated with decreases in N uptake, stomatal conductance, summer, proportionally more of the potential yield is lost 1 3 37 Page 12 of 19 J. Macholdt et al. Table 2 Overview results for treatment-specific mean yield, tempo- 2010–2012; species number: 1974, 1991–2000, 2010–2012; propor- ral yield variance, and plant species diversity.  n/a: no data available. tion legumes: 1991–2000). Treatment explanations are provided in Mean values based on available years of surveys (soil pH: 1998– Table  A3 Supplementary material (FYM/PM  =  farmyard/poultry 2014; yield: 1965–2018; Shannon’s diversity index: 1991–2000, manure; N* = sodium nitrate). Treatment Fertilization Liming Soil pH Mean yield [t/ha] Yield variance [(t/ha) ] Shannon's Species Proportion diversity number legumes index [%] Parameter Standard error Parameter Standard error estimate estimate error 3 Nil pH 7 7.2 3.14 0.16 1.11 0.21 2.6 31 5.8 pH 6 6.3 3.52 0.18 1.51 0.35 2.7 30 3.6 pH 5 5.1 2.12 0.23 0.93 0.26 2.1 27 1.1 no chalk 5.2 2.67 0.29 1.40 0.12 2.2 26 0.3 7/2 P K Na Mg pH 7 7.0 7.02 0.26 2.10 0.17 2.4 23 28.6 pH 6 6.2 7.46 0.16 1.97 0.15 2.3 23 16.5 pH 5 5.2 5.94 0.35 2.59 0.19 2.3 22 25.7 no chalk 5.0 4.15 0.32 2.63 0.16 2.1 23 17.4 6 N1 P K Na Mg pH 7 7.0 7.16 0.24 1.94 0.41 2.5 22 17.1 pH 6 5.9 7.06 0.24 2.01 0.42 2.0 20 11.3 pH 5 n/a n/a n/a n/a n/a n/a n/a n/a no chalk n/a n/a n/a n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 7.1 7.32 0.26 1.62 0.21 2.1 17 2.7 pH 6 6.2 7.37 0.26 1.82 0.19 1.8 17 1.6 pH 5 5.0 6.45 0.22 2.69 0.22 1.5 15 10.8 no chalk 3.6 5.39 0.35 2.38 0.29 0.6 3 0.0 11/1 N3 P K Na Mg pH 7 7.0 8.45 0.38 2.18 0.31 1.5 11 0.0 pH 6 6.1 7.54 0.30 1.99 0.19 1.6 12 0.0 pH 5 5.1 7.16 0.36 2.75 0.27 0.8 10 0.0 no chalk 3.5 6.54 0.27 3.38 0.41 0.0 2 0.0 13/2 FYM/PM pH 7 6.9 7.07 0.21 1.87 0.18 2.3 20 0.6 pH 6 6.0 7.75 0.19 2.27 0.23 2.3 20 4.6 pH 5 5.3 6.86 0.35 2.79 0.39 2.3 20 0.2 no chalk 5.1 6.17 0.40 2.85 0.27 2.0 19 0.5 17 N*1 pH 7 7.1 3.95 0.22 1.07 0.25 2.2 23 1.6 pH 6 6.3 4.03 0.18 1.25 0.22 2.1 24 0.1 pH 5 5.8 4.11 0.14 1.24 0.17 2.0 22 0.0 no chalk 5.8 3.92 0.12 0.82 0.11 2.2 24 0.0 (Table 3 and Fig. 5d). In comparison, for treatment ‘N2 P K In addition to slopes, the intercept of regression lines is the Na Mg—pH 6/7’, lower slopes (b = 0.95/0.99) of regression second important criterion of Finlay–Wilkinson regression lines were found, suggesting less yield sensitivity (or respon- analysis. The intercept provides information about the overall siveness) to climatic perturbation (Table 3 and Fig. 5a, b). yield level of a treatment (Table 3 and Fig. 5), with a higher This confirms findings of short-term grassland studies show - intercept indicating higher yield performance (e.g. treatment ing that higher/intensive nutrient inputs often increase yield ‘N3 P K Na Mg—pH 6/7’) and a smaller (or more negative) but can destabilize productivity (Hautier et al. 2014, 2020; value referring to treatments with low yield level (e.g. unfer- Zhang et al. 2016, 2019; Crawley et al. 2005). Moreover, the tilized treatment ‘Nil—no chalk’). Overall, the combined results showed that plots with greater fertilizer inputs were evaluation of slope, intercept, and variance of deviation can less resistant to climatic perturbation (higher yield variance), point to treatments that show a favourable combination of which is in agreement with the findings of Deguines et al. resistance to climatic perturbation (slope < 1) together with (2014), who reported that relative adaptability decreased a high and stable yield level (high intercept, low variance of with increasing land use intensity. deviation). In this study, liming (pH 7/6) and moderate N 1 3 Long-term trends in yield variance of temperate managed grassland Page 13 of 19 37 Table 3 Treatment-specific criss-cross regression analyses (see predicted mean yields for each year (=environmental mean, x-axis) Fig.  5; detailed explanation Table  A7 Supplementary material).  The and regressed on the fertilization  ×  liming treatment yields (y-axis); calculation of the environmental mean considered the three selected getting regression lines with treatment-specific intercepts and slopes. climatic drivers (x1 ‘accumulated number of water stress days’; x2 Treatment explanations are provided in Table  A3  Supplementary ‘mean air temperature from May to June’; x3 ‘mean air temperature material (FYM/PM = farmyard/poultry manure; N* = sodium nitrate). from July to August’ (Fig.  2 and Table  1), and were used to obtain Treatment Fertilization Liming Criss-cross regression analyses Slope Intercept Variance of deviation Parameter Standard error Parameter Standard error Parameter estimate estimate estimate 3 Nil pH 7 0.70 0.14 -1.21 0.96 0.04 pH 6 0.76 0.14 -1.18 0.95 0.03 pH 5 0.69 0.15 -1.93 1.01 0.12 no chalk 0.88 0.16 -2.74 1.06 0.20 7/2 P K Na Mg pH 7 1.18 0.17 -0.21 1.11 0.29 pH 6 1.99 0.17 1.14 1.13 0.33 pH 5 1.17 0.18 -1.16 1.16 0.39 no chalk 1.02 0.17 -1.60 1.14 0.36 6 N1 P K Na Mg pH 7 1.18 0.17 -0.11 1.12 0.30 pH 6 1.07 0.18 0.38 1.14 0.33 pH 5 n/a n/a n/a n/a n/a no chalk n/a n/a n/a n/a n/a 9/2 N2 P K Na Mg pH 7 0.93 0.17 1.86 1.14 0.34 pH 6 0.98 0.18 1.55 1.17 0.42 pH 5 1.21 0.17 -1.17 1.11 0.29 no chalk 1.19 0.19 -1.90 1.27 0.62 11/1 N3 P K Na Mg pH 7 1.02 0.21 2.21 1.35 0.83 pH 6 0.88 0.17 2.41 1.14 0.36 pH 5 1.09 0.20 0.76 1.29 0.69 no chalk 1.36 0.25 -1.79 1.60 1.50 13/2 FYM/PM pH 7 1.11 0.21 -0.01 1.37 0.87 pH 6 1.19 0.22 0.31 1.43 1.02 pH 5 1.20 0.20 -0.33 1.32 0.76 no chalk 1.09 0.21 -0.32 1.36 0.83 17 N*1 pH 7 0.86 0.15 -1.41 0.98 0.08 pH 6 0.77 0.15 -0.67 0.97 0.05 pH 5 0.87 0.17 -1.20 1.09 0.25 no chalk 0.62 0.16 0.19 1.07 0.22 Predicited x1 (Water stress days) -0.02 0.00 14.66 1.45 n/a environmental x2 (Temperature May-June) -0.18 0.13 14.66 1.45 n/a mean x3 (Temperature July-August) -0.64 0.11 14.66 1.45 n/a supply (e.g. treatment ‘N2 P K Na Mg—pH 6/7’) were identi- suggesting a similar environmental adaptability and yield fied as the most important agronomic practices for sustaining performance. The largest yield susceptibility to lower pH yield under changing climatic conditions (Table 3 and Fig. 5). values was observed in the ‘P K Na Mg’ treatment (possibly A striking finding was that the position and slope of explained by the specific response for legumes, see above), regression lines, except ‘Nil’ (unfertilized) and ‘N*1’ (with and the least reactivity was shown by the ‘N3 P K Na Mg’ −1 48 kg N  ha sodium nitrate), showed the least differentiation treatment, followed by the ‘N2 P K Na Mg’ and FYM/PM’ between treatments in the liming variant ‘pH 6’ (Fig. 5b), treatments (Table 3 and Fig. 5). 1 3 37 Page 14 of 19 J. Macholdt et al. Fig. 5 Graphical visualization Treatment specific yield of the treatment-specific criss- [t/ha] cross regression analyses (see Table 2; detailed explanation (a) pH 7 (b) pH 6 provided in Table A7 Supple- mentary material) depending on liming. a Treatments with Environmental pH 7. b Treatments with pH mean 6. c Treatments with pH 5. d N3PKNaMg Treatments with no chalk. The N2PKNaMg calculation of the environmen- N1PKNaMg tal mean considered the three PKNaMg selected climatic drivers: the 4 FYM/PM ‘accumulated number of water stress days’ and ‘mean air N*1 temperature from May–June Nil and July–August’ (Fig. 2 and Table 1) were used to obtain predicted mean yields for each year (=environmental mean, x-axis) and regressed on the fertilization × liming treat- (c) pH 5 (d) no chalk ment yields (y-axis); getting regression lines with treatment- specific intercepts and slopes. Treatment explanations are provided in Table A3 Sup- plementary material (FYM/ PM = farmyard/poultry manure; N* = sodium nitrate). 0246 810 0246 810 Environmental mean yield [t/ha] Key message for ‘Research question III—part A’: Higher 3.3 Correlation between plant species diversity plant species diversity was correlated not only with more and yield variance stable grassland yields, but also with lower yield levels. In addition, the pCCA of the effect of year on plant com- All correlations were assessed via the non-parametric Kend- munity composition including plot as a covariate (Fig. A12 all’s tau, which provides a robust measure of association. In Supplementary material) explained 16.4% of total variance the PGE, an increased plant species diversity, expressed as (P < 0.001). Years tended to cluster in discrete periods char- Shannon’s diversity index, number of species, and proportion acterized by the dominance of different species supporting of legumes, was found to be associated with lower, but more the idea of environmental perturbation promoting species stable, biomass yields and vice versa (Table 2 and Table A10 coexistence. For example, Crepis capillaris appears to have Supplementary material), which is in line with findings of been favoured by the environmental conditions in 1991 and recent grassland studies (Haughey et al. 2018; Sanderson 1992 and Leontodon autumnalis in 1999 and 2000. The first 2010). For mean yield and temporal yield variance, the cor- axis discriminated between years dominated by grasses and relations (Kendall’s tau) with plant species diversity indices those dominated by forbs. These results imply that the pro- across treatments were significantly negative ( P < 0.05): for ductivity of diverse plots with a more balanced community yield variance and Shannon’s diversity index (τ = −0.30), in terms of the ratio of grasses to forbs will be more resilient species number (τ  =  −0.30), and proportion of legumes to variability in the environment. The pCCA of the effect of (τ = −0.20); for mean yield and Shannon’s diversity index ‘fertilizer × liming treatment’ (sub-plot) on plant community (τ = −0.15), species number (τ = −0.50), and proportion of composition including year as a covariate explained 64.2% legumes (τ = −0.18) (Table A10 Supplementary material). 1 3 Long-term trends in yield variance of temperate managed grassland Page 15 of 19 37 of total variance (P < 0.001) and discriminated between from aluminum toxicity effects on root growth (Kohler et al. unlimed plots with higher N-input (e.g. 9/2d: ‘N2 P K Na 2016), which may have magnified temporal yield variance. Mg—no chalk’; 11/1d: ‘N3 P K Na Mg—no chalk’) that Hence, liming is a very important management factor for are dominated by grasses and the unfertilized plots (3b/a: stabilizing yields and maintaining legumes in soils prone to ‘Nil—pH 6/7’) that had a higher proportion of forbs (Fig. acidification (Fornara et al. 2011; Storkey et al. 2015), which A13 Supplementary material). in this study has also been proven relevant for supporting sta- Although our results are correlative and plant biodiversity ble grassland productivity over time (Table 1 and Fig. 4) and data were only available for a sub-set of years, they are sup- under varying climatic conditions, particularly limiting the ported by evidence showing that greater plant species rich- water balance and higher temperatures from July to August ness and phylogenetic diversity in managed grasslands may (Table 3 and Fig. 5). enhance their resilience to climate change via enhanced asyn- chrony in the performance of co-occurring species and result in more stable biomass production in response to disturbance 4 Conclusion (Isbell et al. 2017; Eisenhauer et al. 2019; Hautier et al. 2020). This assumes that a large species pool is likely, by chance, Overall, our analysis led to the conclusion that liming, fol- to possess one or two stress-tolerant species that are able to lowed by moderate nutrient supply, promoted plant species resist abiotic stress (e.g. drought) (Kahmen et al. 2005). Such diversity, yield stability, and environmental adaptability stress-tolerant species could compensate for less tolerant spe- and enhanced the long-term sustainability of grassland pro- cies and thus help stabilize productivity in grasslands, and this duction (in terms of stable productivity, biodiversity, and could become more important with respect to climate change reduced synthetic fertilizer inputs). The three research ques- (Loreau and de Mazancourt 2013; Trnka et al. 2021; Haughey tions can be answered as the following: et al. 2018). In the PGE, treatments with high proportions of legumes (‘P K Na Mg’ and ‘N1 P K Na Mg’) had mean yields I. Yes, there were temporal trends in interannual yield close to those in plots with higher rates of inorganic N fertilizer variance over the study period and these varied in while also maintaining higher species richness. However, the relation to fertilizer and lime applications. In particu- yield variance in those plots was similar to that in plots with lar, there were consistent unimodal trends observed moderate fertilizer inputs and sufficient liming (treatment ‘N2 in yield variance in plots with low to moderate or no P K Na Mg – pH 7’; Table 2). These results could be explained nitrogen fertilizer additions, with a peak in the 1990s, by the diversity of fast vs slow functional traits (Storkey and after which variability declined. Yield was most vari- Macdonald 2022; Reich and Cornelissen 2014). Grassland able in plots with higher nutrient inputs and lower communities dominated by slow species were found to stabi- soil pH. lize biomass productivity by increasing mean yield relative to II. The accumulated days of water stress from March temporal yield variance (Craven et al. 2018). to October and mean air temperature from July to Key message for ‘Research question III–part B’: Yield August were identified as the most important climatic variance increased, and plant species diversity decreased with drivers, explaining approximately one-third of the greater fertilizer inputs and reduced liming. Treatments with observed interannual yield variance. Yes, yield sen- high proportions of legumes had mean yields close to those sitivity to these climatic drivers depended on agro- in plots with higher rates of inorganic N fertilizer while also nomic management. Liming and moderate N supply maintaining higher species richness. reduced yield sensitivity to abiotic stresses. The relationship between stable productivity (low yield III. Yes, there was a correlation between plant species variance) and plant species diversity strongly varies depend- diversity, mean yield, and yield variance. Higher ing on agronomic management and soil conditions (Bullock plant species diversity was correlated not only with et al. 2001; Hector et al. 1999; Tracy and Sanderson 2004; more stable grassland yields but also with lower yield Crawley et al. 2005; Storkey et al. 2015). For example, spe- levels. Yield variance increased, and plant species cies diversity decreased sharply with reduced soil pH, par- diversity decreased with greater fertilizer inputs and ticularly under enhanced N supply in the form of ammonium reduced liming. sulfate (Crawley et al. 2005). Acidification due to long-term application of ammonium sulfate on the PGE has selected As a limitation of this study, it should be noted that there for a few grass species tolerant of low pH and has signifi- is a lack of replication in the PGE and that plot errors and cantly reduced the proportion of legumes on plot 11/1 (see ‘fertilizer × liming’ interactions (=residual) could not be Table 2, ‘N3 PK Na Mg—no chalk’, pH 3.5). This acidic separated. Furthermore, biodiversity data were only avail- plot was probably subject to physiological stresses imposed able for a sub-set of years and plots. For this reason, our by low pH (Dodd et al. 1994) and drought effects resulting findings should be interpreted carefully and validated by 1 3 37 Page 16 of 19 J. Macholdt et al. Data Availability The yield dataset of the PGE and the mean annual air further detailed analyses including prospective yield data, temperatures from 1965 to 2018 at Rothamsted were provided by the biodiversity surveys, and soil–climate measurements. electronic Rothamsted Archive and are available from the e-RA web- As possible features of future studies, meta-analyses site (see https://doi. or g/10. 23637/ r pg5-yield s1965- 2018- 01 and https:// of various grassland LTEs under different climate and site doi.or g/10. 23637/ r ms-RMAA tem p-02 , respectively). The precipitation chemistry data from 1992 to 2015 were provided by the UK Environ- conditions may provide further valuable information about mental Change Network (ECN) and are publicly available (see https:// temporal trends in yield variance depending on agronomic doi. or g/ 10. 5285/ 18b7c 387- 037d- 4949- 98bc- e8db5 ef426 4c). The management. In addition to retrospective analyses, grassland annual emissions of sulfur dioxide were provided by the UK National LTEs are also a valuable source for application of agroecosys- Atmospheric Emissions Inventory and are publicly available (https:// naei.beis. go v.uk/ dat a/dat a-selec t or?vie w=air pollut ants ). The species of tem models to simulate grassland responses under contrasting grasses, forbs, and legumes comprising at least 5% of the aboveground soil conditions and under future climate scenarios (Qi et al. biomass found in the surveys of the PGE are available from the e-RA 2018), which should be addressed more in upcoming studies. website (see https://doi. or g/10. 23637/ r pg5-speci es_ 1991- 2000- 01 and Overall, the analysis of long-term grassland experiments, http:// www. era. rotha msted. ac. uk/ datas et/ rpg5/ 01- OAPGs pecies). like this study based on the PGE, could help to improve the Code availability The R-scripts for the statistical analyses generated climate resilience and sustainability of grassland systems by during the current study are available in Tables A6–A7 Supplementary identifying climatic drivers and optimizing the agronomic material or from the corresponding author on reasonable request. management accordingly. In particular, the application of the new designed criss-cross regression approach, in which the Declarations environmental mean was modeled using the selected environ- Ethics approval and consent to participate Not applicable. mental abiotic covariates, allows the assessment of the yield sensitivity (or responsiveness) to changes in climatic condi- Consent for publication Not applicable. tions. The application of this criss-cross regression approach in other agro-ecological trials could help to identify climatic Conflict of interest The authors declare no competing interests. drivers of production risk and to derive agronomic manage- Open Access This article is licensed under a Creative Commons Attri- ment strategies for improving the climate resilience of crop- bution 4.0 International License, which permits use, sharing, adapta- ping systems. This will become increasingly important for tion, distribution and reproduction in any medium or format, as long stable agricultural production in the face of climate change as you give appropriate credit to the original author(s) and the source, and the associated growing risk for abiotic stresses. provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are Supplementary Information The online version contains supplementary included in the article's Creative Commons licence, unless indicated material available at https:// doi. org/ 10. 1007/ s13593- 023- 00885-w. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not Acknowledgements We thank the Rothamsted farm staff, who man- permitted by statutory regulation or exceeds the permitted use, you will aged the PGE and the large teams of people who were involved in the need to obtain permission directly from the copyright holder. To view a fieldwork, vegetation sampling, and sorting between 1965 and 2018. copy of this licence, visit http://cr eativ ecommons. or g/licen ses/ b y/4.0/ . We thank the UK Biotechnology and Biological Sciences Research Council (BBS/E/C/000J0300) and the Lawes Agricultural Trust for supporting the PGE. We thank the Lawes Agricultural Trust and Rothamsted Research for data from the e-RA database. We thank the References UK Environmental Change Network (ECN) for the Precipitation Chem- istry Data. Further thanks are extended to the anonymous reviewers and Addy JWG, Ellis RH, Macdonald AJ, Semenov MA, Mead A (2021) editors for their feedback on earlier versions of this paper. Changes in agricultural climate in South-Eastern England from 1892 to 2016 and differences in cereal and permanent grassland Authors' contributions Janna Macholdt initiated and conceived the yield. Agr Forest Meteorol 308–309:108560. https:// doi. org/ 10. study. All the authors contributed to the further research process. Sarah 1016/j. agrfo rmet. 2021. 108560 Perryman and Tony Scott compiled the datasets. Janna Macholdt, Stef- Addy JWG, Ellis RH, MacLaren C, Macdonald AJ, Semenov MA, fen Hadasch, Hans-Peter Piepho, and Merete Styczen designed and Mead A (2022) A heteroskedastic model of Park Grass spring carried out the statistical analysis. The first draft of the manuscript was hay yields in response to weather suggests continuing yield written by Janna Macholdt, and all the authors commented on previous decline with climate change in future decades. J R Soc Interface versions of the manuscript. All the authors have read and approved the 19(193):20220361. https:// doi. org/ 10. 1098/ rsif. 2022. 0361 final manuscript. Anon (2006) Guide to the classical and other long-term experiments, datasets and sample archive. Lawes Agricultural Trust Ltd, Funding Open Access funding enabled and organized by Projekt DEAL. Harpenden UK. https:// doi. org/ 10. 23637/ rotha msted- long- term- The PGE is supported by the UK Biotechnology and Biological Sciences exper iments- guide- 2006 Research Council (BBS/E/C/000J0300) and the Lawes Agricultural Arata L, Fabrizi E, Sckokai P (2020) A worldwide analysis of trend in Trust. The authors Janna Macholdt and Hans-Peter Piepho acknowledge crop yields and yield variability: evidence from FAO data. Econ support by DFG (Deutsche Forschungsgemeinschaft) grants MA Model 90:190–208. https://doi. or g/10. 1016/j. econm od. 2020. 05. 006 7094/1-1 and PI 377/20-2, respectively. The author Jonathan Storkey Avery BW, Catt JA (1995) The soil at Rothamsted. 44 e-RA. Lawes acknowledges the support by the Biological Sciences Research Council Agricultural Trust Co. Ltd, Harpenden UK. https:// doi. org/ 10. (BBSRC) funded Soils to Nutrition project (BBS/E/C/000I0320). 23637/ ERADOC- 1- 143 1 3 Long-term trends in yield variance of temperate managed grassland Page 17 of 19 37 Baca Cabrera JC, Hirl RT, Schaufele R, Macdonald A, Schnyder H O, Freschet GT, Giling DP, Hattenschwiler S, Hillebrand H, (2021) Stomatal conductance limited the CO response of grass- Hines J, Isbell F, Koller-France E, Konig-Ries B, de Kroon H, land in the last century. BMC Biol 19(1):50. https:// doi. org/ 10. Meyer ST, Milcu A, Muller J, Nock CA, Petermann JS, Roscher 1186/ s12915- 021- 00988-4 C, Scherber C, Scherer-Lorenzen M, Schmid B, Schnitzer SA, Bengtsson J, Bullock JM, Egoh B, Everson C, Everson T, O'Connor Schuldt A, Tscharntke T, Turke M, van Dam NM, van der Plas T, O'Farrell PJ, Smith HG, Lindborg R (2019) Grasslands-more F, Vogel A, Wagg C, Wardle DA, Weigelt A, Weisser WW, important for ecosystem services than you might think. Ecosphere Wirth C, Jochum M (2019) A multitrophic perspective on biodi- 10(2). https:// doi. org/ 10. 1002/ ecs2. 2582 versity-ecosystem functioning research. Adv Ecol Res 61:1–54. Berti A, Marta AD, Mazzoncini M, Tei F (2016) An overview on long-term https:// doi. org/ 10. 1016/ bs. aecr. 2019. 06. 001 agro-ecosystem experiments: present situation and future potential. Fornara DA, Steinbeiss S, McNamara NP, Gleixner G, Oakley S, Eur J Agron 77:236–241. https:// doi. org/ 10. 1016/j. eja. 2016. 01. 004 Poulton PR, Macdonald AJ, Bardgett RD (2011) Increases in Bommarco R, Kleijn D, Potts SG (2013) Ecological intensification: soil organic carbon sequestration can reduce the global warm- harnessing ecosystem services for food security. Trends Ecol Evol ing potential of long-term liming to permanent grassland. Glob 28(4):230–238. https:// doi. org/ 10. 1016/j. tree. 2012. 10. 012 Chang Biol 17(5):1925–1934. https:// doi. org/ 10. 1111/j. 1365- Brookshire EN, Weaver T (2015) Long-term decline in grassland 2486. 2010. 02328.x productivity driven by increasing dryness. Nat Commun 6:7148. Graux A-I, Bellocchi G, Lardy R, Soussana J-F (2013) Ensemble https:// doi. org/ 10. 1038/ ncomm s8148 modelling of climate change risks and opportunities for man- Bullock JM, Pywell RF, Burke MJW, Walker KJ (2001) Restoration of aged grasslands in France. Agric for Meteorol 170:114–131. biodiversity enhances agricultural production. Ecol Lett 4(3):185–https:// doi. org/ 10. 1016/j. agrfo rmet. 2012. 06. 010 189. https:// doi. org/ 10. 1046/j. 1461- 0248. 2001. 00215.x Hadasch S, Laidig F, Macholdt J, Bönecke E, Piepho HP (2020) Carswell AM, Gongadze K, Misselbrook TH, Wu L (2019) Impact of Trends in mean performance and stability of winter wheat and transition from permanent pasture to new swards on the nitrogen winter rye yields in a long-term series of variety trials. Field use efficiency, nitrogen and carbon budgets of beef and sheep Crops Res 252. https:// doi. org/ 10. 1016/j. fcr. 2020. 107792 production. Agric Ecosyst Environ 283:106572. https:// doi. org/ Hall DO, Scurlock JMO (1991) Climate change and productivity 10. 1016/j. agee. 2019. 106572 of natural grassland. Ann Bot 67:49–55. https:// www. jstor. org/ Chang J, Ciais P, Viovy N, Soussana JF, Klumpp K, Sultan B (2017) stable/ 42758 390 Future productivity and phenology changes in European grass- Haughey E, Suter M, Hofer D, Hoekstra NJ, McElwain JC, Luscher A, lands for different warming levels: implications for grassland Finn JA (2018) Higher species richness enhances yield stability management and carbon balance. Carbon Balance Manag in intensively managed grasslands with experimental disturbance. 12(1):11. https:// doi. org/ 10. 1186/ s13021- 017- 0079-8 Sci Rep 8(1):15047. https://doi. or g/10. 1038/ s41598- 018- 33262-9 Craven D, Eisenhauer N, Pearse WD, Hautier Y, Isbell F, Roscher Hautier Y, Seabloom EW, Borer ET, Adler PB, Harpole WS, Hillebrand C, Bahn M, Beierkuhnlein C, Bonisch G, Buchmann N, Byun H, Lind EM, MacDougall AS, Stevens CJ, Bakker JD, Buckley C, Catford JA, Cerabolini BEL, Cornelissen JHC, Craine JM, YM, Chu C, Collins SL, Daleo P, Damschen EI, Davies KF, Fay De Luca E, Ebeling A, Griffin JN, Hector A, Hines J, Jentsch A, PA, Firn J, Gruner DS, Jin VL, Klein JA, Knops JM, La Pierre Kattge J, Kreyling J, Lanta V, Lemoine N, Meyer ST, Minden V, KJ, Li W, McCulley RL, Melbourne BA, Moore JL, O’Halloran Onipchenko V, Polley HW, Reich PB, van Ruijven J, Schamp B, LR, Prober SM, Risch AC, Sankaran M, Schuetz M, Hector A Smith MD, Soudzilovskaia NA, Tilman D, Weigelt A, Wilsey (2014) Eutrophication weakens stabilizing effects of diversity in B, Manning P (2018) Multiple facets of biodiversity drive the natural grasslands. Nature 508(7497):521–525. https://doi. or g/10. diversity-stability relationship. Nat Ecol Evol 2(10):1579–1587. 1038/ natur e13014 https:// doi. org/ 10. 1038/ s41559- 018- 0647-7 Hautier Y, Zhang P, Loreau M, Wilcox KR, Seabloom EW, Borer Crawley MJ, Johnston AE, Silvertown J, Dodd M, de Mazancourt ET, Byrnes JEK, Koerner SE, Komatsu KJ, Lefcheck JS, Hector C, Heard MS, Henman DF, Edwards GR (2005) Determinants A, Adler PB, Alberti J, Arnillas CA, Bakker JD, Brudvig LA, of species richness in the Park Grass Experiment. Am Nat Bugalho MN, Cadotte M, Caldeira MC, Carroll O, Crawley M, 165(2):179–192. https:// doi. org/ 10. 1086/ 427270 Collins SL, Daleo P, Dee LE, Eisenhauer N, Eskelinen A, Fay Deguines N, Jono C, Baude M, Henry M, Julliard R, Fontaine C PA, Gilbert B, Hansar A, Isbell F, Knops JMH, MacDougall AS, (2014) Large-scale trade-off between agricultural intensification McCulley RL, Moore JL, Morgan JW, Mori AS, Peri PL, Pos ET, and crop pollination services. Front Ecol Environ 12(4):212– Power SA, Price JN, Reich PB, Risch AC, Roscher C, Sankaran 217. https:// doi. org/ 10. 1890/ 130054 M, Schutz M, Smith M, Stevens C, Tognetti PM, Virtanen R, Digby PGN (2009) Modified joint regression analysis for incom - Wardle GM, Wilfahrt PA, Wang S (2020) General destabilizing plete variety × environment data. The J Agric Sci 93(1):81–86. effects of eutrophication on grassland productivity at multiple https:// doi. org/ 10. 1017/ S0021 85960 00861 59 spatial scales. Nat Commun 11(1):5375. https:// doi. org/ 10. 1038/ Dodd ME, Silvertown J, McConway K, Potts J, Crawley M (1994) s41467- 020- 19252-4 Application of the British national vegetation classification to Hector A, Schmid B, Beierkuhnlein C, Caldeira MC, Diemer M, Dimi- the communities of the park grass experiment through time. trakopoulos PG, Finn JA, Freitas H, Giller PS, Good J, Harris R, Folia Geobot Phytotx 29(3):321–334. https:// doi. org/ 10. 1007/ Hogberg P, Huss-Danell K, Joshi J, Jumpponen A, Korner C, Lead- bf028 82911 ley PW, Loreau M, Minns A, Mulder CP, O’Donovan G, Otway SJ, Dodd M, Silvertown J, McConway K, Potts J, Crawley M (1997) Sta- Pereira JS, Prinz A, Read DJ, Scherer-Lorenzen M, Schulze ED, bility in the plant communities of the Park Grass Experiment: Siamantziouras ASD, Spehn EM, Terry AC, Troumbis AY, Wood- the relationships between species richness, soil pH and biomass ward FI, Yachi S, Lawton JH (1999) Plant diversity and productiv- variability. Philos Trans R Soc Lond B Biol Sci 346(1316):185– ity experiments in European grasslands. Science 286(5442):1123– 193. https:// doi. org/ 10. 1098/ rstb. 1994. 0140 1127. https:// doi. org/ 10. 1126/ scien ce. 286. 5442. 1123 Eberhart SA, Russell WA (1966) Stability Parameters for Comparing Höglind M, Thorsen SM, Semenov MA (2013) Assessing uncertain- Varieties. Crop Sci 6(1):36–40. https:// doi. org/ 10. 2135/ crops ties in impact of climate change on grass production in North- ci1966. 00111 83X00 06000 10011x ern Europe using ensembles of global climate models. Agric for Eisenhauer N, Schielzeth H, Barnes AD, Barry K, Bonn A, Brose Meteorol 170:103–113. https://doi. or g/10. 1016/j. ag rfor met.2012. 02. 010 U, Bruelheide H, Buchmann N, Buscot F, Ebeling A, Ferlian 1 3 37 Page 18 of 19 J. Macholdt et al. Holland JE, Bennett AE, Newton AC, White PJ, McKenzie BM, George Perryman S, Scott T, Hall C (2019) Rothamsted 30-year meteorological TS, Pakeman RJ, Bailey JS, Fornara DA, Hayes RC (2018) Liming means 1981–2010. Electronic Rothamsted Archive, Rothamsted impacts on soils, crops and biodiversity in the UK: a review. Sci Research. https:// doi. org/ 10. 23637/ OARES 30YrM eans8 180 Total Environ 610–611:316–332. https:// doi. org/ 10. 1016/j. scito Perryman S, Crawley M, Ostler R, Storkey J (2021) Dataset: Park Grass tenv. 2017. 08. 020 species, fertilizer and lime treatments 1991–2000. Electronic Isbell F, Adler PR, Eisenhauer N, Fornara D, Kimmel K, Kremen C, Rothamsted Archive, Rothamsted Research. https:// doi. org/ 10. Letourneau DK, Liebman M, Polley HW, Quijas S, Scherer-Lor-23637/ rpg5- speci es_ 1991- 2000- 01 enzen M, Bardgett R (2017) Benefits of increasing plant diversity Piepho HP (1998) Methods for comparing the yield stability of crop- in sustainable agroecosystems. J Ecol 105(4):871–879. https://doi. ping systems—a review. J Agron Crop Sci 180(4):193–213. org/ 10. 1111/ 1365- 2745. 12789https:// doi. org/ 10. 1111/j. 1439- 037X. 1998. tb005 26.x Kahmen A, Perner J, Buchmann N (2005) Diversity-dependent produc- Piseddu F, Bellocchi G, Picon-Cochard C (2021) Mowing and warm- tivity in semi-natural grasslands following climate perturbations. ing effects on grassland species richness and harvested biomass: Funct Ecol 19(4):594–601. https:// doi.or g/ 10. 1111/j.13 65-24 35. meta-analyses. Agron Sustain Dev 41(6). https://d oi.o rg/1 0.10 07/ 2005. 01001.xs13593- 021- 00722-y Kipling RP, Virkajarvi P, Breitsameter L, Curnel Y, De Swaef T, Preissel S, Reckling M, Schläfke N, Zander P (2015) Magnitude and Gustavsson AM, Hennart S, Hoglind M, Jarvenranta K, Minet farm-economic value of grain legume pre-crop benefits in Europe: J, Nendel C, Persson T, Picon-Cochard C, Rolinski S, Sandars a review. Field Crops Res 175:64–79. https:// doi. org/ 10. 1016/j. DL, Scollan ND, Sebek L, Seddaiu G, Topp CFE, Twardy S, Van fcr. 2015. 01. 012 Middelkoop J, Wu L, Bellocchi G (2016) Key challenges and pri- Prieto I, Violle C, Barre P, Durand JL, Ghesquiere M, Litrico I (2015) orities for modelling European grasslands under climate change. Complementary effects of species and genetic diversity on pro - Sci Total Environ 566–567:851–864. https:// doi. org/ 10. 1016/j. ductivity and stability of sown grasslands. Nat Plants 1(4):15033. scito tenv. 2016. 05. 144https:// doi. org/ 10. 1038/ NPLAN TS. 2015. 33 Kohler IH, Macdonald AJ, Schnyder H (2016) Last-century increases Qi A, Holland RA, Taylor G, Richter GM (2018) Grassland futures in in intrinsic water-use efficiency of grassland communities fave Great Britain—productivity assessment and scenarios for land use occurred over a wide range of vegetation composition, nutrient change opportunities. Sci Total Environ 634:1108–1118. https:// inputs, and soil pH. Plant Physiol 170(2):881–890. https:// doi. doi. org/ 10. 1016/j. scito tenv. 2018. 03. 395 org/ 10. 1104/ pp. 15. 01472 Raman A, Ladha JK, Kumar V, Sharma S, Piepho HP (2011) Stability Lawes JB, Gilbert JH (1859) Report of experiments with different analysis of farmer participatory trials for conservation agriculture manures on permanent meadow land: part i: produce of hay per using mixed models. Field Crops Res 121(3):450–459. https://doi. acre. J Roy Agr Soc Engl 19:552–573org/ 10. 1016/j. fcr. 2011. 02. 001 Loreau M, de Mazancourt C (2013) Biodiversity and ecosystem stabil- Ray DK, Gerber JS, MacDonald GK, West PC (2015) Climate varia- ity: a synthesis of underlying mechanisms. Ecol Lett 16:106–115. tion explains a third of global crop yield variability. Nat Commun https:// doi. org/ 10. 1111/ ele. 12073 6:5989. https:// doi. org/ 10. 1038/ ncomm s6989 Macdonald A, Poulton P, Clark I, Scott T, Glendining M, Perryman Ray DK, West PC, Clark M, Gerber JS, Prishchepov AV, Chatterjee S, Storkey J, Bell J, Shield I, McMillan V, Hawkins J (2018) S (2019) Climate change has likely already affected global food Rothamsted long-term experiments—guide to the classical production. Plos One 14(5):e0217148. https:// doi. org/ 10. 1371/ and other long-term experiments, datasets and sample archive. journ al. pone. 02171 48 Rothamsted Research (UK). https:// www. rotha msted. ac. uk/ sites/ Reckling M, Ahrends H, Chen T-W, Eugster W, Hadasch S, Knapp defau lt/ files/ RRes_ LTE% 20Gui debook_ 2018_% 20web% 20AW. S, Laidig F, Linstädter A, Macholdt J, Piepho H-P, Schiffers K, pdf. Accessed 24 Feb 2023 Döring TF (2021) Methods of yield stability analysis in long-term Macholdt J, Hadasch S, Piepho HP, Reckling M, Taghizadeh-Toosi A, field experiments. A review. Agron Sustain Dev 41(2). https://doi. Christensen BT (2021) Yield variability trends of winter wheat org/ 10. 1007/ s13593- 021- 00681-4 and spring barley grown during 1932–2019 in the Askov Long- Reich PB, Cornelissen H (2014) The world-wide ‘fast-slow’ plant term Experiment. Field Crops Res 264:108083. https:// doi. org/ economics spectrum: a traits manifesto. J Ecol 102(2):275–301. 10. 1016/j. fcr. 2021. 108083https:// doi. org/ 10. 1111/ 1365- 2745. 12211 Midolo G, Alkemade R, Schipper AM, Benítez-López A, Perring Robertson GP, Vitousek PM (2009) Nitrogen in agriculture: balancing MP, De Vries W, Xu X (2018) Impacts of nitrogen addition on the cost of an essential resource. Annu Rev Environ Resour 34:97– plant species richness and abundance: a global meta-analysis. 125. https:// doi. org/ 10. 1146/ annur ev. envir on. 032108. 105046 Glob Ecol Biogeogr 28(3):398–413. https:// doi. or g/ 10. 1111/ Römer T (1917) Are high-yielding varieties more yield stable? Com- geb. 12856 mun DLG 31:87–89 Nabugoomu F, Kempton RA, Talbot M (1999) Analysis of series of Sanderson MA (2010) Stability of production and plant species diver- trials where varieties differ in sensitivity to locations. J Agric Biol sity in managed grasslands: a retrospective study. Basic Appl Ecol Environ Stat 4(3). https:// doi. org/ 10. 2307/ 14003 88 11(3):216–224. https:// doi. org/ 10. 1016/j. baae. 2009. 08. 002 Olesen JE, Bindi M (2002) Consequences of climate change for Euro- Schmidhuber J, Tubiello FN (2007) Global food security under climate pean agricultural productivity, land use and policy. Eur J Agron change. Proc Natl Acad Sci U S A 104(50):19703–19708. https:// 16(4):239–262. https:// doi. org/ 10. 1016/ S1161- 0301(02) 00004-7doi. org/ 10. 1073/ pnas. 07019 76104 Onofri A, Seddaiu G, Piepho H-P (2016) Long-Term Experiments Silvertown J, Poulton P, Johnston E, Edwards G, Heard M, Biss PM with cropping systems: case studies on data analysis. Eur J Agron (2006) The Park Grass Experiment 1856–2006: its contribution 77:223–235. https:// doi. org/ 10. 1016/j. eja. 2016. 02. 005 to ecology. J Ecol 94(4):801–814. https://d oi.o rg/1 0.1 111/j.1 365- Payne RW (2015) The design and analysis of long-term rotation experi-2745. 2006. 01145.x ments. Agron J 107(2):772–785. https:// doi. org/ 10. 2134/ agron Storkey J, Macdonald AJ (2022) The role of long-term experiments in j2012. 0411 validating trait-based approaches to achieving multifunctionality Perryman S, Ostler R (2021) Dataset: Park Grass hay yields, ferti- in grasslands. Front Agric Sci Eng 9(2):187–196. https:// doi. org/ lizer and lime treatments 1965–2018. Electronic Rothamsted 10. 15302/j- fase- 20214 38 Storkey J, Macdonald AJ, Poulton PR, Scott T, Kohler IH, Schnyder H, Archive, Rothamsted Research https:// doi. org/ 10. 23637/ r pg5- Goulding KW, Crawley MJ (2015) Grassland biodiversity bounces yield s1965- 2018- 01 1 3 Long-term trends in yield variance of temperate managed grassland Page 19 of 19 37 back from long-term nitrogen addition. Nature 528(7582):401– Verbyla A, Cullis B, Kenward M, Welham S (1999) The analysis of 404. https:// doi. org/ 10. 1038/ natur e16444 designed experiments and longitudinal data by using smoothing Tilman D, Reich PB, Knops J, Wedin D, Mielke T, Lehman C (2001) splines. J R Stat Soc Series C Stat 48(3):269–311. https:// www. Diversity and productivity in a long-term grassland experiment. Sci-jstor. org/ stable/ 26808 26 ence 294(5543):843–845. https:// doi. org/ 10. 1126/ scien ce. 10603 91 Wilcox KR, Shi Z, Gherardi LA, Lemoine NP, Koerner SE, Hoover Tilman D, Reich PB, Knops JMH (2006) Biodiversity and ecosys- DL, Bork E, Byrne KM, Cahill J Jr, Collins SL, Evans S, Gilgen tem stability in a decade-long grassland experiment. Nature AK, Holub P, Jiang L, Knapp AK, LeCain D, Liang J, Garcia- 441(7093):629–632. https:// doi. org/ 10. 1038/ natur e04742 Palacios P, Penuelas J, Pockman WT, Smith MD, Sun S, White Tracy BF, Sanderson MA (2004) Productivity and stability relation- SR, Yahdjian L, Zhu K, Luo Y (2017) Asymmetric responses ships in mowed pasture communities of varying species compo- of primary productivity to precipitation extremes: a synthesis of sition. Crop Sci 44(6):2180–2186. https:// doi. org/ 10. 2135/ crops grassland precipitation manipulation experiments. Glob Chang ci2004. 2180 Biol 23(10):4376–4385. https:// doi. org/ 10. 1111/ gcb. 13706 Trnka M, Olesen JE, Kersebaum KC, SkjelvÅG AO, Eitzinger J, Zhang Y, Loreau M, Lu X, He N, Zhang G, Han X (2016) Nitrogen Seguin B, Peltonen-Sainio P, RÖTter R, Iglesias ANA, Orlandini enrichment weakens ecosystem stability through decreased spe- S, DubrovskÝ M, Hlavinka P, Balek J, Eckersten H, Cloppet E, cies asynchrony and population stability in a temperate grassland. Calanca P, Gobin A, VuČEtiĆ V, Nejedlik P, Kumar S, Lalic B, Glob Chang Biol 22(4):1445–1455. https:// doi. org/ 10. 1111/ gcb. Mestre A, Rossi F, Kozyra J, Alexandrov V, SemerÁDovÁ D, 13140 ŽAlud Z (2011) Agroclimatic conditions in Europe under climate Zhang Y, Feng J, Loreau M, He N, Han X, Jiang L (2019) Nitrogen change. Glob Chang Biol 17(7):2298–2318. https:// doi. org/ 10. addition does not reduce the role of spatial asynchrony in stabilis- 1111/j. 1365- 2486. 2011. 02396.x ing grassland communities. Ecol Lett 22(4):563–571. https://doi. Trnka M, Balek J, Semenov MA, Semeradova D, Belinova M, Hlavinka org/ 10. 1111/ ele. 13212 P, Olesen JE, Eitzinger J, Schaumberger A, Zahradnicek P, Kopecky D, Zalud Z (2021) Future agroclimatic conditions and implications Publisher's note Springer Nature remains neutral with regard to for European grasslands. Biol Plant 64:865–880. https://doi. or g/10. jurisdictional claims in published maps and institutional affiliations. 32615/ bp. 2021. 005 1 3

Journal

Agronomy for Sustainable Development – Springer Journals

Published: Jun 1, 2023

Keywords: Agronomic management; Biomass production; Climate resilience; Fertilizer input; Food security; Liming; Plant species diversity; Soil pH; Temperature; Water stress

There are no references for this article.