Ecosystem stability at the landscape scale is primarily associated with climatic history

1. There is an increasing interest in landscape-scale perspectives of ecosystem functioning to inform policy and conservation decisions. However, we need a better understanding of the stability of ecosystem functioning (e.g. plant productivity) at the landscape scale to inform policy around topics such as global food security. 2. We investigate the role of the ecological and environmental context on landscape-scale stability of plant productivity in agricultural pasture using remotely sensed enhanced vegetation index data. We determine whether four measures of stability (variability, magnitude of extreme anomalies, recovery


| INTRODUC TI ON
Maintaining stable ecosystem functioning under global environmental change is vital to the continued provision of ecosystems services.To predict ecosystem functioning, such as plant productivity, under future environmental conditions, it is key to identify what underpins their stability so that land managers can be prepared and mitigate the impacts of environmental change.The stability of a system (i.e. its capacity to absorb environmental disturbances, maintain ecosystem functioning and recover from these disturbances; Figure 1A) varies in space and time with differences in the biotic composition and environmental conditions of the system (e.g.Craven et al., 2018).Ecosystem functioning stability in response to experimental disturbances in grassland ecosystems has often been studied at the field scale (e.g.Kreyling et al., 2017;van Ruijven & Berendse, 2010).Across a landscape, however, neither biodiversity nor land cover are constant; meta-communities and habitat patches combine and interact to form complex, heterogeneous landscapes (Oehri et al., 2020).Small-scale field experiments, therefore, are not sufficient to fully capture the processes that generate and maintain ecosystem functioning across a landscape (Bullock et al., 2017;Chalcraft, 2013).The investigations into landscape-scale stability of ecosystem functioning, which incorporate multiple habitat or land cover types, are still relatively uncommon (Oehri et al., 2017) despite their potential to provide insights into large-scale responses of natural systems to realistic disturbance events associated with global change (Huang & Xia, 2019;Mazzochini et al., 2019).
Remotely sensed vegetation indices allow the study of landscapescale ecosystem stability under global change in natural systems (Kéfi et al., 2019), avoiding the restrictions of a single disturbance type (e.g.Kreyling et al., 2017;Vogel et al., 2012) or a limited species set (e.g.van Ruijven & Berendse, 2010).Additionally, land management, conservation practices and policy implementation are often implemented at the landscape scale (Chalcraft, 2013), requiring knowledge on stability at the relevant scale, particularly for ecosystem functions that are important for services such as food security.
We focus here on the stability of plant productivity in agriculture grasslands or pastures.This land cover type is vital for food production and, therefore, it is crucial to understand the dynamics of plant productivity in this system to maintain food security under global climate change.The reduction of plant productivity in agriculture pastures can have serious economic impacts at large scales.For example, in 2012-2013, Ireland experienced a fodder crisis where ecosystem functioning dropped substantially due to a combination biodiversity, climate, ecosystem function, grassland, productivity, remote sensing, resilience, stability F I G U R E 1 (A) The four stability measures investigated are (a) variabilitythe temporal variability in functioning, (b) magnitude of extreme anomalies-the initial change in functioning levels during a perturbation, (c) recovery time-the time taken to return to a baseline level of functioning after a perturbation and (d) recovery rate-the rate of this return, that is engineering resilience.The processes contributing to ecosystem functioning stability emerge at different spatial scales (B), where surrounding landscape heterogeneity shown in the left-hand box can influence a central area (centre).The enhanced vegetation index (EVI) signal of ecosystem functioning will come from the highly productive land covers such as pasture (light green) and other agricultural land (orange) while high biodiversity is likely to be harboured in fragments of marginal land cover types such as forest (dark green) and field boundaries of a poor growing season followed by a cold winter (DAFM, 2017;Green et al., 2018) and resulted in the country having to source additional feed for livestock internationally at a high price (DAFM, 2019).
Identifying the factors that will promote ecosystem functioning stability in these systems, therefore, is vital for sustainable fodder production.
Biodiversity, frequently in the form of species richness, has been shown to promote ecosystem stability at the field scale (e.g.Craven et al., 2018).This relationship supports the insurance hypothesis, where an ecosystem with increased species richness will hold more species capable of maintaining functioning in the face of environmental fluctuations (Yachi & Loreau, 1999).The biodiversitystability relationship, however, shows some inconsistency across studies (Cottingham et al., 2001).For example, neither recovery (De Boeck et al., 2018) nor resistance (Kreyling et al., 2017;van Ruijven & Berendse, 2010) is consistently promoted by plant biodiversity.
At the landscape scale, biodiversity may be less important in determining some ecosystem processes than at small spatial scales.
For example, the role of plant-plant interactions on plant productivity is reduced (Chisholm et al., 2013).Additionally, as spatial scale increases, landscape heterogeneity may mean that biodiversity and the stability of landscape ecosystem functioning each represent different fragments of the landscape (Figure 1B).For example, the dynamics of plant productivity in a landscape containing improved agricultural grassland will predominantly reflect the managed grassland fields, while the biodiversity will predominantly reflect the fragments of habitat such as hedgerows and field margins.
The surrounding landscape can contribute to an ecosystem's functioning and stability (Gunderson, 2000).For example, restored grasslands lying within human-modified landscapes support lower levels of ecosystem functioning than those surrounded by forests and natural grasslands (Zirbel et al., 2019), and the richness of land covers across a landscape can promote inter-annual stability of plant productivity (Oehri et al., 2020).The role of landscape heterogeneity on local stability supports the spatial form (landscape moderated) of the insurance hypothesis (Gonzalez et al., 2009;Tscharntke et al., 2012), which assumes heterogeneous landscapes harbour biodiversity and functionally redundant species disperse to new environments that they are adapted to (Loreau et al., 2003).
Simulations demonstrate that increasing species turnover can increase ecosystem stability through spatial insurance effects (Wang & Loreau, 2016); however, empirical studies testing this hypothesis are lacking.
The large-scale impacts of climate on plant productivity stability at landscape scales need to be considered along with the effect of biodiversity (Craven et al., 2018;García-Palacios et al., 2018;Mazzochini et al., 2019).Firstly, climatic conditions during critical times in a plant's development can dictate plant biomass production (La Pierre et al., 2011), and secondly, the climatic history of an area may influence system processes (Gao et al., 2016), impacting the response of productivity to climatic disturbances.Species occurring in areas that regularly experience extreme climatic events can show a climatic signature, which can leave legacies for ecosystem functioning, that is acclimation, under future extreme events (Craine et al., 2013;Walter et al., 2011).Populations in these areas may show adaption as a result of prior selective pressure and, hence, the ecosystem becomes better able to maintain its functioning at a constant level.Alternatively, climatic variability promotes temporal partitioning of species and, therefore, higher numbers of species can coexist, which may contribute to ecosystem stability through the insurance hypothesis (Hautier et al., 2014).Conversely, if variable regions exhibit low stability of ecosystem functioning, then the system may show no community adaptation and continue to respond to climatic fluctuations similarly over time.This response may mediate, or even outweigh, the relationship of ecosystem stability with other factors, such as biodiversity (García-Palacios et al., 2018;Mazzochini et al., 2019).
The predictions of the impacts of environmental disturbances on ecosystem functioning are challenging due to: the rarity and unpredictability of extreme climatic events (Beniston et al., 2007;Smith, 2011); the synergistic or antagonistic effects of several environmental factors or disturbances (Donohue et al., 2016;Green et al., 2018;Kéfi et al., 2019); and the occurrence of long-term 'press' events (Ives & Carpenter, 2007).Investigating the responses of productivity to specific environmental events in isolation, therefore, can obscure true ecosystem dynamics.Here, we use satellite imagery to investigate the predictors of four stability components of plant productivity at the hectad (10 × 10 km) scale, regardless of the underlying environmental disturbance across the island of Ireland.We hypothesize that the stability of plant productivity is positively associated with: (a) species richness of vascular plants, for example, through the insurance hypothesis; (b) landscape heterogeneity through the landscape-moderated insurance hypothesis; and (c) a history of a variable climate and/or one with many extreme climatic events.We provide new information on how the ecological and environmental context of a landscape impact the stability of plant productivity, which can be used to identify landscapes particularly vulnerable to environmental perturbations associated with global change and be incorporated into ecosystem management.

| Data
We used the enhanced vegetation index (EVI), which measures 'greenness' of an area, to estimate ecosystem functioning as it shows a positive relationship with plant productivity (Sims et al., 2006).
Remotely sensed EVI data for the years 2000-2019 covering the island of Ireland were downloaded from the Moderate Resolution Imaging Spectrometer (MODIS; http://modis.gsfc.nasa.gov/) on NASA 's TERRA and AQUA satellites (products MYD13Q1 and MOD13Q1).Combining these products gives a temporal resolution of approximately 8 days.EVI is particularly superior for grasslands, which dominate in Ireland, to alternative vegetation indices that saturate in high biomass conditions (Huete, 1988;Huete et al., 2002).We relabelled negative values of EVI as 0, as this represents the absence of vegetation within that pixel.The MODIS EVI product has a resolution of 250 × 250 m, which we aggregated to the 10 × 10 km (hectad) scale following the methods in the study by White et al. (2020), so that we could investigate stability at the landscape scale, while also minimizing the impact of spurious pixels on our measures of stability.
We calculated the proportion of each hectad that was covered by each of the 15 mid-level CORINE 2012 land cover classes (http:// www.eea.europa.eu),provided that a hectad was more than 50% land.Pasture-dominated hectads (522 hectads) were defined as having more than 50% of their area classed as pasture (CORINE grid code 18) and were the focus for subsequent analyses.Results for analyses carried out on heterogeneous hectads (defined as having no single land cover accounting for more than 50% of their area, 225 hectads) and all hectads together (830 hectads) are presented in Appendix S5.
To test the insurance hypothesis of ecosystem stability, we estimated species richness using records of vascular plant observations for the years 1970-2017 from the National Biodiversity Data Centre, and the Centre for Environmental Data Recording, which hold biological records from the Republic of Ireland and Northern Ireland respectively.Subspecies were classified at the species level, and hybrids and neophytes (determined using the Online Atlas of British and Irish Flora, https://www.brc.ac.uk/plant atlas/) were excluded from the analyses.We calculated effort-corrected species richness within each hectad using the Frescalo program (Hill, 2012).
See Appendix S1 for details.
The landscape-moderated insurance hypothesis was tested using data on landscape heterogeneity.To calculate the heterogeneity of the landscape surrounding each focal hectad, the Shannon diversity of the 15 mid-level CORINE land covers in the Moore neighbourhood (surrounding eight hectads) was calculated using the package 'vegan' (Oksanen et al., 2020) in R v3.5.1 (R Core Team, 2018).
The association between stability and climate history was investigated using climate data from the E-OBS version 17.0 gridded data product (Haylock et al., 2008; https://www.ecad.eu/download/ ensem bles/downl oad.php).We extracted daily mean temperature and daily total precipitation at the 0.1 degree resolution (approximately 10 km) for Ireland covering the years 1950-1999 and applied nearest neighbour resampling to obtain estimates corresponding to the EVI hectads.The 1950-1999 time period was used to avoid overlap with the MODIS data, thus focusing on the impact of climatic history rather than contemporary climate.We calculated a fat tail measure of climatic extremes and the temporal variance of temperature and precipitation in each hectad.The fat tail measure was calculated as (Q 0.975 − Q 0.025 )/(Q 0.875 − Q 0.125 ), where Q x represents the x quantile of the distribution of either temperature or precipitation in a hectad (Schmid & Trede, 2003).This quantity represents the range of extreme climates relative to the central portion of the data.We also recalculated these climate history metrics using the time period 1980-1999 and 1990-2009, which we used for models of EVI stability calculated from 2000 to 2019 and from 2010 to 2019, respectively, to determine how the length of the time period may affect the results concerning the relationship between stability and climate history.The results for these models are reported in Appendices S2 and S3 respectively.

| Stability calculations
We used EVI to calculate four measures of stability (Figure 1A): magnitude of extreme anomalies (the deviation of plant productivity from 'normal' functioning); recovery time (the time-scale for plant productivity to recover from a perturbation); recovery rate (decrease in plant productivity during a perturbation divided by the corresponding recovery time); and temporal variability (variance through time in plant productivity).
We calculated the magnitude of extreme anomalies, recovery time and recovery rate using a scaled EVI anomaly to remove unwanted seasonal variation from the EVI time series.The scaled EVI anomaly, ∆ EVI (i, t), for hectad i at date t was calculated as: where EVI(i, t) is the EVI of hectad i at date t, m is a month of the year and mean u ∊ m [EVI(i, u)] and sd u ∊ m [EVI(i, u)] are the mean and standard deviation of EVI for hectad i over all dates, u, across the entire period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) falling within month m respectively.This anomaly represents the degree to which any observed EVI value deviates from the expected (mean) EVI of that particular hectad within that particular month.
The magnitude of extreme anomalies, recovery time and recovery rate are based on perturbation events in the EVI anomaly time series.For each hectad, we used an iterative process to identify nonoverlapping perturbations that fell below −2 (i.e. 2 standard deviations below the baseline), which we call 'two-sigma events'.
The magnitude of extreme anomalies was calculated as the absolute value of the EVI perturbation magnitude of each two-sigma event, and is therefore unitless.This measure can be interpreted intuitively as hectads with a high magnitude value represent those with large deviations from the EVI baseline.Recovery time from each two-sigma event was calculated using a moving window average of EVI anomalies spanning approximately 48 days (six MODIS time points) to reduce the noise within the time-series data.
Recovery time was then calculated as the minimum time taken for the moving window average to become non-negative.Recovery rate was calculated as the magnitude of a two-sigma event divided by the recovery time from the two-sigma event.We summarized the average magnitude of extreme anomalies, recovery rate and recovery time in a hectad by their respective means across all twosigma events.
All four measures of stability are based on the methodology presented in the study by White et al. (2020), and more details of these calculations can be found in Appendix S4.A decrease in stability corresponds to either a decrease in recovery rate or an increase in variability, magnitude of extreme anomalies or recovery time.All stability metrics were also calculated using anomalies and measures for the period 2010-2019 to investigate whether varying the climatic history and EVI time period impacted the overall conclusions.

| Statistical analyses
We used generalized least squares models to investigate the predictors of stability.Mean magnitude of extreme anomalies, mean recovery rate and the natural logarithms of variability and mean recovery time for each hectad were taken as response variables.To each of these, we fitted the null intercept model as well as a restricted set of seven models (Table 1) selected to test our hypotheses relating to the importance of biodiversity, land cover heterogeneity and climatic history in explaining spatial variation in the stability of plant productivity.All covariates were scaled and centred to aid with convergence and enable comparisons of effect sizes post hoc.The greater the divergence of the effect size from zero, the more important that variable is in predicting a particular stability measure.
A spatial error term was included within each model by fitting coordinates with an exponential function of the Euclidean distances.
Residuals for all models were normally distributed.A multimodel information theoretic approach (Burnham & Anderson, 2002) was taken using the Akaike information criterion (AIC) to compare model performance, where a difference in AIC of less than 2 indicates models perform equally well.Models were fit across hectads dominated by agricultural pasture.All analyses were carried out in R v3.5.1 (R Core Team, 2018) using the gls function in the r package 'nlme' (Pinheiro et al., 2018).
Before interpreting our results for the effects of climate history, we investigated the effect of temporal autocorrelation in our climate variables (i.e.areas with historical high levels of variation are likely to still be highly variable).We used the psem function in the r package 'piecewiseSEM' (Lefcheck, 2016) for piecewise structural equation models comprised of generalized least squares models to confirm that any effects of climate we captured were those of climatic history as opposed to contemporary effects.The results of these models can be found in Figures S11-S14.

| RE SULTS
For each stability measure, all 'equally performing' models (i.e.models with AICs within two of the models with the lowest AIC, Δ AIC < 2) were retained.The analysis of each stability measure gave results that were consistent across all the retained models (Table 2; The null model for variability performed equally well as the models testing each of our hypotheses (Table 2; Δ AIC < 2 between top performing models).This means there is little support for an association with any of our measured covariates in pasture-dominated hectads (Figure 2).None of our hypotheses, therefore, were supported for predictors of temporal variability of plant productivity.
The only covariate with strong support for an association with the magnitude of extreme anomalies was climatic extremes in temperature, which showed a positive effect.This indicates that areas that have experienced a high level of climatic extremes in temperature in the past experience larger perturbations in plant productivity (effect size = 0.0343, confidence intervals = [0.0107,0.0580]; Figure 3C,E).
The covariate with the strongest effect size on recovery time in pasture-dominated hectads was the positive association with climatic extremes in precipitation (effect size on log recovery time = 0.1022, CI = [0.0608,0.1436]; Figure 4E,F).A history of high variability in  Figure 5D,E).However, the point estimates for the effect size of species richness are consistent with an increase in species richness facilitating faster recovery rates (Figure 5E).
We did not find any support for the landscape-moderated insurance hypothesis as heterogeneity in the land cover of the surrounding landscape was not associated with any measure of stability in pasture-dominated hectads (Figures 2-5) and excluded from many of the best performing models.This result was consistent for heterogeneous hectads (Figure S9) and for recovery time and rate when all hectads across the island of Ireland were considered (Figure S10).
There was some support for a role of landscape heterogeneity in decreasing the magnitude of extreme anomalies across all hectads although it was not included as a covariate in all supported models (Table 2; Figure S10).
Models of all hectads showed strong support for a stabilizing effect of species richness through decreasing temporal variability (Figure S10a) and shorter recovery times (Figure S10c).As with pasture-dominated hectads, past climatic extremes in precipitation conditions had the strongest effect size for recovery time and recovery rate across all hectads as well as in hectads classified as heterogeneous.High temporal variability in past precipitation, however, decreased recovery time when all hectads were considered (Figure S10c), suggesting increased stability and lending support to our climatic history hypothesis of a climatic signature in ecosystem functioning.

| D ISCUSS I ON
Our results show that climatic history can impact multiple measures of ecosystem stability at the landscape scale.This result holds separately for pasture-dominated hectads and heterogeneous hectads, and across all hectads irrespective of the land cover.Despite the suggestion from field experiments that multispecies mixtures can climate-proof agricultural grasslands (Hofer et al., 2016;Mason et al., 2017), at the landscape scale, only recovery time and rate showed any support for an effect of biodiversity in pasturedominated hectads, and was weak at that.Across all land cover types, there was stronger support for a consistent effect of biodiversity across stability measures.The role of the ecological versus the environmental context for stability of productivity, therefore, can depend on whether a specific land cover type is being considered or not.
Our lack of support for a landscape-scale effect of biodiversity in pasture-dominated areas on variability and the magnitude of extreme anomalies, and only weak support for an effect on recovery, is inconsistent with previous biodiversity-stability investigations at the field scale (e.g.Isbell et al., 2015).However, as shown in Figure 1, much of the biodiversity in this study is held in habitat patches outside the agricultural pasture, where only a few species (primarily Lolium perenne) contribute to a high level of productivity.
Therefore, insights gathered from small-scale field experiments on stability mechanisms may not always be sufficient to explain stability in 'real-world' systems at larger spatial scales (Chalcraft, 2013), and the biodiversity-stability relationship is likely scale dependent (Chalcraft, 2013;Delsol et al., 2018), due to heterogeneity in diversity and functioning between habitats within a landscape.This is supported by the role of species richness on multiple dimensions of stability when all land covers are considered (i.e.all hectads).
Additionally, species richness decreased recovery times in pasture when this measure was calculated using only EVI data from 2010 to 2019, although there was no support for an effect of species richness for any other measure for this time period.
Using biological records means that we estimate the impact of realistic communities on ecosystem functioning as opposed to directly manipulating community compositions (Oehri et al., 2017).
However, with this type of data, we were only able to use species richness as a predictor of stability, whereas other biodiversityrelated variables may have a role in maintaining ecosystem functioning including compensatory dynamics, mean-variance scaling and species dominance (Grman et al., 2010).To investigate these mechanisms further, information on species abundances is necessary; however, this is a trade-off of broadscale analyses where presence-only data are often a limiting factor.Functional diversity and phylogenetic diversity are two other components of biodiversity that may also promote ecosystem stability (García-Palacios et al., 2018;Mazzochini et al., 2019).These diversities require accurate community-level data that are difficult to obtain across large geographical extents (but see Mazzochini et al., 2019).
Additionally, our study only looks at ecosystem functioning in terms of plant productivity; however, the multifunctional stability of systems is likely to depend on multiple factors, for example, different biodiversity components and management approaches (Allan et al., 2015;Sasaki et al., 2019).The mechanisms underlying multifunctional stability have begun to be investigated (Sasaki et al., 2019), and have provided further evidence for the need to look at multiple facets of biodiversity.
Our results do not support the landscape-moderated insurance hypothesis as we found no effect of the broader landscape on stability in our study, contrasting with the concept of spatial resilience where there is a hierarchical influence of processes on stability across scales based on location, connectivity and context (Cumming, 2011).The lack of an association with our measures of stability, however, may indicate a limit to the scale at which broader landscape effects operate on these stability components.
Additionally, the detection of land cover heterogeneity will vary with spatial scale, and our approach of using CORINE land cover data may limit the detection of small natural features such as hedgerows, which have been suggested to have disproportionate ecological roles (Hunter, 2017).In this case, small natural features may harbour biodiversity which can contribute to the spatial insurance inherent to the landscape-moderated insurance hypothesis both at and above the hectad scale.
Our study is one of the first to look at historical climate variability and its impact on ecosystem stability.We hypothesized that areas with historical climates that underwent frequent perturbations (either a high variability or relatively frequent extreme events) would lead a system to become adapted to this stress and therefore show greater ecosystem stability in the face of climate perturbations, as previous exposure to extreme droughts can maintain plant productivity through regulation of plant physiological processes (Jentsch et al., 2011).Our results do not consistently support this hypothesis.
For example, a history of extreme precipitation events decreases stability in terms of recovery (increasing recovery time and decreasing recovery rate) across all land cover classes (Figure 2; Figure S2c).
However, in heterogeneous hectads, a history of extreme precipitation events decreases the magnitude of extreme anomalies, and, when all hectads were considered, past variability in precipitation decreased recovery time suggesting that in some instances there may be some sort of legacy in terms of stability.Such an effect, however, is not observed for any stability measure in pasture-dominated hectads, where the destabilizing effect of climatic perturbations F I G U R E 2 Standardized partial regression coefficients showing the relative effect sizes with 95% confidence intervals from gls models of log variability from the best performing models (Δ AIC < 2, Table 2).The null model was the best performing model overall for variability.See Table 1 for abbreviations on productivity may outweigh any positive biodiversity effects (Mazzochini et al., 2019), and already stressed systems appear particularly vulnerable to environmental perturbations.
Our result that a long-term history of extreme events, particularly in precipitation, can lead to slower recovery (longer recovery times and slower recovery rates) may be interpreted as a continued response to perturbations, if areas with a large number of climatic perturbations in the past are assumed to continue to experience similar levels in the present day.To test this, climatic history and contemporary climate were included in a path analysis of stability measures in pasture-dominated hectads showing that the two are related (Appendix S4).Historical extreme event pathways were important in explaining recovery time and rate as well as contemporary climate (Figures S3 and S4), showing that productivity is responding to more than just the magnitude of current perturbations.
By accounting for climate over an extended period of time, similar to the study by Mazzochini et al. (2019), we are able to determine the role of long-term climate in ecosystem functioning as opposed to the more common focus of single extreme weather events (e.g.Jentsch et al., 2011;Vogel et al., 2012).The results for climatic history were similar for when a shorter preceding time period was used (1980-1999;Figures S2-S4).When stability was measured over a shorter time period (2010-2019), the associations with climatic history were lost.Further research is necessary to investigate the impact of climate on stability across multiple temporal scales; however, we believe it is important to maximize the time period that stability is calculated over to avoid characterizing ecosystem stability during periods that experience particularly frequent and large perturbations, for example the multiple fodder crises that occurred between 2010 and 2019 in Ireland (DAFM, 2019;Green et al., 2018).
Investigating how stability varies with temporal scale will also allow the identification of a 'shifting baseline' of stability, that is, are large perturbations to ecosystem functioning becoming more frequent?
Remote sensing data have the advantage of allowing us to estimate ecosystem stability across wide geographical and temporal extents.From this, we can identify features that consistently promote the stability of ecosystem productivity in the face of the multiple types of environmental disturbance, such as increased showing the relative effect sizes with 95% confidence intervals from the best performing models (Δ AIC < 2, Table 2) are shown in (E).See Table 1 for abbreviations extreme climatic events, projected to occur in the future (Beniston et al., 2007).Our approach, using remotely sensed data, builds on the measures of stability from experimental manipulations (e.g.Jentsch et al., 2011;Sasaki et al., 2019), which frequently do not truly represent large-scale or forecasted climatic events (Beier et al., 2012).
Although our measure of the magnitude of extreme anomalies could reflect differences in the strength of perturbations, which we could not measure directly, our measure of recovery rate incorporates the magnitude of the anomaly alongside recovery time, therefore accounting for differences in perturbations.
Our results show some differences between stability within a specific land cover (i.e.pasture) compared to when no specific land cover was considered, particularly for species richness.Pasture covers more than 60% of the island of Ireland and, as an ecosystem, is vital in supporting livestock and ultimately the security of human food supply chains focused on meat and dairy production.Identifying the mechanisms underlying large-scale stability is important for ecological management and sustainability.For example, determining factors that promote the stability of agricultural systems to climate change is vital for sustainable production (Altieri et al., 2015).There has been a call for a paradigm shift in the management of agricultural systems, from a command and control approach for optimal production of natural resources, to one focused on the management of resilience in unstable environments (Altieri et al., 2015).With increased fluctuations in temperature and water availability, agricultural management will be a crucial factor in the stability of grassland production (Vogel et al., 2012).Although land managers cannot change the climate to improve the stability of a landscape's ecosystem functioning, our results can be used to better  et al., 2020), but here we show that different dimensions of stability are associated with different climatic factors.Preparing for temporal variation in production, therefore, will not necessarily decrease the magnitude of or provide faster recovery rates from specific perturbation events.Multiple dimensions of stability need to be considered for sustainable production (Bullock et al., 2017).Additionally, although the effects of biodiversity often had confidence intervals that spanned zero within pasture-dominated hectads, there is a consistent trend towards the positive effect of biodiversity on recovery time and rate, corresponding with results frequently observed at the field plot scale (e.g.Kreyling et al., 2017;van Ruijven & Berendse, 2010).We, therefore, still encourage land managers to maximize the biodiversity within their landscape as this will likely affect both productivity (Oehri et al., 2017) and the recovery of plant productivity following a perturbation at this scale.However, as our results indicate that biodiversity is not the only player, to ensure stable ecosystem functioning across scales and, therefore, sustainable food production, an integrated approach is encouraged incorporating ecological, environmental and socio-economic information.2) are shown in (E).See Table 1 for abbreviations

Figure 2 )
Figure2).In the text, we report only the standardized partial regression coefficients of the best supported model (i.e. the model with the lowest AIC) for pasture-dominated hectads; however results for all of the retained models are shown in the figures.The standardized partial regression coefficients provide information on the sensitivity of the response to individual predictors while controlling for the effects of the other predictors and can be compared to determine their relative importance in explaining the response.See Appendix S5 for the equivalent results from heterogeneous and all land covers.
temperature also increased recovery time (effect size log recovery time = 0.0687, CI = [0.0064,0.1310]; Figure4B,F).There was also some support for an association between recovery time and biodiversity, with richer areas showing shorter recovery times (effect size on log recovery time = −0.0415,CI = [−0.0809,−0.0021]; Figure4A,F).The only covariate with strong support for an association with recovery rate was the negative association with climatic extremes in precipitation (effect size = −0.0132,CI = [−0.0195,−0.0070];

F
The fitted relationships between the magnitude of extreme anomalies and (A) historical variance in temperature, (B) historical variance in precipitation, (C) climatic extremes in temperature and (D) climatic extremes in precipitation, showing the partial regression coefficients with 95% confidence intervals and partial residuals from the climate only gls.The standardized partial regression coefficients Perturbations in plant productivity, such as the 2012-2013 fodder crisis (DAFM, 2017;Green et al., 2018), where ecosystem productivity in Ireland dropped substantially in response to a series of climatic conditions, therefore, have serious economic consequences.By looking at stability within pasture-dominated areas, we can determine the factors which may promote resilient agricultural landscapes and, thus, have important economic applications.

F
The fitted relationships between recovery time (on the scale of the natural logarithm of the number of days to return to the zero baseline of plant productivity) and (A) corrected species richness, (B) historical variance in temperature, (C) historical variance in precipitation, (D) climatic extremes in temperature and (E) climatic extremes in precipitation, showing the partial regression coefficients with 95% confidence intervals and partial residuals from the diversity and climate gls.The standardized partial regression coefficients showing the relative effect sizes with 95% confidence intervals from the best performing models (Δ AIC < 2, Table2) are shown in (F).See Table1for abbreviations prepare for future climate change.Given current climatic conditions, we can identify regions likely to have reduced stability, and this can then be fed into risk analyses and inform management practices such as stocking densities and fodder storage (Department of Agriculture, Food and the Marine, 2017).It also highlights the need for land managers to consider recovery from an environmental perturbation as well as the perturbation itself.Many previous investigations of stability have focused on temporal variation (e.g.van 't Veen This publication has emanated from research conducted with the financial support of Science Foundation Ireland and the Department for the Economy, Northern Ireland under Grant number (15/ IA/2881).The authors thank the many citizen scientists who collected and contributed their data to the National Biodiversity Data Centre and the Centre for Environmental Data and Recording.Open access funding enabled and organized by IRel.F I G U R E 5 The fitted relationships between recovery rate (on the scale of one standard deviation in enhanced vegetation index per day) and (A) historical variance in temperature, (B) historical variance in precipitation, (C) climatic extremes in temperature and (D) climatic extremes in precipitation, showing the partial regression coefficients with 95% confidence intervals and partial residuals from the climate only gls.The standardized partial regression coefficients showing the relative effect sizes with 95% confidence intervals from the best performing models (Δ AIC < 2, Table

652.049 0 Divclim −1,651.828 0.221 Full −1,650.612 1.437 Landclim −1,650.327 1.722
TA B L E 2 AIC values for all models of the four stability measures in pasture-dominated hectads across the island of Ireland.The top performing models with a Δ AIC ≤ 2 are highlighted in bold, indicating that these models perform 'equally well'.See Table1for model abbreviations