Does agricultural intensification cause tipping points in ecosystem services?

Agricultural intensification is being widely pursued as a policy option to improve food security and human development. Yet, there is a need to understand the impact of agricultural intensification on the provision of multiple ecosystem services, and to evaluate the possible occurrence of tipping points. To quantify and assess the long-term spatial dynamics of ecosystem service (ES) provision in a landscape undergoing agricultural intensification at four time points 1930, 1950, 1980 and 2015. Determine if thresholds or tipping points in ES provision may have occurred and if there are any detectable impacts on economic development and employment. To quantify and assess the long-term spatial dynamics of ecosystem service (ES) provision in a landscape undergoing agricultural intensification at four time points 1930, 1950, 1980 and 2015. Determine if thresholds or tipping points in ES provision may have occurred and if there are any detectable impacts on economic development and employment. We used the InVEST suite of software models together with a time series of historical land cover maps and an Input–Output model to evaluate these dynamics over an 85-year period in the county of Dorset, southern England. Results indicated that trends in ES were often non-linear, highlighting the potential for abrupt changes in ES provision to occur in response to slight changes in underlying drivers. Despite the fluctuations in provision of different ES, overall economic activity increased almost linearly during the study interval, in line with the increase in agricultural productivity. Such non-linear thresholds in ES will need to be avoided in the future by approaches aiming to deliver sustainable agricultural intensification. A number of positive feedback mechanisms are identified that suggest these thresholds could be considered as tipping points. However, further research into these feedbacks is required to fully determine the occurrence of tipping points in agricultural systems.

I. To quantify and assess the long-term spatial dynamics of ecosystem service (ES) provision in a landscape undergoing agricultural intensification at four time points 1930, 1950, 1980 and 2015. II. Determine if thresholds or tipping points in ES provision may have occurred and if there are any detectable impacts on economic development and employment.
Methods We used the InVEST suite of software models together with a time series of historical land cover maps and an Input-Output model to evaluate these dynamics over an 85-year period in the county of Dorset, southern England.
Results Results indicated that trends in ES were often non-linear, highlighting the potential for abrupt changes in ES provision to occur in response to slight changes in underlying drivers. Despite the fluctuations in provision of different ES, overall economic activity increased almost linearly during the study interval, in line with the increase in agricultural productivity.
Conclusions Such non-linear thresholds in ES will need to be avoided in the future by approaches aiming to deliver sustainable agricultural intensification. A number of positive feedback mechanisms are identified that suggest these thresholds could be considered as tipping points. However, further research into these feedbacks is required to fully determine the occurrence of tipping points in agricultural systems.

