Assessing the potential delivery of ecosystem services by farmlands under contrasting management intensities

Farming systems under contrasting management practices can contribute differently to the delivery of bundles of ecosystem services (ES) in agricultural landscapes. Low intensity farming systems, such as High Nature Value farmlands, are expected to deliver a wider range of ES, whereas landscapes under more intensive management are expected to deliver mainly provisioning services. Understanding the management practices associated with desirable outcomes in terms of biodiversity and ES in agricultural landscapes is needed. Our research aimed to understand the links between the delivery of ES bundles associated with agricultural landscapes, and their socio-ecological drivers, using a region in northern Portugal as a case study. Based on publicly available data on ecosystems services and drivers, we analyzed ES associations, delineated ES bundles, and investigated their relationship with socioecological drivers. Overall, our results suggested spatial trade-offs between landscapes delivering provisioning services of high economic value, and landscapes delivering a more balanced set of multiple ES. Bundle analysis highlighted an association between higher landscape multifunctionality and higher values of landscape complexity, higher number of farmers, and farm sizes. Our results reflected the complexity of social and ecological factors operating at the landscape level, pinpointed landscapes with higher multifunctionality and disclosed the conditions underlying their occurrence. The results also highlighted the importance of low-intensity farming systems, namely those supporting High Nature Value farmlands, for the delivery of a wider range of ES at the landscape scale.


