Land use change drives the spatio-temporal variation of ecosystem services and their interactions along an altitudinal gradient in Brazil

Land use and land cover (LULC) changes may affect the provision of ecosystem services. However, little is known how LULC changes influence the spatio-temporal variation in ecosystem service and their interactions along altitudinal gradients. Here we assessed the spatio-temporal variation of eight ecosystem services in an altitudinal gradient between the year of 1986 and 2015, and quantified the effect of LULC transitions on the provision and interactions of ecosystems services. We modelled and mapped eight ecosystem services in an altitudinal gradient characterized by low (< 600 m), middle (600–1200) and high altitudes (> 1200) in Zona da Mata, Minas Gerais state, Brazil. We quantified changes in ecosystem services by contrasting ecosystem service maps between 1986 and 2015, and explored how four common LULC transitions affected the variation and the interactions between the eight ecosystem services. The spatio-temporal analysis indicated that six out of eight ecosystem services increased from 1986 to 2015, while soil erosion control and water flow regulation decreased. In areas above 1200 m, regulating services dominated, while in areas below 1200 m provisioning service were most evidenced. LULC transitions from forest to agricultural areas, and vice versa, resulted in trade-offs between provisioning and regulating ecosystem services. LULC changes drive the spatio-temporal variation of ecosystem services along an altitudinal gradient with contrasting biophysical conditions. Future management of ecosystem services in the landscapes should take into the account the biophysical conditions and the consequences of specific LULC transitions.


Introduction
Human-modified landscapes are complex socio-ecological systems, and in the face of global changes, the future of human wellbeing will be shaped by our capability to manage these complex environments to provide multiple ecosystem services that meet societal demands (Reyers et al. 2013;Vallés-Planells et al. 2014;WWDR 2018). Although ecosystem services research has made significant advances in the conceptualization, quantification and monetary valuation of ecosystem services in the last decades (Costanza et al. 2014;De Groot et al. 2012;Yi et al. 2018), we still have limited understanding of how provision patterns of multiple ecosystem services emerge and change in space and over time (Renard et al. 2015). Spatiotemporal patterns of ecosystem services can be influenced by a range of abiotic and biotic factors, as well socio-economic factors (Pan et al. 1999;La Notte et al. 2017;Mayor et al. 2017). Biophysical conditions, such as temperature, precipitation, evapotranspiration, soil type and ultimately vegetation types, can vary sharply along altitudinal gradients. As a consequence, mountainous areas may reflect a mosaic of contrasting ecosystems and land use types at relatively short distances, ultimately leading to heterogenous spatial patterns of ecosystem services (Dieleman et al. 2013;Körner 2003). Therefore, mountainous areas offer a unique opportunity to study how provision patterns of ecosystem services change over time in heterogenous environments with contrasting biophysical conditions.
With a need to increase global food production by 60% by 2050 to support a growing world population (PRB 2018;WWDR 2018) intense changes in LULC can be anticipated for the near future. This poses a challenge for farmers, scientists and policy makers to develop agricultural systems that produce sufficient and nutritious food, while also delivering other essential ecosystem services. Yet, our understanding of the way in which changes in LULC influence provision levels of different ecosystem services is still limited. In many parts of the world, agricultural policies and land managers have focussed on strengthening provisioning services, often resulting in a decrease of other ecosystem services, such as climate regulation and fresh water supply (MA 2005;Butchart et al. 2010). In Brazil, the deforestation in the Amazon biome for food production and pasture land has decreased the capacity of ecosystems to sequester carbon, nutrient cycling, erosion control and water regulation (Foley et al. 2007;Portela and Rademacher 2001). Over the time, the loss in regulating services can have impacts at the local scale, such as reduced water supply, as well as at the national and global scale, for instance impacts on climate regulation. Understanding how LULC changes influence the associated provision of ecosystem services can provide important insights that are relevant for the sustainable land management.
Contrasting LULC transitions may have specific impacts on ecosystem services, which may unfold as synergies (win-win scenarios), trade-offs (win-lose scenarios), or, in case of lose-lose scenarios, as dissynergies (Bennett et al. 2009). While interactions between ecosystem services have been reported at different spatial and temporal scales (Briner et al. 2013;Lang and Song 2018;Li et al. 2018;Qiao et al. 2019;Qiu et al. 2018), the consequences of LULC transitions for interactions between multiple ecosystem services are still relatively poorly understood (Valujeva et al. 2016). The assessment of interactions between multiple ecosystem services is challenging due to the context dependency and the overwhelming complexity of several possible pairwise interactions (Raudsepp-Hearne et al. 2010;Maes et al. 2012;Howe et al. 2014). Analysing the effect of specific LULC transitions on the provisioning levels of ecosystem services may offer a pathway to explore the interactions between multiple ecosystem services, and to identify management actions that lead to synergies between ecosystem services rather than trade-offs. Quantitative information about the effect of LULC transitions on ecosystem services and on their interactions can guide spatial planning, land management and policy making to better manage the provision ecosystem services in the future.
In Brazil, where intense LULC changes have been taking place and still are on-going, information about LULC changes and ecosystem services can inform future land use planning. The Brazilian Atlantic Forest biome supports 70% of the Brazilian population and is considered the 5th hotspot of biodiversity in the world (Myers et al. 2000). However, due to large-scale deforestation the forest cover in the Atlantic Forest biome has been reduced to only 12.4% of the original forest area (Sosma 2019). In the southeast region of Atlantic Forest biome, the Zona da Mata region consists of heterogeneous landscapes with contrasting climatic conditions and has witnessed intense changes in LULC in the past (Gomes et al. 2020). This region contains the Serra do Brigadeiro state park and Caparaó National park located at higher altitudes, which provide a refuge area for several endemic plant and animal species, and offers opportunities for tourism and recreation (Moura et al. 2012;Machado et al. 2016). These aspects make this region a suitable case study area to analyse the effect of LULC changes on the spatio-temporal variation of ecosystems services in heterogeneous landscapes. The objective of the study is to analyse the spatio-temporal variation of eight ecosystem services in an altitudinal gradient in the southeast region of the Atlantic Forest biome in Brazil. More specifically, we (i) assess the spatio-temporal provision of eight ecosystem services, including cultural services, coffee production, livestock production, carbon storage, water flow regulation, soil erosion control, pollination and habitat quality, from 1986 to 2015, (ii) assess how the spatio-temporal provision of ecosystem services varies along an altitudinal gradient representing contrasting biophysical conditions, and (iii) assess how specific LULC transitions affect ecosystem services and their interactions.