Introduction
Agriculture is one of the most significant drivers of global land use change (Tilman et al. 2011), with around 40% of all ice-free land on earth now being used for food production (Ellis et al. 2010). The intensification of agriculture is one of the most widespread policy options being pursued to improve food security and human development, and typically involves increased use of mechanisation, fertilisers and pesticides to increase agricultural productivity (Dawson et al. 2019). Although such approaches have successfully been used to boost food production in many parts of the world, the negative environmental impacts of intensive farming are well established, and include widespread biodiversity loss; emissions of greenhouse gases (GHG) and ammonia; soil compaction, depletion and erosion; eutrophication; dispersal of toxic herbicides and pesticides; and depletion of freshwater (McLaughlin and Mineau 1995;Green et al. 2005;Balmford et al. 2018). It is likely that these impacts will increase globally in coming decades owing to projected increases in food demand owing to human population growth and dietary shifts towards increasing meat consumption (Clark and Tilman 2017). Agricultural intensification can therefore be considered inconsistent with the sustainable use of land, and will likely undermine achievement of the UN's sustainable development goals (SDGs) (Dudley and Alexander 2017;Dawson et al. 2019). Understanding how the increasing global demand for food can be reconciled with the need for sustainable land use practices represents a significant scientific challenge. Specifically, there is a need to identify how agricultural systems can be developed, in which production is maintained or increased while environmental outcomes are enhanced (Tscharntke et al. 2012;Pretty 2018). In other words, an understanding is required of how 'win-win' outcomes can be achieved from farming where benefits for human well-being are obtained simultaneously with positive outcomes for ecosystems (Rasmussen et al. 2018). The concept of ecosystem services (ES), or the benefits provided by ecosystems to people, can be of particular value in this context. 'Win-Win' outcomes would involve enhanced provision of a range of different ES, such as water and soil conservation, soil carbon storage, nutrient recycling and pest control, as well as food production (Pretty 2018). However, assessments of ES provision in agricultural landscapes have identified widespread trade-offs between food production and other ES, indicating that synergistic outcomes are often difficult to achieve in practice (Power 2010;Newton et al. 2012;Howe et al. 2014;Holt et al. 2016;Bernués et al. 2019).
Development of environmentally sustainable approaches to agriculture therefore requires an understanding of the relationships between different ES and the reasons why trade-offs occur (Howe et al. 2014). Although the scientific understanding of ES has developed rapidly in recent years, much of the research that has been undertaken assumes simple linear relationships and effects in the delivery of services of over different spatial and temporal scales. Yet it is well established that the connections between ecosystem processes, functions and benefits to humans are complex and dynamic, suggesting that relationships between these variables may often be non-linear rather than linear (Costanza et al. 2017). There is therefore a need for research into the spatial dynamics of ES in agricultural landscapes, to identify the factors responsible for changes over time in ES provision. Specifically, there is a need to determine whether nonlinear relationships occur between ES provision and underlying drivers. This reflects growing concern that environmental change might lead to abrupt declines in the condition of ecosystems and associated provision of ES, with consequent negative impacts on human well-being (Mace et al. 2015). The potential for such abrupt declines is based on a growing body of theory relating to critical transitions between alternative ecosystem states (Scheffer et al. 2001(Scheffer et al. , 2012 and empirical observations of ecological thresholds (Huggett 2005); Groffman et al. 2006). In this context, the concept of tipping points has received particular research attention (Moore 2018), for example in relation to the global climate system (Lenton and Ciscar 2013;Lenton et al. 2019) and in marine ecosystems (Selkoe et al., 2015;Watson et al., 2018). Tipping points can be considered as ecological thresholds that are driven by positive feedback loops, in accordance with dynamical systems theory (van Nes et al. 2016). However, abrupt ecosystem changes can also be caused by rapid changes in the underlying drivers, or from an interaction between drivers, without a positive feedback mechanism (Andersen et al. 2009;Newton 2021). There is a need to differentiate between these various possibilities, in order for appropriate management and policy responses to be developed (Watson et al. 2018).
We therefore identify the impact of agricultural intensification on provision of ES and the potential occurrence of tipping points as a key knowledge gap. We are not aware of any prior multidecadal study that has investigated the spatial and temporal dynamics of ES provision in a landscape undergoing agricultural intensification, and has used these data to examine the occurrence of tipping points. Here we address this knowledge gap by examining how ES provisions have changed in the southern English county of Dorset over 85 years, using a unique set of land use maps spanning that period. Today, around 71% of Dorset's land area is farmed, which is a similar value to England as a whole (DEFRA 2019). Since the 1930s, in common with much north-western Europe, agriculture in this county has undergone intensification of land use practices, leading to current conflicts in delivery of multiple ES and biodiversity (Hooftman and Bullock 2012;Jiang et al. 2013). This process has involved an increase in inputs, including greater use of machinery and an increase in application of agrochemicals, together with changing patterns of husbandry of both plants and animals, leading to a concomitant increase in agricultural productivity (Tscharntke et al. 2005(Tscharntke et al. , 2012. These trends are largely attributable to changes in agricultural policy and practice that have occurred in the UK, including those attributable to an increased emphasis on food self-sufficiency in the wake of both World War I and II, and the provision of incentives through the EU's Common Agricultural Policy after accession of the UK in 1973 (Brassley 2000;Robinson and Sutherland 2002;Ollerton et al. 2014). Dorset therefore provides an exceptional opportunity to evaluate long-term dynamics in ES owing to the wealth of environmental datasets that exist, including an adapted land cover map (Hooftman and Bullock, 2012) developed from Stamp (1931) and an extensive botanical survey (Good 1937) dating from the 1930s.
Using these resources, we employ a time series of land cover maps of Dorset for the periods 1930, 1950, 1980 and 2015 based on the availability of detailed habitat data from Ridding et al. (2020aRidding et al. ( , 2020b. We then develop models and proxies of ES provision to evaluate changes in ES delivery at multiple time-steps over an interval of 85 years, during which agricultural intensification took place. Using this approach, we examine whether ES provision has demonstrated nonlinear responses to increasing agricultural intensification, and evaluate whether such responses might represent tipping points. Specifically, this research was designed to test the following hypotheses: (i) when landscapes are subjected to increasing agricultural intensification, the provision of multiple ES will decline as food production increases; (ii) the declines in ES may demonstrate threshold responses and can be driven by positive feedback mechanisms; (iii) thresholds or tipping points in ES provision may have detectable impacts on economic development and employment.

Materials and methods
To investigate the above hypotheses, we used newly digitized British land-use maps from the 1950s and 1980s (Ridding et al., 2020a(Ridding et al., ,2020b combined with previously digitized maps from 1930 and 2015. These specific time points were chosen owing to the unusual availability of time-series and field data describing both biodiversity and land cover. This allowed us to map, at a fine resolution, the extent and spatial details of land-use patterns for the southern English county of Dorset. We applied an ecological production function approach using the InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs) model and a series of proxy indices based on the land cover map categories to link habitat type (Wilson and Hoehn 2006) to ES delivery. Using this approach, we modelled carbon sequestration and storage, water yield, water purification (nitrogen retention and export), flood protection, soil quality, timber production, food production (crops & livestock), recreation, aesthetic value, habitat quality for pollination, changes in BAP (Biodiversity Action Plan) species and suitable habitat for BAP species. We then produced maps at each time step to provide estimates of ES delivery from different land-use classes, which were summed at the landscape scale at each time step, enabling trends in ES provision to be described. Using these trend data, we then used Hierarchical Divisive Estimation, two dimensional Kolmogorov-Smirnov tests and general additive models (GAMs) to investigate potential thresholds of change and the possibility of non-linear relationships between agricultural drivers and ES. Finally, to investigate whether the observed changes in ES flows have had detectable impacts on economic development and employment, we conducted an Economic Impact Analysis (EIA) of the region that provides an estimate of the total employment impacts over the study period as well as estimates of total Gross Value Added (GVA) made by an individual industry or sector.

Study area and land cover maps
The county of Dorset lies in on the south coast of England. The current landscape comprises a mosaic of different land cover types including cropland, grassland, woodland, urban areas and coastal margins. To analyse the impacts of land use changes on flows of ES, we used a time series of land cover maps for 1930for , 1950for , 1980for . The 1930 map was produced by Hooftman and Bullock (2012) based on that of Stamp (1931), whereas the 2015 map was derived from the 2015 CEH Land Cover Map (Rowland et al. 2017). The intermediate maps (1950 and 1980) were produced using the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) Rule Based Scenario Generator tool (Sharp et al. 2016). This employed transition probabilities between land cover types that were cross-validated against historical Dorset habitat survey data and demonstrated a high level of accuracy (87 and 84%, respectively) and low levels of model uncertainty. Full details of how these maps were produced are provided by Ridding et al. (2020aRidding et al. ( , 2020b. Dorset is currently ca. 2653km 2 in area. Prior to an extension of county boundaries in 1974, where the towns of Bournemouth, Poole and Christchurch were added, the area of the county was ca. 2500 km 2 (Hooftman and Bullock 2012). Unless stated otherwise, this latter area was the focus of the current study, enabling changes over time to be assessed.

Ecosystem service assessment
The InVEST suite of tools has been developed to explore how changes in ecosystems can lead to changes in the flows of many different benefits to people (Sharp et al. 2016). The models are based on production functions covering a wide variety of ES, and have successfully been applied in a large number of different contexts (Bagstad et al. 2013). Here we used InVEST to model carbon sequestration and storage, water yield, nitrogen retention and export, crop production and recreation. The water yield and nutrient retention InVEST models were specifically selected because they have recently been tested and validated using contemporary data in the Dorset region (Redhead et al. 2016(Redhead et al. , 2018. Carbon sequestration data required for the Carbon Storage and Sequestration InVEST model were also readily available for Dorset (Jing et al. 2013). Where an InVEST model was unavailable for a particular ES or considered unsuitable, an extended benefit transfer approach was utilized incorporating indices based on the land cover map categories linked to ES delivery. The ES that were mapped using these proxy values included: flood regulation, timber production, livestock production, soil quality, aesthetic value, habitat suitability for pollinators and biodiversity. The proxy values for each land cover type were derived from Newton et al. (2012) and Hodder et al. (2014) and who previously assessed spatial patterns in ES across Dorset. All ES maps were produced at a resolution of 100 9 100 m, consistent with recommendations of other studies (e.g., Redhead et al. 2018;Ridding et al. 2020b) and processed using ArcGIS v 10.1 (ESRI, Redlands, California, USA). We briefly describe the major features and data inputs for the ES models below; additional detail is provided in Appendix 1.

InVEST models
Carbon sequestration and storage The InVEST (v3.4.4) Carbon Storage and Sequestration: climate Regulation model was used to map carbon storage within each land cover type. Carbon pools data for input into the model were extracted from Jiang et al.
(2013) for above-ground stored carbon in biomass (Mg ha -1 ), below-ground stored carbon in biomass (Mg ha -1 ), carbon stored in soil and carbon stored in dead organic matter in Mg ha -1 (see Appendix 1). The sum of these categories provides an estimate of the biophysical amount of carbon stock per 1-ha grid cell. Market or social costs of carbon were not considered in this study.
Water yield, nutrient retention and export The InVEST (v3.4.4) Annual Water Yield model was used to identify the relative contributions of water from different parts of the landscape of Dorset, offering insight into how changes in land use patterns have affected annual surface water yield. The InVEST model requires five biophysical parameters, in addition to the baseline land cover map. These are root restricting layer depth (mm), plant-available water content (PAWC, as a proportion), average annual precipitation (mm), average annual potential evapotranspiration (PET, mm) and a watershed vector layer. We obtained these data from a variety of sources (see Appendix 1). A seasonality constant and the biophysical table of plant evaporation coefficients (Kc) were also included by matching class descriptions in Redhead et al. (2016), who validated the InVEST Annual Water Yield model for the UK and provided parameter values. The InVEST (v.3.4.4) Nutrient Delivery Ratio (NDR) model was then used to calculate the retention and export of nitrogen (N) in terrestrial vegetation across the catchments of Dorset. This model had the same data needs as the Annual Water Yield model, with the addition of a biophysical table, with water quality coefficients including nitrogen loading, and vegetation filtering values for each pixel (see Appendix 1). Additionally, the NDR model also required the Flow Accumulation rate from InVEST (v.3.5.0) RouteDEM. We used the OS 50 m digital terrain model (DTM) (Ordnance Survey (GB) 2015) for the digital elevation model (DEM).
Crop production and livestock production The InVEST (v3.5.0) crop production regression model was used to provide estimates of historical wheat and barley yields (tonnes/ha) in Dorset. To run this model, estimates for average nitrogen, phosphate and potash application rates for each crop were taken from Nix (2018). It is key to note yield data were converted to modelled crop yields (tonnes/ha) using crop production percentiles from the year 2015. Hence, changes in production values reflect land-use change and/or intensification but not how improvements in crop productivity may have been influenced by more modern innovations, e.g. better monitoring technologies, genetically engineered crop-based agriculture or improved methods of pest control. We consider changes in relative crop productivity in the Discussion. Furthermore, the 2015 land cover map uses specific data about agricultural land use, often including the exact crop planted; but the 1930s, 1950s and 1980s maps do not provide this level of detail. To address this, we weighted the output of the 1930-1980 regressions by the proportion of land given over to each crop type and historical fertiliser application data that were obtained from Agricultural Census survey data (Tavener (1952) and DEFRA 2015), Wheat and barley comprised over 70% of the arable land use in all four time periods. Recorded livestock numbers in Dorset (cattle, pigs, sheep and poultry) were also sourced for the 1930s and 1950s from Tavener (1952) and from UK June Agricultural Census survey data (DEFRA 2015) for the 1980s and 2015, as a proxy for livestock production.
Recreation The InVEST (v3.4.4) Visitation: Recreation and Tourism model was used to map potential recreational value across the study area. We parameterized the model using a proxy for visitation, namely geotagged photographs posted across Dorset in 2015 to the website Flickr (Yahoo 2018), following Wood et al. (2013). As this method could only be applied for the 2015 land cover map, owing to the lack of historical data, these estimates were combined with the Scenario Analysis component of the model to explore how historic changes to the landscape could have affected visitation rate based on current patterns of recreational use, using a least squares regression method (see Appendix 1). It is important to note that the model assumes that people's responses to broad habitats (that serve as predictors in the model) will not change over time. In other words, in the past, people were assumed to be attracted to, or repelled by, the predictors in the same way that they are currently. This of course is an oversimplification and such assumptions are discussed further in the Discussion section. As household gardens and urban parks could not be distinguished from other urban areas they were excluded from further analysis. Improved grassland and arable land were also assumed to have little or no recreational value, following Newton et al. (2012).

Indicator and proxy-based models
Flood regulation The InVEST Urban Flood Risk Mitigation model (v3.4.4) was considered for this ES but was not implemented owing to its narrow focus on built infrastructure, which has limited relevance to this study area. Instead, the capacity of vegetation to mitigate flood risk was assessed utilising a proxybased scoring approach developed by Hodder et al. (2014). Differences in land cover will affect flood risk through effects on surface roughness or infiltration capacity, which will affect water retention rates, and hence the volume and timing of flow (Nelson et al. 2009). Land cover classes were each given a score (1-10) based on modelling evidence and expert judgement (see Appendix 1).
Timber production The sustainable maximum yield of timber that could be produced from the management of woodlands in Dorset was estimated using a set of coefficients produced by the Forestry Commission (2008). This value represents sustainable wood yield (tonnes) that can be locally harvested and is available for conventional markets. The average productivity of timber on public and private woodlands was estimated by multiplying the total area of coniferous and broadleaved woodland in private and public ownership by the approximate biomass yield derived from the mean yield class for trees in England (see Appendix 1).
Soil quality Estimated soil erosion rates in terms of degree (intensity) and extent (spatial distribution) of soil degradation were used to reclassify land cover categories using a proxy method based on Graves et al. (2015) who previously estimated soil erosion rates in England and Wales by land use/soil type category. Erosion rates for four soilscapes (clay, silt, sand and peat) were averaged and used to create a single ''soil quality'' index for each land cover type (see Appendix 1).
Aesthetic value ('naturalness') Aesthetic Value was assessed using Campaign to Protect Rural England (CPRE) Tranquility Mapping data (Jackson et al. 2008). The CPRE 'naturalness' indicator is a proxy based on the sense of tranquillity people feel at a given location (see Appendix 1). Here CPRE 'naturalness' scores were obtained for each land cover type, from existing data (Jackson et al. 2008).
Habitat quality for pollinators Spatial values for nectar productivity, species nectar diversity and functional nectar diversity for different land cover types were extracted from Baude et al. (2016). All three measures were normalised and weighted equally to create an index of habitat quality for pollinators (see Appendix 1). To include urban areas within this study, pollinator abundance was scaled from farmland values using results presented by Baldock et al. (2015).
Biodiversity Two approaches were used to examine the potential impacts of land use change on biodiversity, involving (i) a species richness indicator developed by Newton et al. (2012) and (ii) a measure of UK Biodiversity Action Plan (BAP) 'priority' habitat area available to species of conservation concern (UK BAP species). Values of both BAP species number and density were based on compiled field observations that were subsequently normalised to provide measures of the relative biodiversity value of each land cover type (see Appendix 1). Secondly, the total habitat area available to BAP species was recorded by totalling the area of suitable land cover types that each species was associated with (see Appendix 1).

Identification of thresholds in ecosystem services
To examine changes in the flow of ES over time, the ES pixel values for each year (1930, 1950, 1980 and 2015) were calculated for each land cover type and then summed across all land cover types to give an overall value for each ES at each time step, at the county scale (Appendix 2). Hierarchical Divisive Estimation was then applied to each of the ES over the 85-year interval (1930-2015) to detect non-linear trends. This nonparametric statistical technique (James and Matteson 2013) allows for the estimation of multiple change points (or breakpoints) partitioning the time series into contiguous segments, based on an iterative procedure to locate each single break-point, with no distributional assumptions other than the existence of certain absolute datapoints. The statistical significance of an estimated breakpoint is determined through a permutation test. Normal distribution of data was verified by the Kolmogorov-Smirnov test (Massey 1951). As breakpoint analysis often considers an interpolation function between successive points (e.g. Samhouri et al. 2017), we used a stepwise linear interpolation function to normalise each ES time series prior to the analysis to create a single annual time series for each variable. This method was implemented in R (R Core Team 2016) using the ecp (nonparametric multiple change-point analysis of multivariate data) package (James and Matteson 2013). The method was applied to all possible cut off lengths for each time series to test the sensitivity of results obtained from ecp analysis. A cut of length of one for each of the regimes to be tested was chosen as it was the highest sensitivity available and thresholds are often associated with short periods of variability (Watson et al. 2018). Quantitative estimates of nonlinear thresholds were defined as the point of inflection where the second derivative changes sign (Large et al. 2015).

Relationship between ES trends and agricultural drivers
To determine appropriate pressure-state relationships between agricultural intensification and ES flows, time series data for five agricultural metrics were sourced from the literature and environmental databases (see Appendix 3). These metrics were used as proxies for the process of agricultural intensification in Dorset. Data for the total amount of nitrogen-based fertilizer applied in Britain were sourced from Robinson and Sutherland (2002). Total income from farming (TIFF) in the UK and total factor productivity of UK agriculture (TFPA) were sourced from national statistics (DEFRA 2012). Data describing mechanisation (tractor numbers) were taken from FAOSTAT (2018), while data for the area sprayed with insecticide were obtained from Potts et al. (2010). We then fitted generalised additive models (GAMs) to test for potential linear and non-linear relationships between the time-series variables. We tested for non-linearities in all possible pressure-state relationships, but subsequently excluded those without plausible mechanistic relationships. While this approach increased the possibility of detecting non-linearities and tipping points, our interest was on the precautionary identification of thresholds rather than statistical significance, especially given the limitations of inferences based on p-values (e.g., Samhouri et al. 2017). An eigenvalue optimisation process was carried out to prevent overfitting using the ''mgcv'' package (Wood 2011). Generalised cross validation (GCV) was used to estimate a smoothing parameter for each term. Through this eigenvalue optimisation process, smoothing terms with linear functions in response to pressure variables were removed from the model if they did not improve the fit (Wood and Augustin 2002).

Analysis of economic trends
To examine the relationship between ES dynamics and economic trends, we used two approaches to assess economic trends for Dorset. Firstly, we used a combination of office for national statistics (ONS) data and the Cambridge econometrics local economic forecasting model (LEFM; Cambridge Econometrics 2015) to create baseline projections for the Dorset area between 1970 and 2015 using traditional economic metrics of Gross Value Added (GVA) and full-time equivalent (FTE) employment. Second, we used a Standard Industrial Classification (SIC) sector approach (Hughes 2008) to examine Dorset's environmental services and goods sectors in terms of a flow of value. For this part of the analysis we included the towns of Bournemouth and Poole owing to their disproportionately large contribution to the economy of Dorset. We identified 18 industry sectors as having relatively strong links to the environment (Appendix 4). We then conducted an Economic Impact Analysis (EIA) of Dorset's economy for the period 1981-2015 using an 'Extended Input-Output Model' (Appendix 4; Newton et al. 2021). An industrial breakdown back to 1970 was not possible because sector definitions have changed markedly over this period. The model is instead based on a set of more recent  economic and social accounts (see Appendix 4) that allowed analysis of the structure of SIC relationships within the economy. The model was used to calculate historic trends in GVA and FTE employment.

Trends in land cover
In the 1930s, Dorset was dominated by semi-natural (neutral unimproved) grassland pastures (* 41% of area) with considerable swathes of calcareous grassland (19%) and smaller patches of acid grassland (1.7%). Heathland (5.5%) and broadleaved woodland (7.7%) were also significant land cover types (Fig. 1). By the 1950s there was a considerable shift towards arable cropland and improved grassland, which increased by approximately 28 and 6.3% respectively. Neutral grassland and calcareous grassland habitats both declined by [ 50% of their values in the 1930s as large areas were converted to intensive agriculture. By 1980 neutral grassland and calcareous grassland had almost been eliminated in the county (\ 2% of total habitat area), accompanied by marked declines in heathland and acid grassland, which together represented \ 3% of total habitat area. During this time step, there was also a sharp transition from broadleaved woodland to coniferous woodland, which increased from 0.1 to 3% of habitat area. Areas of different land cover types in 2015 were similar to those of 1980, but in this period, there were continuing declines in arable cropland, broadleaved woodland and calcareous grassland at the expense of improved grasslands and coniferous woodlands. All other habitat types remained fairly constant over the study period.

Trends in ecosystem services
Breakpoint analysis identified significant abrupt thresholds in all but two of the ES over the 85 years studied ( Table 1). Numbers of livestock increased significantly between 1930 and 1960 (breakpoint P = 0.025, Table 1) and continued to increase at a higher rate between 1960 and 1980 (breakpoint P = 0.047, Table 1), then gradually declined thereafter to 2015 (Fig. 2). Arable crop production increased substantially (? 178 M tonnes) between 1930 and 1950, but then declined thereafter, the rate of decline reducing after 1980. The 1930-1950 increase is most noticeable in those areas that were already arable (e.g. the central chalklands) and in the westerly parts of the heathlands, where a considerable amount of land was reclaimed or recultivated for cereals (Appendix 5). However, no breakpoint was detected between 1930 and 2015. Production of broadleaved timber declined steadily from 1930 to 1980, with a significant breakpoint detected in 1960 (P = 0.011, Table 1). The downward trend in broadleaved timber harvested was paralleled with increases in coniferous timber production, with a significant (P = 0.011, Table 1) breakpoint also detected in the in the 1960s.
Total carbon stocks displayed a similar pattern to broadleaved woodland, showing a steady decline between 1930 and 1980, with a significant breakpoint detected in the 1970s (P = 0.044, Table 1). Overall, we calculated stocks of 22.1 M tonnes in 1930 and 18 M tonnes in 2015 (Fig. 2). By 2015 approximately 51% of the carbon stocks of Dorset were stored in either broadleaved or coniferous woodland (Appendix 2). The distribution of carbon also changed noticeably over the four periods (Appendix 5) especially in areas such as west and north Dorset where areas of high carbon sequestration (e.g. by unimproved grasslands) were converted to land uses containing lower carbon stocks, such as improved grassland or arable land. No breakpoint was detected for flood risk mitigation, which largely mirrored the trends in arable crops, with a pronounced decline from 1930 to 1950, values increasing thereafter following a curvilinear pattern. By 2015, values were similar to those recorded at the outset. Water quality in the form of nitrogen retention and export showed similar trends with values increasing significantly between 1930 and1950, then demonstrating more rapid increases between 1950 and 1980, with a highly significant breakpoint in the 1960s (P = 0.005, Table 1), before declining gradually thereafter. Negative trends between the 1930 and 1980 were also apparent for the ES of soil quality (Fig. 2), with two separate breakpoints detected in the 1960s (P = 0.021, Table 1) and again in the 1980s (P = 0.023, Table 1). Most of the erosion was associated with areas of arable and improved grassland, with existing arable areas becoming more degraded by the 1950's (Appendix 5).
Modelled total annual water yield of the study's catchments equated to a steady 81% decrease in water yield between 1930 (244 Million m 3 /ha) and 1980 (45 Million m 3 /ha), with the greatest change occurring in the 1950s (P = 0.026, Table 1). These trends are most visible in the higher central chalklands of Dorset and the poorer soils to the north and west of Dorchester (Appendix 5). Nutrient loads and retention efficiencies were shown to vary greatly across Dorset's catchments (Fig. 2). For example, the highest retention and export values for nitrogen were consistently recorded in the lower Frome catchment across all four periods, with nitrogen export levels showing an increasing trend in the rivers and groundwater transitioning into Poole Harbour (Appendix 5). Aesthetic value and visitation rates also showed a strong initial decline with a breakpoint in the 1950s (breakpoint P = 0.004, Table 1) and 1960s (breakpoint P = 0.028, Table 1) respectively, followed by a partial increase in vitiation rates by the 1980s (Fig. 2). By the 1980s and 2015 recreation values were generally concentrated in hotspots (Appendix 5) around coastal areas (e.g., The Jurassic Coast World Heritage Site), fragments of woodland and designated heathland nature reserves (e.g., Studland Heath National Nature Reserve). The remaining three ES, representing biodiversity and habitat quality (for BAP species and for pollinators), all showed a broadly similar pattern, with relatively rapid declines observed prior to 1980 and much more gradual declines recorded between 1980 and 2015 (Fig. 2). Significant (P B 0.05) breakpoints were observed for all three of these services between 1950 and 1980 (Table 1).
Of the 75 possible GAM models among all driverresponse combinations, 14 were statistically significant (P B 0.05) with the smoothing function included ( Table 2). The full trajectory of each significant GAM is given in the supplementary information (Appendix  6). The deviance explained by the models that include a significant smoothing function ranged from 58 to 98% (Table 2). Changes in livestock numbers, broadleaved woodland, coniferous woodland, nitrogen export, nitrogen retention and habitat quality for BAP species were all significantly correlated (P \ 0.05) with the increased application of fertilizer. A significant (P \ 0.05) negative non-linear relationship was found between pesticide application and habitat quality for pollinators (Table 2; see also Appendix 6). A significant (P \ 0.05) linear relationship (Appendix 6) was also found between total farming income, aesthetic value of landscapes and arable crop production. The loss of carbon and its relationship to agricultural intensification was reflected by significant negative relationships recorded with both agricultural productivity (TFPA) and number of tractors (Table 2), the latter being a proxy for agricultural mechanisation. Soil quality was also negatively related to this proxy, as were recreation value and habitat quality for BAP species (Table 2). Application of pesticides was a significant driver  increase of 15%. In all cases, these trends were broadly linear (Fig. 3). Twelve of the 15 environmentally attributable sectors increased in GVA value between 1981 and 2015 (Appendix 4). The largest increases ([ 100%) were registered in the 'Chemicals', 'Pharmaceuticals' and 'Recreational services' sectors. 'Agriculture, forestry and fishing', which is a significant industry for Dorset, also increased by 55.2% in terms of GVA during this period. However, this sector was also one of only two industries where employment fell during this period, with values declining by 61%. The gross output of agriculture to other industries also declined, by more than half; by 2015, the value was only 5.3% of total outputs (Appendix 4).

Discussion
The data presented here illustrate trends in multiple ES in a landscape undergoing agricultural intensification. Arable crop yields in the UK increased by a factor of approximately three between 1930 and 2015 (Ritchie and Roser 2013), and overall agricultural productivity increased in a near linear fashion during much of this period (Appendix 3). These productivity gains were associated with increased mechanisation of farming, as illustrated by the rapid increase in the number of tractors that occurred after 1945, reaching a peak in the 1980s. Fertilizer applications also increased steadily after 1940 to reach a maximum in the 1980s, after which they underwent a slight decline. Pesticide use increased markedly after the late 1980s, and is currently near an all-time high (Appendix 3). These trends illustrate how the process of intensification has evolved over time, with sequential increases in mechanisation, fertilizer use and pesticide application. Other associated trends during this period include a tendency for farms to increase in size and become more specialised, and a marked decline in farm labour (Robinson and Sutherland 2002; Appendix 3). The changes in patterns of agricultural land use observed in this study are largely attributable to the shifting provision of capital grants and subsidies, resulting from changes in government policies such as the 1947 Agriculture Act and the Common Agricultural Policy (Robinson and Sutherland 2002). These policy developments led to some changes in farming practice, such as the widespread switch from cultivation of arable crops to livestock husbandry observed here between 1950 and 1980. However, as the proportion of land area devoted to agriculture has remained approximately constant in Dorset during the past 85 years, at around 80%, the changes in land use that have occurred largely reflect a process of intensification of agriculture rather than its expansion.
The land cover changes documented here provide clear evidence of biodiversity loss, as indicated by the substantial area declines of habitats with high conservation value, such as calcareous grassland, unimproved neutral grassland and lowland heathland. The current findings are consistent with those obtained from field monitoring data; Ridding et al. (2020a) recorded losses of 97% of neutral grassland sites and 70% of calcareous grassland sites in Dorset between 1930 and 2015. The pronounced declines in biodiversity indices observed here are also consistent with a broad literature that has documented the negative environmental impacts of agricultural intensification, such as the widespread declines observed in habitat b total employment in Dorset (thick black line) and total employment including only those sectors with strong links to the environment (thin dotted line) calculated using the Cambridge econometrics local economic forecasting model (LEFM) diversity and quality, and declines in populations of particular groups of species such as seed-eating birds (Robinson and Sutherland 2002) and plants adapted to arable habitats (Storkey et al. 2012); Emmerson et al. 2016). The steep decline in habitat for insect pollinators recorded here between 1930 and 1950 closely accords with the analyses presented by Ollerton et al. (2014), who documented a sustained period of extinctions of pollinating insect species in Britain from the late 1920s to the late 1950s. These losses were attributed to agricultural intensification during this period (Ollerton et al. 2014). Negative impacts of agricultural intensification on different groups of species have been attributed to increased application of fertilizers, insecticides and herbicides, together with a range of other factors including loss of food resources and reduction of habitat diversity at the landscape scale (Firbank et al. 2008). In particular, we found confirmation of non-linear thresholds in BAP habitat area and condition linked to agricultural drivers such increased fertilizer rates, agricultural productivity and increased mechanisation. The spatial distributions of these measures also changed markedly, indicating the remaining habitat was more fragmented and spatially segregated into hotspots around designated protected areas of Dorset. The loss of biodiversity and increase in agricultural production reflect the UK-wide trends reported in the UKNEA (2011) and more globally (Chaudhary et al. 2016;Molotoks et al. 2017).
Provision of a number of other ES also declined markedly in Dorset over the past 85 years. Recreation value, soil quality and carbon storage all displayed similar responses to the biodiversity measures, declining in a curvilinear pattern throughout much of the study interval. Other ES were characterised by more complex trajectories; water yield, mitigation of flood risk and aesthetic value each declined sharply to a minimum value, after which some recovery was observed. These patterns are largely attributable to the peak in arable crop cultivation in 1950, which was followed by a subsequent shift from arable to livestock farming in parts of the study area. Cultivation of arable crops is known to increase the risk of flooding owing to reduced water retention capacity of the soil, which can increase run-off (Deasy et al. 2014). This can be attributed to soil compaction caused by the use of heavy machinery; the creation of bare soil from which runoff is unchecked; concentration of overland flow in plough tines and tyre tracks; and rapid transfer of runoff to water courses via vehicle tracks (O'Connell et al. 2007). According to the results of stakeholder surveys, arable cropland also has a relatively low aesthetic value (Newton et al. 2012;Gosal et al. 2018).
Despite these more complex trends, overall, these results provide further evidence of the trade-offs that have been widely documented between provisioning ES, such as food and timber production, and many other ES (Power 2010;Newton et al. 2012;Howe et al. 2014;Holt et al. 2016;Bernués et al. 2019). However, the time-series data presented here indicate that these trade-offs have varied over time. For example, the initial increase in arable crop production  was associated with marked declines in some measures (notably habitat quality for pollinators, aesthetic value, water yield and flood protection) but less pronounced changes in others. Conversely, the subsequent period  witnessed a pronounced growth in livestock farming and a decline in arable cropland, which coincided with a more rapid decrease in one service (carbon storage) but a reduced rate of decline in others; conversely a third group (e.g. flood protection and aesthetic value) displayed some recovery. These findings are relevant to the concept of ES ''bundles'', or sets of ES that repeatedly co-occur across space or time. It has been suggested that identification of such bundles could help inform development of appropriate environmental policies or land management approaches, by enabling tradeoffs between ES to be identified (Raudsepp-Hearne et al. 2010a). While the concept has been widely investigated (Saidi and Spray 2018), the majority of studies have analysed spatial coincidence of services at a single time; as noted by Spake et al. (2017), this cannot provide a mechanistic understanding of ES dynamics and relationships. For this, time-series data are required, such as those presented here. The current results demonstrate that the relationships between ES can change over time, highlighting the potential pitfalls associated with identifying ES bundles at a single time point.
Relatively few other studies have examined the temporal dynamics of ES provision ; examples include Renard et al. (2015) in Ncube et al. (2018) in southern Scotland. In the former study, spanning 35 years, Renard et al. (2015) found that provision of most ES increased through time, although the rate of increase varied spatially. As in the current study, changes in agricultural policy were identified as a key driver of the observed trends. However, while trade-offs between provisioning and cultural ES were identified in 1971, some of these had shifted to either a neutral or synergistic relationship by 2006 (Renard et al. 2015), in contrast to the current results. Developments in agricultural and forest policy were similarly found to be the principal cause of changes in ES provision observed in southern Scotland over a 60-year interval (1946-2009). Key changes during this period included a substantial increase in livestock and timber production and a decline in arable crop production, which were accompanied by declines in most regulating ES as well as biodiversity (Ncube et al. 2018). These general trends and the trade-offs observed closely parallel those observed here, although in contrast to the current results, an increase in both carbon storage and flood control were observed in the Scottish study. Collectively these studies demonstrate the value of analysing trends over time for detailing the changing relationships between ES, and for understanding the mechanisms underlying ES dynamics, supporting suggestions made previously (Bürgi et al. 2015;Dallimer et al. 2015;).
Here we build upon the previous research into temporal dynamics by examining the occurrence of thresholds in ES provision. Breakpoint analysis indicated non-linear trends in 13 of the 15 ES considered here, suggesting that thresholds in ES provision may be widespread. Sutherland et al. (2016) similarly documented non-linear trajectories of all ES examined in a Canadian forest following logging. These findings are consistent with theoretical expectations (Bullock et al. 2011) and meta-analyses of empirical data (Martin et al. 2013;Spake et al. 2015) indicating that ES often display non-linear trajectories of recovery following disturbance events. However, less evidence is available regarding ecosystems undergoing continuous degradation or decline, such as those examined here. Peng et al. (2017) provide an example describing the effects of urbanization in China; above thresholds of population density and urbanisation intensity, total ES provision was found to decline rapidly. In a similar way, our results suggest that continuing agricultural intensification may result in thresholds in ES provision, characterised by relatively abrupt shifts in the rate of change. Such phenomena could undermine efforts to develop environmentally sustainable approaches to agricultural intensification (Rockström et al. 2017;Pretty et al. 2018).
Might these thresholds represent tipping points? This depends critically on how tipping points are defined. Here we follow van Nes et al. (2016) and Milkoreit et al. (2018) in considering tipping points to be ecological thresholds that are driven by positive feedback loops. As noted by Dallimer et al. (2015), such feedbacks have largely been ignored by the literature on ES, which has major implications for our understanding of how ES delivery will alter in future as a result of land-use/cover change. This knowledge gap partly reflects the difficulty of identifying potential feedback mechanisms, and evaluating their impacts on ES provision. Some insights into these mechanisms are provided by our analysis of the relationships between changes in the provision of ES and the underlying drivers. Most of these relationships were linear, indicating that variation in the ES appeared to track variation in the driver (see Appendix 6). While some departures from linearity were observed in a minority of cases, these could have been caused by interactions between drivers, or from an abrupt change in the state of the ecosystem with a small change in a driver, rather than positive feedback loops (Andersen et al. 2009;Watson et al. 2018).
The current results therefore provide preliminary evidence that these thresholds could be considered as tipping points. However, further research into positive feedback mechanisms is required to fully determine the occurrence of tipping points in agricultural systems. This reflects limitations in the data, both in terms of the temporal frequency of sample points, and the resolution of the driver data. Potential feedback mechanisms can be hypothesized from our GAM analysis; for example, it is conceivable that progressive soil degradation caused by intensification of crop cultivation could lead to increasing amounts of fertilizer being applied, in order to maintain crop productivity. In particular, intensification of agriculture is known to reduce soil organic carbon (SOC), which can lead to losses in soil structure, waterholding capacity and filtration, leading to a decline in crop yields (Edmondson et al. 2014;Haygarth and Ritz 2009). To mitigate such losses, farmers might potentially apply progressively more fertilizer, but information on this behaviour is lacking. Positive feedbacks could also account for the rapid declines occurring in many elements of biodiversity associated with agricultural intensification, as a result of cascading extinctions and the consequent collapse of food webs (Newton 2021). Other possible tipping point mechanisms include possible feedbacks between pesticide uses and outcomes (Sponsler et al. 2019), the over-saturation of soil nutrients leading to accelerated leaching to water courses, and the local extinction of Rhizobium bacteria as a result of sewage sludge applications (Haygarth and Ritz 2009). Moreover, information is lacking on all of the positive feedback mechanisms associated with agricultural land use, especially those that link both the environmental and economic elements of the socio-ecological system (Benton et al. 2017). Their role in driving tipping points is therefore uncertain, and merits further research attention.
Other imitations of this investigation include the use of spatial modelling tools and proxies to assess ES provision. While rapid progress has recently been made in developing methods for assessing ES values (Costanza et al. 2017), there is still considerable debate regarding the relative effectiveness of different approaches. Although InVEST is one of the most widely used tools for mapping ES, and is considered to be relatively robust (Bagstad et al. 2013;Vorstius and Spray 2015), we recommend the following uncertainty assessment analyses: exploration of alternative input datasets for the study region, sensitivity analyses of all model input parameters (e.g. historical crop yields, nutrient input efficiencies, carbon storage estimates and water yield coefficients) and a thorough exploration of the model outputs before using them to inform decisions. This reflects the recommendations of previous InVEST studies across a number of ES (e.g. Sánchez-Canales et al. 2012;Sharps et al. 2017;Redhead et al. 2018;Bagstad et al. 2018). Although the use of proxies such as land cover for assessing ES is also very widespread, this similarly has limitations (Eigenbrod et al. 2010). Development of appropriate methods for assessment of cultural ES has proved to be particularly challenging; it has been argued that current approaches fail to capture the full value of places to people's lives and provision of such services may not be directly linked to features of ecosystems that can be readily measured (Gosal et al. 2018). For example, the InVEST (v3.4.4) Visitation: Recreation and Tourism model used in this study might have overestimated visitor rates before 2015, owing to a possible model constraint associated with only using Flickr uploads as a data source. Additional linkages with other traditional sources of data, such as recreation surveys, or official tourism data could be used in future to improve model outputs. Given these limitations, the historical cultural ES results should clearly be viewed with caution. Another caveat is that we analysed ES change using only four time points in time prior to time-series interpolation. Simply, no similar land cover maps were generally available for the 1940, 1960-1970or 1990-2010periods. However, Ridding et al. (2020b has recently synthesised maps for the 1990s and showed that some habitats such as improved grassland had a roughly linear trend in Dorset over the last century. Nonetheless, patterns of semi-natural habitats have fluctuated in a non-linear fashion over the study period and our snapshots do not fully capture these temporal subtilties or how this may have affected ES. Despite these caveats, the current results demonstrate that the process of agricultural intensification can lead to relatively abrupt changes in the provision of ES. This is consistent with suggestions that agricultural land use could be associated with tipping points, which have the potential cause food shortfalls and price spikes (Benton et al. 2017). To test this hypothesis more fully, a deeper understanding is required of feedback processes occurring in both the environmental and socio-economic components of agricultural systems, together with interactions between the two (Benton et al. 2017). The economic data included here provide some indication of these linkages. Strikingly, despite the fluctuations in provision of different ES, overall economic activity increased almost linearly during the study interval, in line with the increase in agricultural productivity. This relates to the paradox identified by Raudsepp-Hearne et al. (2010b), namely that human well-being (which is linked to economic development) has increased in many areas, while provision of many ES has declined. This is the converse of what would be expected if human well-being is dependent on ES, as many researchers have posited (Costanza et al. 2017). Raudsepp-Hearne et al. (2010b) attribute this paradox to the fact that technology and social innovation have decoupled human well-being from ecosystem degradation; technological advancement has enabled humans to exploit some ES at the expense of others. Agricultural intensification arguably provides a leading example of this process. Efforts at achieving sustainable intensification of agriculture (Rockström et al. 2017;Pretty et al. 2018) will therefore need to address this trade-off if they are to be successful.