INTRODUCTION
Farming systems occurring in agricultural landscapes drive the delivery of key ecosystem services (ES) e.g., food provision, pollination, or touristic value (Bignal and McCracken 1996, DeClerck et al. 2016, Wood et al. 2018. Farming systems result from farmers' choices, biophysical and socioeconomic factors, that altogether affect the management practices and ultimately shape agricultural landscapes (Plieninger et al. 2016, Santos et al. 2021. Intensification of farming systems has caused an increase in farm size and landscape homogenization due to specialization, i.e., the production of fewer products (either animal or crop production) with increasing use of machinery and external inputs (Emmerson et al. 2016). Changes in farming systems and associated social-ecological drivers influence the potential of agricultural landscapes to deliver different ecosystem services. Consequently, intensively managed farming systems deliver mostly provisioning services (e.g., food and fiber), while low-intensity farming systems can support a wider range of ES and high levels of biodiversity (Foley et al. 2005, Swinton et al. 2007, Rockström et al. 2017. Known as High Nature Value farmlands (HNVf), such low-intensity agricultural landscapes are characterized by the prevalence of high levels of natural and seminatural habitats, e.g., heathlands and grasslands (HNVf type 1, hereafter HNVf1), and/or by the occurrence of mosaics of crop fields intermingled with small-scale linear elements, such as field margins and hedgerows (HNVf type 2, hereafter HNVf2) (Andersen et al. 2003, Lomba et al. 2014, Mäkeläinen et al. 2019. To better manage agricultural landscapes, it is therefore essential to understand the relationship between contrasting farmland practices, resulting landscape patterns and biodiversity, and ES delivery. Ecosystem services associations can be identified when changes in one service alter the provision of another service, or when two or more services respond to the same driver of change (Bennett et al. 2009, Spake et al. 2017. ES bundles correspond to the mapping exercise of ES associations and allow investigation of how different ecosystem services occur together across space and time (Raudsepp-Hearne et al. 2010). Hence, ES bundles delineate spatial units delivering similar magnitude and types of ES (Spake et al. 2017) and allow the exploration of synergies and trade-offs between services in distinct social-ecological contexts (Turner et al. 2014, Spake et al. 2017. ES synergies occur when ES are fostered by other ES, while trade-offs occur when the provision of one service is reduced through increased delivery of another (Spake et al. 2017, Sylla et al. 2020. Synergies and trade-offs are typically inferred by positive and negative spatial overlaps, respectively (Spake et al. 2017, Sylla et al. 2020. Ultimately, understanding the link between ES associations and supporting landscapes can provide stakeholders with essential information for the management of complex landscapes . However, understanding how ES and ES bundles relate to, and vary with, specific management practices in human-shaped landscapes needs further exploration (Emmerson et al. 2016).
Here, we seek to fill this knowledge gap by investigating the relationship between ES delivery and social-ecological drivers. We used the Entre-Douro-e-Minho (EDM) region in northern Portugal, a heterogeneous area characterized by the occurrence Fig. 1. The study area, the Entre-Douro-e-Minho (EDM) agrarian region, in the European (a), and regional (b) contexts. The location of the dairy farming region (DFR), Natura 2000 areas, High Nature Value farmlands type 1 and 2 (HNVF1, HNVf2, respectively), and civil parishes are also presented.
The EDM is an agrarian region characterized by heterogeneous landscapes, with larger farms under intensive agricultural practices located primarily in the valleys and on gentle slopes where suitable soils for agriculture are more frequent, and smaller and scattered low-intensity farms prevailing in the mountains (Lomba et al. 2013). The agricultural landscapes consist of farmed areas intermingled with heathlands and sparse vegetation often used for grazing and forest patches (Lomba et al. 2013). A dairy farming region occurs in the coastal area of the EDM associated with the river basins of Cávado, Este, and Ave (Brito et al. 2011, Lomba et al. 2013. Wine production is also prominent in the region and is especially associated with the river basins. Overall, ~ 67% of the EDM is designated as mountain/hill Less Favoured Area for agriculture (Directive 75/268/EEC) and 23% of the area is within EU Natura 2000 Network areas. HNVf occur mainly in the central-east and mountainous areas of the EDM both within and outside Natura 2000 areas (Lomba et al. 2020a).

Data collection and management
Selection of ES indicators and drivers is a key step for the evaluation and understanding of ecosystem services delivery (Spake et al. 2017). Thus, we targeted ES commonly associated with agricultural landscapes and relevant for the study area. The targeted ES are listed in Table 1, for which we provided the rationale underlying their selection (for detailed information see Table A1.1, Appendix 1). Overall, we followed the ES classification defined by CICES (Common International Classification of Ecosystem Services, Haines-Young and Potschin 2018). A total of eight indicators reflecting the supply of different ES were selected, three provisioning, three regulation and maintenance (hereafter regulation), one cultural, and one reflecting support to biodiversity. Although the latter is not defined as an ES under the CICES classification (Haines-Young and Potschin 2018), we considered support to biodiversity conservation an important potential benefit provided by HNV farmlands and related to landscape multifunctionality (Kremen and Merenlender 2018). We then selected proximate drivers of ES supply related to agricultural landscapes and their management practices (Table 1, for detailed information see Table A1.1). The selected drivers reflected farming intensity, specialization, and landscape composition, and provided information on both the intensity and complexity of agricultural landscapes (Persson et al. 2010, Silva et al. 2020, Ribeiro et al. 2021, Santos et al. 2021. All data used in this research were obtained from publicly available databases (Table 1), with the spatial resolution ranging from the parish level to 1x1 km, for the period between 2000 and 2010.
We performed analyses at the scale of the civil parish since this represents the administrative unit at which most data were available. Data available at higher resolutions (e.g., 1 km 2 ) were upscaled, by calculating the median value of each indicator at the parish level. Whenever required, data were divided by the area of each parish, to allow comparisons across parishes with different sizes (Castro et al. 2015). Subsequently, all variables were scaled (z-score method), to allow direct comparisons between variables with different original units (Spake et al. 2017).
We used ArcMap 10.3 and R (R Core Team 2013) to implement all statistical and spatial analyses. Collinearity between drivers was analyzed through the Spearman correlation and the Variance Inflation Factor (VIF), using the R package faraway (Faraway  (Borcard et al. 2011, Dormann et al. 2013; for detailed information see Table A1.2).

Data analysis
To explore ES associations and their relationships with drivers, we performed a multivariate analysis of ES supply and drivers across parishes. This allowed the identification of ES synergies and trade-offs (referring to positive and negative associations of ES, respectively). Then we looked at how ES associated in space within the EDM, by performing bundles analyses delineating groups of parishes supplying similar levels of ES, and identified the most relevant drivers associated with each ES bundle. Finally, to assess the relationships between ES bundles, landscape multifunctionality, and HNVf, we analyzed the variation of multifunctionality and the percentage of HNVf area across ES bundles.

ES supply and drivers
A Principal Component Analysis (PCA) was applied to the selected ES indicators to explore how ES associate in the EDM region by identifying the main factors explaining the variability and distribution of the ES across parishes. Then we performed a Redundancy Analysis (RDA) to analyze the relationship between the social-ecological drivers considered and ES. Specifically, we implemented a stepwise procedure within the RDA analysis to assess drivers' importance, and only those showing statistical significance were retained. The RDA model was evaluated for significance through a permutation test. Both PCA and RDA were conducted using the vegan package for R (Oksanen et al. 2018).

ES bundles and drivers
To analyze ES and drivers at the landscape scale (here represented by groups of parishes), we first applied a clustering procedure using scores from significant PCA axes and derived clusters of parishes, based on their associated ES (Mouchet et al. 2014).
Using PCA axes has the advantage that they represent uncorrelated components of ES data (Legendre and Legendre 2012). We followed the Kaiser-Guttman criterion (eigenvalue > 1) to select significant PCA axes (Legendre and Legendre 2012).
To determine the most appropriate clustering method and optimal number of clusters, we tested different clustering algorithms, including: agglomerative clustering ("hierarchical" and "agnes"); partitioning clustering ("k-means", "pam" and "clara") divisive clustering ("diana" and self-organising tree algorithm -"sota"); and normal mixture model-based clustering ("model"); where two to six clusters were tested (Brock et al. 2008). We used the Euclidian distance to calculate the distance matrix and Ward as the agglomeration method. We then selected the best combination of clustering algorithms and number of clusters based on best values of cluster validation measures and interpretability of bundles of ES for the region. Cluster validation included analysis of internal (silhouette coefficient) and stability (figure of merit -FOM) measures appropriate to determine the best method and the optimal number of clusters for a given dataset (Kassambara 2017). We implemented cluster analyses using the Cluster (Maechler et al. 2018) and clValid (Brock et al. 2008) R packages. We mapped the selected clusters to depict the spatially explicit patterns of bundles of ES. To identify the most important drivers of the spatial patterns observed, we ran a multinomial logit model on ES bundles, using ES drivers as predictors of bundles. We performed the multinomial logit model in R using the mlogit package (Croissant 2013). https://www.ecologyandsociety.org/vol27/iss1/art5/

ES bundles, multifunctionality and HNVf landscapes
To relate bundles with the ES multifunctionality, we estimated landscape multifunctionality as the diversity of ES associations provided for each parish using the transformation (H) of the Gini-Simpson's index (Renard et al. 2015, Spake et al. 2017. This was based on minimum-maximum normalization of ES values. Such transformation provides a measure of the "effective number of services", facilitating comparisons across studies, independently of the diversity index used (Jost 2006, Renard et al. 2015. The values range from 1 to infinity, with higher values depicting higher multifunctionality (Jost, 2006). To summarize multifunctionality for each bundle, we calculated the average and standard deviation of the diversity index for bundles constituting parishes. Finally, to relate bundles with HNVf presence, we calculated the percentage of HNVf1 and 2 occurring within each individual ES bundle. The percentage of HNVf1 and HNVf2 in each parish were available for the agrarian region of Entre-Douro-e-Minho from Lomba et al. (2020a).

Fig. 2.
Principal Component Analysis of the ecosystem services analysed across the study-area. Ecosystem services are presented in red, whereas civil parishes are represented as grey circles.
RDA analysis of ES and drivers produced two gradients explaining most of the data variability (cf. Figure 3 and Table  A2.2, Appendix 2), with RDA1 explaining 12.66 % and RDA2 4.69 % of the variance. RDA1 reflected a gradient of low to high number of farmers (loadings: 0.93) and other drivers related to landscape composition, such as edge density (0.50), land cover diversity (0.39), but also production value (0.39). Overall, parishes with high number of farmers and more complex landscapes are related to provisioning services, such as food crops and wine (2.36, 1.87, respectively) and regulation services (carbon sequestration: 1.25, erosion prevention: 1.04). RDA2 depicted a gradient of farming systems from high to low production value (-0.89), specialization index (-0.21), and associated mainly with cattle (-1.91). A gradient of land cover diversity (0.44) from low to high was also depicted by RDA2, associated with delivery of biodiversity (farmland birds: 0.59) and regulation services (pollination: 0.54).

Assessing bundles of Ecosystem Services and associated socialecological drivers in the EDM
To assess bundles of ES, we used the scores from the first three PCA components in the cluster analysis. Overall, the best performance (reflected as the best compromise between validation measures and interpretability) was achieved through k-means with four clusters (Silhouette: 0.30, FOM: 1.18, see Table  A3.1 in Appendix 3 for detailed results). Below, we present a summary of mean ES values for each bundle, and relate each bundle to drivers and HNVf occurrence, following increasing values of landscape multifunctionality (for full results see Table  A3.2 and A3.4, Appendix 3). Bundle 4 (466 parishes; 37.26% of the EDM) corresponds to a coastal strip, most lowlands, and some mountains in the interior area of EDM (Figure 1 and 4). It includes most urban areas. Bundle 4 delivered the lowest values for most ES, particularly regulation, provisioning, and biodiversity related services, and represented the second lowest delivery of carbon sequestration, cattle, and nature tourism. Number of farmers, diversity of land cover and farm size, significantly decreased the log-odds of bundle 4 occurrence. Bundle 4 presented the lowest value for the multifunctionality index (H:3.97 ± 0.68), a low percentage of HNVf1 (6.45%), and the lowest percentage of HNVf2 (1.18%).
Bundle 3 (381 parishes; 39.40% of EDM), occurred mostly in the lowlands located in the central area of EDM and in Natura 2000 sites located in mountains areas (Figures 1 and 4). Bundle 3 showed the highest delivery of support to biodiversity, second highest delivery of cultural services, and the highest contribution to pollination ( Figure 4). Diversity of land cover significantly increases the log-odds of Bundle 3 occurrence, while production value, edge density, and the number of farmers significantly decrease those log-odds. Bundle 3 had the third highest multifunctionality (H: 4.37 ±0.68), the highest percentage of HNVf1 (17.33%) and a high percentage of HNVf2 (2.65%).
Bundle 2 (212 parishes;10.46% of EDM) mostly overlapped the dairy farming region (Figures 1 and 4). Bundle 2 showed the highest contribution for production of cattle and regulation services, such as erosion prevention and carbon sequestration ( Figure 4). Production value and to a lesser extent the number of farmers significantly increased the log-odds of the Bundle 2 being present, while diversity of land cover significantly decreased the log-odds of Bundle 2 occurrence. The multifunctionality value was the second highest (H: 5.02 ±0.54). Bundle 2 had the lowest percentage of HNVf1 (2.29%) and low percentage of HNVf2 (1.35%).
Bundle 1 (270 parishes; 12.88% of the EDM), occurred in the central area of the EDM region ( Figure 4). Bundle 1 presented the highest values of provisioning (food crops and wine), cultural (nature tourism), and regulation services (pollination) and intermediate values of other regulation services (erosion prevention and carbon sequestration) and support to biodiversity ( Figure 4). The number of farmers, edge density, and to a lesser extent the diversity of land cover and farm size were associated with the occurrence of Bundle 1. Bundle 1 presented the highest value of multifunctionality (H: 5.40± 0.52) and percentage of HNVf2 (5.77%) while presenting a lower percentage of HNVf1 (8.62%).

DISCUSSION
Management of farmlands to deliver multiple ES has been suggested as a strategy to foster sustainability of agricultural landscapes (Lescourret et al. 2015). In the EDM, different farming systems have produced agricultural landscapes contributing differently to the supply of targeted ES in the region. We analyzed ES associations and related ES to farm management drivers to scrutinize patterns of ES delivery in EDM. Overall, we found potential synergies between multiple ES within complex landscapes, and trade-offs between provisioning and other services within intensively managed agricultural landscapes.
Our results showed positive associations between multiple ES, specifically farmland birds, wine, nature tourism and pollination Ecology and Society 27(1): 5 https://www.ecologyandsociety.org/vol27/iss1/art5/ services (cf. Figure 2 and 4). The relationship between pollination and higher biodiversity and rural aesthetics has been linked previously to the conservation of natural vegetation within farming systems (Wratten et al. 2012). Vineyards in the EDM have also been related to high scenic value and are in themselves rural tourism attractions (Marques 2006, Thiele et al. 2019. Such an association was illustrated spatially by Bundle 1, with the highest multifunctionality, and in part by Bundle 3. Bundle 1 reflected potential synergies between support to biodiversity (farmland birds), production (e.g. wine), cultural (nature tourism), and regulation (pollination) ES and occurred scattered across the central region of the EDM and associated with river basins (cf. Figure 1 and 4). It overlapped with HNVf2 areas, supporting previous research highlighting that HNVf landscapes also deliver multiple ES (Plieninger et al. 2019, Lomba et al. 2014, 2020b. Bundle 3 occupied a large extent of the EDM region including mountainous areas and Natura 2000 areas. Here the supply of regulation services (pollination) and biodiversity support (farmland birds) had higher values and related to the occurrence of HNVf1 landscapes. HNVf1 are associated with high amounts of natural and seminatural vegetation supporting the occurrence of many species and habitats of conservation concern (Doxa et al. 2010, Lomba et al. 2020b).
Overall, the observed association of ecosystem services highlighted a spatial trade-off with both provisioning services (cattle production) and regulation services (carbon sequestration and erosion prevention) (depicted spatially by Bundle 2) (cf. Figure 2 and 4). Such trade-offs between provisioning ES and cultural and regulation ES and support to biodiversity have been reported previously (e.g., see Maes et al. 2012, Castro et al. 2015, Lee and Lautenbach 2016, Quintas-Soriano et al. 2016, Sylla et al. 2020. For example, Sylla et al. (2020) reported trade-offs between food provisioning and nutrient cycling and aesthetic features in a peri-urban area of Poland. Castro et al. (2015) and Quintas-Soriano et al. (2016) reported that intensive horticulture production did not co-occur with other ES supplies, such as water regulation or habitat quality, in a study performed in greenhouses in Spain. Bundle 2, characterized by high levels of cattle production, matched the dairy farming region of EDM, an area of high economic importance (Brito et al. 2011). As in other landscapes under intensive farming, Bundle 2 was observed in lowlands associated with river basins where the soils most suitable to farming occur (Raudsepp-Hearne et al. 2010, Turner et al. 2014. Interestingly, Bundle 2, along with Bundle 1, delivered some of the highest values for regulation services, namely carbon sequestration and erosion prevention. This seems to be related to topography combined with the occurrence of tree plantations dominated by eucalyptus and pine. Typically, higher erosion potential is more associated with higher slopes, such as mountains, while tree plantations have been associated with better soil erosion control in the EDM (cf, Figure 1, , Carvalho-Santos et al. 2016, Alegria et al. 2020. For carbon sequestration, higher levels may relate to the high potential for biomass accumulation within tree plantations (Verkerk et al. 2015(Verkerk et al. , 2019. Importantly, the dairy farming region and tree plantations have been associated with negative environmental impacts such as biodiversity loss, pollution of soils, water, and air and streamflow reduction in the EDM (Brito et al. 2011, Carvalho-Santos et al. 2016). Contrastingly, Bundle 4 presented the lowest supply of ES and was found to match areas dominated by urban settlements. While urban areas have been reported to deliver ES, particularly cultural services (Ernstson et al. 2008, Andersson et al. 2014, Queiroz et al. 2015, our analysis focused on ES delivered by agricultural landscapes, which may explain the results obtained.
Our research allowed a better understanding of the socialecological drivers underlying ES delivery. We found that parishes where farmland birds, wine, food crops, nature tourism, and pollination co-occurred were characterized by higher number of farmers, and higher edge density and diversity of land uses, but negatively related to farm size, specialization index and production value (cf. Figure 3). While previous studies have reported both the existence and lack of a relationship between landscape complexity and multifunctionality, our results (Bundle 1) highlight the potential role of more complex and less intensive landscapes, such as HNVf2, in the multifunctional supply of ES and support to biodiversity (Bignal and McCracken 1996, Finney and Kaye 2017, Birkhofer et al. 2018. Bundle 3, characterized by more diverse landscapes, low production values and low number of farmers, showed intermediate delivery of ES. Specifically, Bundle 3 delivery of lower values of provisioning services and higher values of support to biodiversity and pollination, might relate to the presence of HNVf1. HNVf1 occurred mostly in mountain areas and less favored areas for agriculture (Directive 75/268/EEC) which due to their natural constraints and difficulty of access, are less prone to agricultural intensification. Such characteristics may explain higher levels of support to biodiversity as the intensification of landscapes has been linked to the loss of biodiversity, particularly farmland birds (Donald et al. 2001, Reif and Vermouzek 2019). Simultaneously, our results are consistent with other reports of lower supply of ES in mountainous areas, and importantly, such areas have been linked to farmland abandonment (Beilin et al. 2014, Quintas-Soriano et al. 2019, Lomba et al. 2020a. Conversely, results from the RDA showed drivers such as production value, specialization index, and farm size being associated with the provisioning service, cattle. This is also substantiated by the results of drivers associated with Bundle 2, indicating a homogenized landscaped with an industrialized focused production. Finally, socialecological drivers associated with Bundle 4 show that these areas may have the lowest presence of agricultural areas with low number of farmers and low average farm size contributing to low supply of analyzed ES. Overall, our results allow an increased understanding of spatially explicit ES associations across EDM agricultural landscapes. Such results can be used to identify farmlands under similar socioeconomic and ecological conditions and thus to identify priority areas for targeted management actions toward the delivery of ES (Verhagen et al. 2018). Our results represent the first spatial assessment for the EDM that is easy to share and to visualize. This information can be used to inform land managers, support (participatory) decision-making and planning (e.g., at the subregional scale) and, ultimately, support the design of tailored policy instruments at several scales of decision-making (Lescourret et al. 2015, Verhagen et al. 2018, Sylla et al. 2020. Moreover, aligning ES supply and demands at the local and regional scales, and linking results with the policy and socioeconomic factors that underly the management of such https://www.ecologyandsociety.org/vol27/iss1/art5/ landscapes, would represent the next step toward implementing ES-based management in EDM landscapes (Laca et al. 2021). For example, production landscapes that may profit from pollinators were not within the same bundle where pollination services were the highest and fostering the introduction and/or maintenance of hedgerows could improve the delivery of this ES also in more intensively managed farmlands . Importantly, as in most agricultural landscapes in Europe, the EDM farming systems face two contrasting processes: intensification, such as in the case of the dairy farming systems, and abandonment, in the case of HNVf landscapes (Lomba et al. 2020a). The Common Agricultural Policy has failed to support farmers to promote sustainable landscapes and hence to promote the multifunctional potential of rural landscapes (Pe'er et al. 2019, Lomba et al. 2020b). If promoting multifunctional landscapes is considered a desirable outcome for agriculture landscapes, multiple policy instruments should be used, including land management plans, river basin management, forest management, and Common Agricultural Policy (CAP) tools. The latter can be particularly important for promoting heterogeneous landscapes (e.g., keeping field margins) or providing farmers with support to maintain more multifunctional farming systems and prevent land abandonment (Santos et al. 2021).
Our study was performed using publicly available data at the parish level, so that replicability could be assured. This represents a cost-effective approach providing a basis for comparative analysis with other regions, with other time periods, or under scenarios of change (Raudsepp-Hearne et al. 2010, Queiroz et al. 2015. However, the selection of ES and drivers was limited by available data and thus results are limited by the ES selected and scale of analysis (Spake et al. 2017). Within such a complex region as the EDM, these limitations contributed to the lack of variance explained by our PCA and RDA, yet this is not unusual (Turner et al. 2014, Quintas-Soriano et al. 2016, Spake et al. 2017. In fact, analysis of ES bundles at the parish level is limited by documented scale effects, where ES relationships will principally reflect land use distribution (Spake et al. 2017). For example, we can seek to understand the supply at the landscape level but not the contribution of each specific land use. While it was not within the scope of this study, and in any case data was not available at the time of this research, a multiscale analysis would allow an investigation of cross-scale effects of ES associations (Scholes et al. 2013, Spake et al. 2017. Thus, while our results inform on pattern-based multifunctionality they do not inform on processbased multifunctionality which would further improve the understanding of ES relationships (Spake et al. 2017). However, we selected ES on which data were available and proximate drivers considered relevant to describe the agricultural landscapes in the EDM region (Plieninger et al. 2016, Santos et al. 2020. Therefore, our findings represent an essential first assessment of the main ES associations in the EDM.

CONCLUSIONS
Our research aimed to achieve a better understanding of ES and their relationships within agricultural landscapes. This work identified regions that provide a more multifunctional supply of ES, e.g., complex and less intensified HNVf landscapes, and the potential ES synergies and trade-offs occurring at the regional level. It also demonstrated the spatial trade-offs occurring within this region, where the former multifunctional landscapes trade-off with highly specialized farming system landscapes producing provisioning services. Being dynamic entities, agricultural landscapes in the study region are continuously changing with some farming systems specializing, such as the dairy sector, and others decreasing in representation (due to e.g., abandonment), such HNVf landscapes. Such farming trajectories will impact the potential delivery of ES of these landscapes in the future.
The dairy sector in the EDM is of high economic importance for the region. But the analyses presented here show how that supply is associated with much lower delivery of other ES such as nature tourism and farmland birds. These dairy areas do occur in association with some other ES, specifically carbon sequestration and erosion prevention, though these are more likely to reflect historic land use features in the wider landscape (in the form of tree planting) and not the characteristics of the dairy farms themselves. Conversely, a greater range of additional ES is associated with the complex landscapes where HNVf primarily occurs, and in most cases, those ES will be linked to the characteristics of the farms themselves. However, the limited income achievable from the agricultural products from these HNVf together with the current lack of any payment/reward system for these additional ES means that the continued economic viability of those farms is low, leading to abandonment and the loss of those ES.
The knowledge produced here highlights the conflicts that need to be considered when seeking to inform the planning and management of complex social-ecological systems. The results also highlight the importance of low-intensity farming systems, such as occur on High Nature Value farmlands, as they provide a wider range of ES delivery within those landscapes.
Responses to this article can be read online at: https://www.ecologyandsociety.org/issues/responses. php/12947

Appendix 1 | ES and drivers' indicators.
We based the selection of ES and drivers on their relevance for the agricultural landscapes of the Entre Douro Minho (EDM) region. In appendix 1, we present the ES and drivers ( Figure A1.1), their description and rationale (Table A1.1) and analyses conducted to evaluate correlation and collinearity within services and drivers (Table A1.2, A1.3, A1.4). We also present the results from the calculation of ES multifunctionality for the parishes of the EDM (Figure A1.2).    Table A1.3 -Results from the correlation analysis of pairs of ecosystem services with the Spearman correlation test. Significance at p < 0.05 denoted by * and at p < 0.01 by **.  Table A1.4 -Results from the correlation analysis of pairs of drivers with the Spearman correlation test.

Appendix 2| ES supply and drivers
To explore ES associations and analyse their potential drivers, we perform a PCA and RDA of ES and ES and drivers, respectively. In appendix 2, we present the PCA and RDA results. For the PCA, we present the eigenvalues and variance explained of the eigenvalues above 1 following Kaiser-Guttman criterion to select significant PCA axes (Legendre and Legendre, 2012), ES loadings (Table A 2.1) and the mapping of scores for the significant PCA axes in the EDM ( Figure   A 2.1). For the RDA, we present the variance explained by the significant constrained axes, and ES and drivers' scores (

ES bundles
To understand how ES and their drivers spatially associate at the landscape scale, we performed a clustering analysis with PCA scores of significant PCA axes. For this we tested different clustering algorithms and numbers of cluster and we presented validation measures results (Table A3.1). For the selected cluster results we present the mean values of each ES for each bundle (Table A3.2) and the mean values of the PC scores for each bundle (Table A3.3).