Study area
The study area is located in the Zona da Mata of Minas Gerais state and covers an area of 11.300 km 2 in the Brazilian Atlantic Forest biome (Fig. 1). This region can be considered as a complex socio-ecological system with the predominance of family farmers, and have been subject to LULC changes in the last three decades (Cardoso et al. 2001;Jackson et al. 2012;Gomes et al. 2020). The study region has an altitudinal gradient ranging from 27 m to almost 2900 m above sea level. This altitudinal gradient gives rise to a heterogeneous landscape mosaic with strong gradients in temperature, precipitation, geomorphology and soil type at short distances ( Fig. 1, Online Appendix 1). According to the gradient of temperature and precipitation (Appendix 1) we distinguished three main altitudinal zones in the study region: Low (\ 600 m), Middle (600-1200m) and High altitudes ([ 1200 m). Altitudes higher than 600 m offer optimal climatic conditions for coffee production, a key cash crop in the area (Alègre, 1959;Fick and Hijmans, 2017;Souza et al. 2012), and protected nature conservation areas prevail at 1200 m and higher. The Low altitude zone has a mean annual precipitation of 1242 mm, a mean annual temperature of 22°C and a potential evapotranspiration of 1515 mm. The main soil type is an Acrisol with a low soil fertility and accumulation of clay around 50 cm of the soil depth, which makes this soil unsuitable for crops with a deep rooting system. The Low altitude zone is characterized by large plains, where the main rural activity is dairy farming and cattle ranging, and the population density in urban area increased from 20 to 62 persons per km 2 from 1980 to 2010. The Middle altitude zone has an annual precipitation of 1333 mm, a mean annual temperature of 19°C and potential evapotranspiration of 1410 mm. The main soil type is a Ferrassol, representing deep and weathered soils with a low nutrient availability, and a granular soil structure that is favourable for water infiltration and retention. The Middle zone is characterized by the predominance of hills and the dominant land use is coffee production and cattle ranging. The population density in urban areas increased from 15 to 65 persons per km 2 from 1980 to 2010. The High altitude zone has a mean annual precipitation of 1510 mm, a mean annual temperature of 168C and a potential evapotranspiration of 1245 mm. The main soil types are Cambisols and Regosols, which are shallow, but can have a high nutrient availability. The geomorphology is characterized by sharp hills with rock formations covered with Campos Rupestres (grassy-scrub vegetation) at the mountain tops. This area is not permanently inhabited.

LULC transitions
In this study we used the LULC maps from 1986 to 2015 developed by Gomes (2020, Online Appendix 2). The LULC maps were developed using the images from Landsat satellite 5 and 8 and classified using the Random Forest algorithm with an average performance of 0.90 for kappa and 0.93 for accuracy. In contrast to many other tropical areas, the study region witnessed an increase in forest cover from 18 to 24%, an increase in the area of coffee plantations from 3 to 11%, and a decrease of pasture area from 76 to 58% from 1986 to 2015. These changes have been fostered by governmental supported investments in agriculture and protection of the environment (Gomes et al. 2020). We focussed on LULC transitions associated with pasture, forest, coffee and urban area because these four LULC types represent approximately 95% of the area in the study region. Using the maps of LULC from 1986 to 2015, we derived the four most prevalent LULC transitions: ''pasture to forest'', ''forest to coffee'', ''pasture to coffee'' and ''pasture to urban area'' (Fig. 2).

Ecosystem services and InVEST
We focussed on eight ecosystem services: cultural services, coffee production, livestock production, carbon storage, water flow regulation, soil erosion control, pollination and habitat quality. These ecosystem services were chosen because these are considered important by the local population (Teixeira et al. 2018) and relevant data sources were available. More specifically, we used nature conservation areas as an indicator of cultural services because these reflect opportunities for tourism, recreation and natural history education. Coffee production and livestock production were selected as these represent major cash generating activities for farmers in the region (Cardoso et al. 2001;Souza et al. 2012). Water regulation was selected since it was identified as the main ecosystem service perceived by farmers in the region (Teixeira et al. 2018). We selected erosion control because it is a common source of soil degradation in this mountainous region, and pollination and habitat quality were included as proxies for biodiversity. Using the LULC maps of 1986 and 2015 we assessed the cultural ecosystem services, coffee production and livestock. Cultural service provision was based on the area of natural conservation areas (ha), and coffee production was informed by the area of coffee plantations (ha). To assess livestock production, we used cattle stocking rate data (animals/ha) obtained from the Brazilian bureau of statistics from 1986 to 2015 (IBGE, 2018) for each municipality in the study area and combined these with the pasture areas from the LULC maps from 1986 to 2015. Carbon storage (Mg), water flow regulation (mm), soil erosion control (tons ha -1 year -1 ), pollination and habitat quality (index from 0 to 1) were assessed using the Integrated Valuation and of Ecosystem Services and Trade-offs (InVEST) model. The InVEST model comprises a suite of free software models developed by Natural Capital project to map and value the ecosystem services (Nelson et al. 2009;Sharp et al. 2018). InVEST uses maps and information of LULC types as input and generates maps of ecosystem services. The quantitative and spatially explicit outputs can inform management and policy making by, for instance, evaluating alternative management options and assess the likely impacts of human activities on ecosystem services. We used LULC maps from 1986 to 2015 at 30 9 30 m resolution and associated information of the land use types as model input (Online Appendix 1) and the InVEST outputs consisted of maps of the respective ecosystem services. To enable the comparison of the output maps of the various ecosystem services, output maps were normalized, generating maps with values ranging from 0 to 1 for all ecosystem services.

Carbon storage
Carbon storage was assessed using the LULC maps and the carbon density per LULC class, considering four carbon pools: above ground biomass, below ground biomass, soil and dead organic matter (Online Appendix 3).
where carbon storage CA x is the carbon stored in each pixel x (Mg). PA is the pixel area (30 9 30 m = 900 m 2 or 0.09 ha), C A is the aboveground carbon density (Mg ha -1 ); C B is the belowground carbon density (Mg ha -1 ); C S is the soil carbon density (Mg ha -1 ); and C D is the dead mass carbon density (Mg ha -1 ). We used carbon pool data from literature based on local studies for all LULC types (Online Appendix 3), and the output consisted of a map of carbon storage expressed as Mg of carbon per pixel.

Water flow regulation
The water flow regulation was expressed as the water yield index, which is defined as the amount of water that runs off from each pixel in the landscape (Tallis 2011). The average annual water yield was calculated using a water balance based on precipitation and evapotranspiration data: where Y x is the average annual water yield, P x is the average annual precipitation (mm), and AET x is the annual actual evapotranspiration (mm) in pixel x (Sharp et al. 2018). We used spatial data of average annual precipitation of the study area from 1970 to 2000 from the WorldClim database (Fick and Hijmans 2017), the reference annual evapotranspiration from Dias (2018) and the soil depth data from the Brazilian soil database (Cooper et al. 2005). The output map indicated water yield expressed in mm per pixel.

Soil erosion control
The soil erosion control was assessed by the average annual rate of soil loss (ARSL), which was calculated using the revised Universal Soil Loss Equation for each pixel x (Sonneveld and Nearing, 2003) (Eq. 3).
where A x is the annual rate of soil loss (tons ha -1year -1 ) in pixel x, R x is the rainfall erosivity (MJ mm (ha h) -1 ); K x is the soil erodibility factor (ton ha hr (MJ ha mm) -1 ); LS x is the slope gradient (dimensionless); C x is the crop management factor (dimensionless); and P x is the support practice factor (dimensionless). The digital elevation model (DEM) was obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) dataset with a resolution of 30 9 30 m (http://earthexplorer.usgs.gov/). We obtained the rainfall erosivity parameter R x using a multivariate equation based in altitude, longitude and latitude developed for Brazilian territory (Mello et al. 2013). The soil erodibility parameter K x was based in the soil class, and the LS x factor using the equation from Moore and Burch (1986). The crop management factor C x and the support practice P x for each LULC class were derived from the literature (Online Appendix 4, 5). The output consisted of a map of soil erosion expressed in tons ha -1 year -1 . For the temporal analysis, we considered the inverse of relative changes of soil erosion and referred to this as the soil erosion control index.

Pollination
Coffee production is an important cash crop in the study region and there is a diversity of wild bees that pollinate the coffee plants, especially in agroforestry coffee systems (Ferreira 2008). Here, we focus on the unmanaged honey bee Apis mellifera and the wild bee Trigona spinapis, the most abundant bee species in coffee plantations (Ferreira 2008;Malerbo-Souza and Halak 2012). We assessed the potential pollination services provided by these two bee species using the InVEST pollination model. The procedure entails the scoring of land use types for their potential to support floral resources and nesting sites and generates an index of the relative abundance of pollinators. Species specific estimates for foraging distance, habitat for nesting and foraging likelihood in each land cover class were obtained from literature (Online Appendix 6, 7).

Habitat quality
Habitat quality refers to the ability of ecosystems to provide suitable conditions for species population viability in terms of anthropogenic threat/disturbance levels (Yan et al. 2018). We considered four threats: coffee fields, pastures, urban areas, and paved roads (Table A.6; Duarte et al. 2016). Data on the location of paved roads of the study area were obtained from IBGE website (IBGE 2019). Habitat quality was assessed based on the relative impact of threats, the sensitivity of the habitat to threats, and the distance between habitats and location of threats. The impact of a threat decreases with increasing distance from the location of the threat, and an impact map is generated by integrating the impact zones around the land use types considered as a threat. Habitat quality is then derived as maps with a relative metric ranging between 0 and 1, with low values for high impact zones and high values for low impact zones (Online Appendix 8).

Interactions between ecosystem services
We analysed the interactions for pairwise combinations of ecosystem services in the study region from 1986 (T1) to 2015 (T2) (Haase et al. 2012; (Fig. 3). First, we created a map of temporal changes (DES j ) for each ecosystem service j in 1986 and 2015 (DES j = ES1 j,T2 -ES1 j,T1 ) using the raster calculator in ArcGIS 10 (Fig. 3a). Then, we normalized the DES maps values generating a new map (NDES j ) with values ranging between -1 and 1 (Fig. 3b). Next, we clipped the NDES j using a map mask of Land Use Transitions (LUT j ) that contained only the pixels that were converted from specific transitions between 1986 and 2015 (Fig. 3c). This resulted in a map (ES j,LUT ) representing the NDES j pixels values only for the areas where LULC transitions (LUT) took place between 1986 and 2015 (Fig. 3d). Finally, we visualised the interactions between ecosystem services by plotting the mean values of the map of ecosystem services (lESj LUT ) for each pairwise combination (Fig. 3e). We expressed the interactions between pairwise ecosystem services in terms of synergies (win-win), trade-offs (lose-win; win-lose) or dis-synergies (lose-lose) (Bennett et al. 2009;Haase et al. 2012). This procedure was followed for each of the four LULC transitions: ''pasture to forest'', ''forest to coffee'', ''pasture to coffee'' and ''pasture to urban area''.

Results
Across the whole study region, six out of eight ecosystem services increased between 1986 and 2015: cultural services (? 33.3%), coffee production (? 266.6%), habitat quality (? 40%), pollination (? 29.2%), livestock stocking rates (? 26.3%) and carbon storage (? 1.8%), while water flow regulation (-2.2%) and soil erosion control (-3.7%; Fig. 4) decreased. However, the three altitude zones showed distinct spatio-temporal patterns of ecosystem services (Fig. 5). While the livestock stocking rate in the Low zone increased from 0.75 to 1.24 animals/ha between 1986 and 2015, the overall herd size in the study region decreased by 10% because pasture areas decreased from 76% in 1986 to 58% in 2015. On the other hand, the Low Zone showed a strong increase in habitat quality (? 76%), water flow regulation (? 60%), pollination (? 54%) and soil erosion control (? 27%) between 1986 and 2015. In the Middle zone the area for coffee production increased more than 260%, followed by an increase of habitat quality (? 75%) and pollination (? 45%) and a decrease in soil erosion control (-3.4%). In the High zone there was an increase in water yield from 0.3 to 1986 to 0.34 in 2015 (normalized values), and cultural services increased by 33.3% due to the establishment of the Serra do Brigadeiro state park. Habitat quality (? 8.5%) and carbon storage (? 0.4%) increased as well.
The LULC transitions affected the provision of ecosystem services and their interactions in different ways (Fig. 6). The conversion of pasture to forest lead to an estimated 100% decrease of livestock production, but a 516% increase of pollination services, an approximate 98% increase of habitat quality and water flow regulation, and a 55% increase of carbon storage (Fig. 6a). This LULC transition lead to joint increases for ten out of fourteen pairwise comparisons, with strong synergies for carbon storage and habitat quality (0.55 and 0.98, respectively) and a weaker synergy of water flow regulation and soil erosion control (0.04 and 0.002, respectively; Fig. 6b). The conversion of forest to coffee had a positive impact on coffee production (? 100%), but led to reduction of soil erosion control (-3229%), water flow regulation (-397%), habitat quality (-100%), carbon storage (-88%), and pollination services (-79%; Fig. 6a). This conversion led to ten dis-synergies and five trade-off responses between ecosystem services, with strong dis-synergies interactions of carbon storage and habitat quality (-0.88 and -0.98), carbon storage and pollination (-0.88 and -0.51), and habitat quality and pollination (-0.98 and -0.51; Fig. 6b). The conversion of pasture to coffee led to a 95% and 100% increase in water flow regulation and coffee production, respectively, and had a negative impact on livestock production (-100%), soil erosion control (-237%) and carbon storage (-32%; Fig. 6a). This conversion led mostly to weak dis-synergy and tradeoff interactions (Fig. 6b). The conversion of pasture to urban area had a limited impact on most ecosystem services, but a strong negative impact on water flow regulation (-3254%), livestock production (-100%) and carbon storage (-43.7%), and a positive impact on soil erosion control (? 83%). Fig. 3 Schematic overview of the methodological approach to assess the impact of specific land use land cover transitions on multiple ecosystem services and their interactions. DESj: map of the changes of ecosystem service j from 1986 to 2015 (a); NDESj: map with normalized values of the changes of ecosystem service j (b); ESj LUT : map of ecosystem service j containing only the values for the pixels from specific LULC transitions (e.g., pasture to forest) (c); lESj LUT : the mean pixel value from the ESj LUT map (d); and visualisation of interactions between pairwise ecosystem services in synergy, trade-off and dis-synergy domains (e) Fig. 4 Maps of the spatio-temporal variation of eight ecosystem services (carbon, habitat, pollination, water flow regulation (water yield), soil erosion control, coffee production, livestock production and cultural services) from 1986 to 2015 in a study region in the Atlantic Forest Biome, Brazil

Discussion
In this study we assessed the spatio-temporal variation of ecosystem services and the impact of four LULC transitions on eight ecosystem services and their interactions from 1986 to 2015 in a region of the Atlantic Forest biome, Brazil. Key findings of our study are that (i) overall, the provision of six out of eight ecosystem services increased in the study area, (ii) the spatio-temporal variation of ecosystem services showed contrasting responses in the three altitude zones due to different biophysical conditions, and (iii) conversion of forest to coffee or pasture has strong negative impacts on erosion control and water flow regulation, with mostly trade-offs and dissynergies between ecosystem services. In contrast, conversion of pasture or coffee to forest has a positive impact on most ecosystem services, except livestock production and coffee production, with the predominance of synergies between regulating ecosystem services.
Our analysis indicates that between 1986 and 2015 there was an increase in coffee production (? 266,6%), habitat quality (40%), pollination (29.2%) and carbon storage (1.8%) in the study region, which were associated with increases in forest cover and coffee production. Worldwide, the expansion of agricultural land for food, feed and raw material production in the last five decades has been accompanied by a reduction up to 40% in forest areas, and has resulted in a decrease of the provision of regulating ecosystem services (Liu et al. 2017;Rukundo et al. 2018). Our study area reflects a contrasting trend with the area under forest cover increasing from 18 to 24% and coffee production increasing from 3 to 11%, resulting in a joint increase in provisioning and many regulation services (Figs. 4,5). The joint increase in coffee and forest areas is mainly due to public policies and environmental legislation, which increased the financial support for smallholder farmers to cultivate coffee and the surveillance measures to protect forest areas in the last three decades (Gomes et al. 2020). In parallel, the state government, NGO's and farmer unions created the Serra do Brigadeiro State Park in 1996, a protected nature conservation area of about 11,500 ha, which is used for recreation and educational programs for natural history and environmental protection. Indeed, public policies are considered important drivers of LULC changes and the associated provision of ecosystem services (Guerra et al. 2016;Rounsevell et al. 2012). The analysis of historical socio-economic trajectories, LULC changes and spatio-temporal provision of ecosystem services may generate new insights in how public policies and socio-economic drivers may influence ecosystem service delivery levels via LULC changes. As such, this approach may inform governmental/non-governmental actions to strengthen ecosystem services.
We found that the contrasting biophysical conditions along the altitudinal gradient gave rise to distinct spatio-temporal patterns of ecosystem service provision in the study region. The provision levels of regulating and cultural services, such as carbon storage, pollination, habitat quality and tourism, were higher in the High zone ([ 1200 m) than at lower altitudes. These findings are in line with other studies. Worldwide, mountainous regions are important for water supply, nature conservation and widely used for tourism (Beniston 2003;Catalan et al. 2017). In China, Gao et al. (2018) found that higher altitudinal zones mainly provided regulating and cultural services, while the lower altitude zones supported provisioning services. In the Alps provisioning ecosystem services were the highest in the valleys that provide suitable conditions for agricultural production, and regulating services were the highest in forested areas which are dominant at higher elevations (Egarter Vigl et al. 2017). The suitability of mountains for regulating services areas can in part be explained by their biophysical conditions, such as steep slopes, low temperatures, and shallow soil layers, which make these areas less suitable for agricultural activities.
Our study also revealed that ecosystems service provision levels are less dynamic in the High zone than in lower zones, and that mountain regions can be expected to remain as sources of regulating ecosystem services in the future. In contrast, in the Middle zone there have been strong dynamics in provisioning and regulating ecosystem services due to intense LULC changes in the last three decades. This zone has favorable biophysical conditions to support agriculture, such as deep weathered soils with a high water holding capacity and excellent climatic conditions for coffee production. Areas with favorable biophysical conditions for agricultural production are the most susceptible ecosystems for conversion to agricultural land (Ramankutty et al. 2002), and the associated changes in ecosystem services. The Low zone is mainly used for livestock raising and are not suitable for coffee production. Nevertheless, in the last decades this zone has seen an increase of pollination, habitat quality and carbon storage, which is associated with the increase of forest patches. Therefore, the higher provision level and lower temporal variability of regulating services in the High zone compared with the Low and Middle zones, highlights that the contrasting biophysical conditions along an altitudinal gradient affects LULC transitions and can be considered key determinants to influence the spatio-temporal provision of ecosystem services. While variation in the spatio-temporal distribution of ecosystem services has been associated with differences in socio-ecological systems at municipality level (Raudsepp-Hearne et al. 2010;Andersson et al. 2015;Queiroz et al. 2015), here we show that biophysical conditions can also influence the spatio-temporal distribution of ecosystem services at smaller spatial scales, either directly (e.g. less pollination at low temperatures) or indirectly through LULC changes. Accounting for biophysical conditions and associated LULC changes may be useful to better understand how ecosystem services may develop in the future under scenarios of climate change.
The exploration of the relationship between LULC change and ecosystem service provision levels showed how specific LULC transitions affect the ecosystem services and their interactions. For instance, the conversion of pasture to forest may result in a 55% increase in carbon storage, while converting forest to coffee may lead to declines in habitat quality (-100%), carbon storage (-88%) and pollination (-79%; Fig. 6). These results are in accordance with previous studies showing that forest areas support more regulating services, while agricultural areas deliver more provisioning services (West et al. 2010;Baral et al. 2013). The conversion of pasture to coffee has increased the provision of provisioning and regulating services, but soil erosion is still a challenge in full sun in coffee cultivation. In the Zona da Mata the coffee plantations are mostly conventionally managed unshaded coffee systems, which are prone to soil erosion, and the intensive insecticide use to control pests might undermine essential pollination services by wild and managed pollinators (Goulson et al. 2015). In contrast, agroforestry coffee systems are less prone to water runoff and soil erosion (Franco et al. 2002;Cannavo et al. 2011) and have superior natural pest suppression (Rezende et al. 2014). Maintaining trees and adopting terraces can also improve water infiltration and decrease soil erosion in the degraded pastures of the Zona da Mata region (Chen et al. 2017;Lopes et al. 2019). While the effect of LULC changes on ecosystem services has been studied based on the association between temporal changes in the LULC and ecosystem services provision (Rodríguez-Echeverry et al. 2018), our analysis extends these findings by showing the consequences of four common LULC transitions on the complex interactions between ecosystem services. For instance, the conversion of forest to coffee areas had negative impacts on a suite of ecosystem services, with strong dis-synergies between carbon storage, habitat quality and pollination. Identifying the intensity of interactions between multiple ecosystem services can help to design and manage landscapes to provide a balanced set of ecosystem services (Gong et al. 2019). Anticipating the intensity of interactions between ecosystem services due to LULC transitions enables to better plan the landscape configuration to avoid unwanted and irreversible loss in ecosystem services. Earlier studies on trade-offs and synergies between ecosystem services mainly focused on the spatial or temporal variation of these interactions Sun et al. 2018). Analysing the quantitative effect of LULC transition on the ecosystem services and their interactions can be used to manage landscapes to achieve desired levels of ecosystem services in the future. intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.