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, Power 2010, Halada et al. 2011, 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, Power 2010, 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 (Meacham et al. 2016). 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 of a gradient of contrasting farming systems (Carvalho 2011, Lomba et al. 2020a), as our study area. Our research aimed specifically to: 1) investigate ES spatial associations across the EDM and identify bundles of ES delivery; and, 2) contribute to a better understanding of the relation of social-ecological drivers with ES associations. We discuss our results in the context of landscape multifunctionality, the occurrence of HNVf in the study area, and in the broader context of planning and managing agricultural landscapes.

METHODS

Study Area

The EDM is located in northwest Portugal (Figure 1) and consists of 53 municipalities and 1326 civil parishes (local administrative units, Eurostat, http://epp.eurostat.ec.europa.eu/). The EDM includes a mixture of lowland areas, large valleys, and mountain massifs, with altitudes ranging from sea level to ca. 1500m. The topography closely follows the basins of the main rivers, Minho, Lima, Cávado, Ave, and Douro, all with an east-northeast-west-southwest orientation.

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²) 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 2011). Indicators were considered for further analyses when exhibiting correlation values under 0.7 and VIF values under 10, i.e. when found to be below the commonly defined thresholds (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).

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).

RESULTS

Patterns of ecosystem services delivery and their drivers

PCA analysis (Figure 2, for detailed information see Table A2.1, Appendix 2) produced three principal components (PCs) explaining 53% of the variance (eigenvalues > 1). PC1 (23.43% of total variance) represented a gradient of parishes with increasing potential for the delivery of four main ES: wine (loading: 0.50), food crops (0.48), erosion prevention (0.43), and carbon sequestration (0.41). PC2 (16.75%) contrasted parishes with high delivery of support to farmland birds, pollination, and nature tourism (loadings: 0.48; 0.36; 0.30, respectively) with parishes with higher supply of provisioning and regulation ES (loadings: cattle, -0.49; carbon sequestration, -0.42; erosion prevention, -0.30). PC3 (12.95%) represented a gradient of parishes with increasing potential for the delivery of food crops, wine, and nature tourism (loadings: 0.29; 0.22; 0.15, respectively) and decreasing potential of delivery of farmland birds and cattle (loadings: -0.69; -0.44, respectively).

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 social-ecological 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 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, Halada et al. 2011, 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, Perez-Soba et al. 2010, 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, 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 social-ecological 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, Power 2010, 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, social-ecological 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 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 (Kremen et al. 2018). 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 process-based 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.

LITERATURE CITED

Alegria, C., N. Roque, T. Albuquerque, S. Gerassis, P. Fernandez, and M. M. Ribeiro. 2020. Species ecological envelopes under climate change scenarios: a case study for the main two wood-production forest species in Portugal. Forests 11:880. https://doi.org/10.3390/f11080880

Andersson, E., S. Barthel, S. Borgström, J. Colding, T. Elmqvist, C. Folke, and Å. Gren. 2014. Reconnecting Cities to the Biosphere: Stewardship of Green Infrastructure and Urban Ecosystem Services. AMBIO 43:445-453. https://doi.org/10.1007/s13280-014-0506-y

Beilin, R., R. Lindborg, M. Stenseke, H. M. Pereira, A. Llausàs, E. Slätmo, Y. Cerqueira, L. Navarro, P. Rodrigues, N. Reichelt, N. Munro, and C. Queiroz. 2014. Analysing how drivers of agricultural land abandonment affect biodiversity and cultural landscapes using case studies from Scandinavia, Iberia and Oceania. Land Use Policy 36:60-72. https://doi.org/10.1016/j.landusepol.2013.07.003

Bennett, E. M., G. D. Peterson, and L. J. Gordon. 2009. Understanding relationships among multiple ecosystem services. Ecology Letters 12:1394-1404. https://doi.org/10.1111/j.1461-0248.2009.01387.x

Bignal, E. M., and D. I. Mccracken. 1996. Low-Intensity Farming Systems in the Conservation of the Countryside. Journal of Applied Ecology 33:413-424. https://doi.org/10.2307/2404973

Birkhofer, K., G. K. S. Andersson, J. Bengtsson, R. Bommarco, J. Dänhardt, B. Ekbom, J. Ekroos, T. Hahn, K. Hedlund, A. M. Jönsson, R. Lindborg, O. Olsson, R. Rader, A. Rusch, M. Stjernman, A. Williams, and H. G. Smith. 2018. Relationships between multiple biodiversity components and ecosystem services along a landscape complexity gradient. Biological Conservation 218:247-253. https://doi.org/10.1016/j.biocon.2017.12.027

Borcard, D., F. Gillet, and P. Legendre. 2011. Numerical ecology with R. Use R! series, Springer Science, New York. https://doi.org/https://doi.org/10.1007/978-1-4419-7976-6

Brito, L., J. Alonso, J. Mamede, and J. Rey-Graña. 2011. Slurry management in the milk farms of Entre Douro e Minho. Revista de Ciências Agrárias (Portugal) 34:80-93.

Brock, G., V. Pihur, S. Datta, and S. Datta. 2008. clValid: An R Package for Cluster Validation. 2008 25:22. https://doi.org/https://doi.org/10.18637/jss.v025.i04

Carvalho, A. 2011. Recenseamento Agrícola. 2009: Análise dos Principiais Resultados. INE-Instituto Nacional de Estatística. Lisboa, Portugal.

Carvalho-Santos, C., R. Sousa-Silva, J. Gonçalves, and J. P. Honrado. 2016. Ecosystem services and biodiversity conservation under forestation scenarios: options to improve management in the Vez watershed, NW Portugal. Regional Environmental Change 16:1557-1570. https://doi.org/10.1007/s10113-015-0892-0

Castro, A. J., B. Martín-López, E. López, T. Plieninger, D. Alcaraz-Segura, C. C. Vaughn , and J. Cabello. 2015. Do protected areas networks ensure the supply of ecosystem services? Spatial patterns of two nature reserve systems in semi-arid Spain. Applied Geography 60:1-9. https://doi.org/10.1016/j.apgeog.2015.02.012

Croissant, Y. 2013. mlogit: multinomial logit model. R package version 0.2-4. R Foundation for Statistical Computing: Vienna, Austria.

Declerck, F. A. J., S. K. Jones, S. Attwood, D. Bossio, E. Girvetz, B. Chaplin-Kramer, E. Enfors, A. K. Fremier, L. J. Gordon, F. Kizito, I. Lopez Noriega, N. Matthews, M. Mccartney, M. Meacham, A. Noble, M. Quintero, R. Remans, R. Soppe, L. Willemen, S. L. R. Wood, and W. Zhang. 2016. Agricultural ecosystems and their services: the vanguard of sustainability? Current Opinion in Environmental Sustainability 23:92-99. https://doi.org/10.1016/j.cosust.2016.11.016

DGT 2007. Carta de Uso e Ocupação do Solo de Portugal Continental para 2007 - (COS2007v1.0). Retrieved December 20, 2016. [online] URL: http://mapas.dgterritorio.pt/.

Donald, P. F., R. E. Green, and M. F. Heath. 2001. Agricultural intensification and the collapse of Europe's farmland bird populations. Proceedings of the Royal Society of London. Series B: Biological Sciences 268:25-29. https://doi.org/10.1098/rspb.2000.1325

Dormann, C. F., J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carré, J. R. G. Marquéz, B. Gruber, B. Lafourcade, P. J. Leitão, T. Münkemüller, C. Mcclean, P. E. Osborne, B. Reineking, B. Schröder, A. K. Skidmore, D. Zurell , and S. Lautenbach. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36:27-46. https://doi.org/10.1111/j.1600-0587.2012.07348.x

Doxa, A., Y. Bas, M. L. Paracchini, P. Pointereau, J. M. Terres, and F. Jiguet. 2010. Low-intensity agriculture increases farmland bird abundances in France. Journal of Applied Ecology 47:1348-1356. https://doi.org/10.1111/j.1365-2664.2010.01869.x

Emmerson, M., M. B. Morales, J. J. Oñate, P. Batáry, F. Berendse, J. Liira, T. Aavik, I. Guerrero, R. Bommarco, S. Eggers, T. Pärt, T. Tscharntke, W. Weisser, L. Clement, and J. Bengtsson. 2016. Chapter Two - How Agricultural Intensification Affects Biodiversity and Ecosystem Services. In: Dumbrell, A. J., Kordas, R. L., and Woodward, G. (eds.) Advances in Ecological Research. Academic Press. https://doi.org/10.1016/bs.aecr.2016.08.005

Ernstson, H., S. Sörlin , and T. Elmqvist. 2008. Social Movements and Ecosystem Services—the Role of Social Network Structure in Protecting and Managing Urban Green Areas in Stockholm. Ecology and Society 13(2):39. https://doi.org/https://doi.org/10.5751/ES-02589-130239

Faraway, J. 2011. faraway: Functions and datasets for books by Julian Faraway. R package version, 1, 54.

Finney, D. M., and J. P. Kaye. 2017. Functional diversity in cover crop polycultures increases multifunctionality of an agricultural system. Journal of Applied Ecology 54:509-517. https://doi.org/10.1111/1365-2664.12765

Foley, J. A., R. Defries, G. P. Asner, C. Barford, G. Bonan, S. R. Carpenter, F. S. Chapin, M. T. Coe, G. C. Daily, H. K. Gibbs, J. H. Helkowski, T. Holloway, E. A. Howard, C. J. Kucharik, C. Monfreda, J. A. Patz, I. C. Prentice, N. Ramankutty , and P. K. Snyder. 2005. Global Consequences of Land Use. Science 309:570-574. https://doi.org/10.1126/science.1111772

Haines-Young, R., and M. B. Potschin. 2018. Common international classification of ecosystem services (CICES) V5. 1 and guidance on the application of the revised structure. European Environment Agency, 53. https://doi.org/https://doi.org/10.3897/oneeco.3.e27108

Halada, L., D. Evans, C. Romão, and J.-E. Petersen. 2011. Which habitats of European importance depend on agricultural practices? Biodiversity and Conservation 20:2365-2378. https://doi.org/10.1007/s10531-011-9989-z

INE 2009. Recenseamento agrícola - séries históricas. Localização geográfica NUTS-2001. Retrieved July 25, 2013 [online] URL: http://ra09.ine.pt

Jost, L. 2006. Entropy and diversity. Oikos 113:363-375. https://doi.org/10.1111/j.2006.0030-1299.14714.x

Kassambara, A. 2017. Practical guide to cluster analysis in R: Unsupervised machine learning, Sthda.

Kremen, C., and A. M. Merenlender. 2018. Landscapes that work for biodiversity and people. Science 362:eaau6020. https://doi.org/10.1126/science.aau6020

Kremen, C., L. K. M'gonigle, and L. C. Ponisio. 2018. Pollinator Community Assembly Tracks Changes in Floral Resources as Restored Hedgerows Mature in Agricultural Landscapes. Frontiers in Ecology and Evolution 6. https://doi.org/10.3389/fevo.2018.00170

Lee, H., and S. Lautenbach. 2016. A quantitative review of relationships between ecosystem services. Ecological Indicators 66:340-351. https://doi.org/10.1016/j.ecolind.2016.02.004

Legendre, P., and L. Legendre. 2012. Numerical ecology, Elsevier Scientific.

Lescourret, F., D. Magda, G. Richard, A.-F. Adam-Blondon, M. Bardy, J. Baudry, I. Doussan, B. Dumont, F. Lefèvre, I. Litrico, R. Martin-Clouaire, B. Montuelle, S. Pellerin, M. Plantegenest, E. Tancoigne, A. Thomas, H. Guyomard, and J.-F. Soussana. 2015. A social-ecological approach to managing multiple agro-ecosystem services. Current Opinion in Environmental Sustainability 14:68-75. https://doi.org/10.1016/j.cosust.2015.04.001

Lomba, A., A. Buchadas, E. Corbelle-Rico, R. Jongman, and D. Mccracken. 2020a. Detecting temporal changes in the extent of High Nature Value farmlands: The case-study of the Entre-Douro-e-Minho Region, Portugal. Landscape and Urban Planning 195:103726. https://doi.org/https://doi.org/10.1016/j.landurbplan.2019.103726

Lomba, A., A. S. Vaz, F. Moreira, and J. Pradinho Honrado. 2013. Hierarchic species-area relationships and the management of forest habitat islands in intensive farmland. Forest Ecology and Management 291:190-198. https://doi.org/10.1016/j.foreco.2012.11.012

Lomba, A., C. Guerra, J. Alonso, J. P. Honrado, R. Jongman, and D. Mccracken. 2014. Mapping and monitoring High Nature Value farmlands: Challenges in European landscapes. Journal of Environmental Management 143:140-150. https://doi.org/10.1016/j.jenvman.2014.04.029

Lomba, A., F. Moreira, S. Klimek, R. H. Jongman, C. Sullivan, J. Moran, X. Poux, J. P. Honrado, T. Pinto-Correia, T. Plieninger, and D. I. Mccracken. 2020b. Back to the future: rethinking socioecological systems underlying high nature value farmlands. Frontiers in Ecology and the Environment 18:36-42. https://doi.org/10.1002/fee.2116

Maechler, M., P. Rousseeuw, A. Struyf, M. Hubert, K. Hornik, M. Studer, and P. Roudier. 2018. Package CLUSTER-Finding groups in data: Cluster analysis. CRAN R studio.

Maes, J., M. L. Paracchini, G. Zulian, M. B. Dunbar, and R. Alkemade. 2012. Synergies and trade-offs between ecosystem service supply, biodiversity, and habitat conservation status in Europe. Biological Conservation: 155:1-12. https://doi.org/10.1016/j.biocon.2012.06.016

Mäkeläinen, S., A. Harlio, R. K. Heikkinen, I. Herzon, M. Kuussaari, K. Lepikkö, A. Maier, T. Seimola, J. Tiainen, and A. Arponen. 2019. Coincidence of High Nature Value farmlands with bird and butterfly diversity. Agriculture, Ecosystems, and Environment 269:224-233. https://doi.org/10.1016/j.agee.2018.09.030

Marques, H. 2006. Research Report: Searching for Complementarities between Agriculture and Tourism — The Demarcated Wine-Producing Regions of Northern Portugal. Tourism Economics 12:147-155. https://doi.org/10.5367/000000006776387141

Meacham, M., C. Queiroz, A. V. Norström, and G. D. Peterson. 2016. Social-ecological drivers of multiple ecosystem services: what variables explain patterns of ecosystem services across the Norrström drainage basin? Ecology and Society 21(1):14. https://doi.org/10.5751/ES-08077-210114

Mouchet, M. A., P. Lamarque, B. Martín-López, E. Crouzat, P. Gos, C. Byczek, and S. Lavorel. 2014. An interdisciplinary methodological guide for quantifying associations between ecosystem services. Global Environmental Change 28:298-308. https://doi.org/10.1016/j.gloenvcha.2014.07.012

Oksanen, F. J., M. Friendly, R. Kindt, P. Legendre, D. Mcglinn, P. R. Minchin, R. O’hara, G. L. Simpson, P. Solymos , and M. H. H. Stevens. 2018. Vegan: community ecology package. R package version, 2. 5-2. [online] URL: https://CRAN.R-project.org/package=vegan

Pe'er, G., Y. Zinngrebe, F. Moreira, C. Sirami, S. Schindler, R. Müller, V. Bontzorlos, D. Clough, P. Bezák, A. Bonn, B. Hansjürgens, A. Lomba, S. Möckel, G. Passoni, C. Schleyer, J. Schmidt, and S. Lakner. 2019. A greener path for the EU Common Agricultural Policy. Science 365:449-451. https://doi.org/10.1126/science.aax3146

Perez-Soba, M., P. Verburg, E. Koomen, M. Hilferink, P. Benito, J. Lesschen, and M. Banse. 2010. Land use modelling-implementation; Preserving and enhancing the environmental benefits of land-user services. Geodan Next, Wageningen.

Persson, A. S., O. Olsson, M. Rundlöf, and H. G. Smith. 2010. Land use intensity and landscape complexity—Analysis of landscape characteristics in an agricultural region in Southern Sweden. Agriculture, Ecosystems, and Environment 136:169-176. https://doi.org/10.1016/j.agee.2009.12.018

Plieninger, T., H. Draux, N. Fagerholm, C. Bieling, M. Bürgi, T. Kizos, T. Kuemmerle, J. Primdahl, and P. H. Verburg. 2016. The driving forces of landscape change in Europe: A systematic review of the evidence. Land Use Policy 57:204-214. https://doi.org/10.1016/j.landusepol.2016.04.040

Plieninger, T., M. Torralba, T. Hartel, and N. Fagerholm. 2019. Perceived ecosystem services synergies, trade-offs, and bundles in European high nature value farming landscapes. Landscape Ecology 34:1565-1581. https://doi.org/10.1007/s10980-019-00775-1

Power, A. G. 2010. Ecosystem services and agriculture: tradeoffs and synergies. Philosophical Transactions of the Royal Society B: Biological Sciences 365:2959-2971. https://doi.org/10.1098/rstb.2010.0143

Queiroz, C., M. Meacham, K. Richter, A. V. Norström, E. Andersson, J. Norberg, and G. Peterson. 2015. Mapping bundles of ecosystem services reveals distinct types of multifunctionality within a Swedish landscape. AMBIO 44:89-101. https://doi.org/10.1007/s13280-014-0601-0

Quintas-Soriano, C., A. J. Castro, H. Castro, and M. García-Llorente. 2016. Impacts of land use change on ecosystem services and implications for human well-being in Spanish drylands. Land Use Policy 54:534-548. https://doi.org/10.1016/j.landusepol.2016.03.011

Quintas-Soriano, C., M. García-Llorente, A. Norström, M. Meacham, G. Peterson, and A. J. Castro. 2019. Integrating supply and demand in ecosystem service bundles characterization across Mediterranean transformed landscapes. Landscape Ecology 34:1619-1633. https://doi.org/10.1007/s10980-019-00826-7

R Core Team. 2013. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. [online] URL: http://www.R-project.org/.

Raudsepp-Hearne, C., G. D. Peterson, and E. M. Bennett. 2010. Ecosystem service bundles for analyzing tradeoffs in diverse landscapes. Proceedings of the National Academy of Sciences 107:5242-5247. https://doi.org/10.1073/pnas.0907284107

Reif, J., and Z. Vermouzek. 2019. Collapse of farmland bird populations in an Eastern European country following its EU accession. Conservation Letters 12:e12585. https://doi.org/10.1111/conl.12585

Renard, D., J. M. Rhemtulla, and E. M. Bennett. 2015. Historical dynamics in ecosystem service bundles. Proceedings of the National Academy of Sciences 112:13411-13416. 10.1073/pnas.1502565112 https://doi.org/10.1073/pnas.1502565112

Ribeiro, P. P. D. C. F., J. M. Osório De Barros De Lima E Santos, M. J. Prudêncio Rafael Canadas, A. M. Contente De Vinha Novais, F. M. Ribeiro Ferraria Moreira, and Â. C. De Araújo Rodrigues Lomba. 2021. Explaining farming systems spatial patterns: A farm-level choice model based on socioeconomic and biophysical drivers. Agricultural Systems 191:103140. https://doi.org/10.1016/j.agsy.2021.103140

Rockström, J., J. Williams, G. Daily, A. Noble, N. Matthews, L. Gordon, H. Wetterstrand, F. Declerck, M. Shah, P. Steduto, C. De Fraiture, N. Hatibu, O. Unver, J. Bird, L. Sibanda, and J. Smith. 2017. Sustainable intensification of agriculture for human prosperity and global sustainability. Ambio 46:4-17. https://doi.org/https://doi.org/10.1007/s13280-016-0793-6

Santos, J. L., F. Moreira, P. F. Ribeiro, M. J. Canadas, A. Novais, and A. Lomba. 2021. A farming systems approach to linking agricultural policies with biodiversity and ecosystem services. Frontiers in Ecology and the Environment 19:168-175. https://doi.org/10.1002/fee.2292

Scholes, R. J., B. Reyers, R. Biggs, M. J. Spierenburg, and A. Duriappah. 2013. Multi-scale and cross-scale assessments of social-ecological systems and their ecosystem services. Current Opinion in Environmental Sustainability 5:16-25. https://doi.org/10.1016/j.cosust.2013.01.004

Schulp, C. J. E., G.-J. Nabuurs, and P. H. Verburg. 2008. Future carbon sequestration in Europe—Effects of land use change. Agriculture, Ecosystems, and Environment 127:251-264. https://doi.org/10.1016/j.agee.2008.04.010

Schulp, C. J. E., S. Lautenbach, and P. H. Verburg. 2014. Quantifying and mapping ecosystem services: Demand and supply of pollination in the European Union. Ecological Indicators 36:131-141. https://doi.org/10.1016/j.ecolind.2013.07.014

Silva, J. F., J. L. Santos, P. F. Ribeiro, M. J. Canadas, A. Novais, A. Lomba, M. R. Magalhães, and F. Moreira. 2020. Identifying and explaining the farming system composition of agricultural landscapes: The role of socioeconomic drivers under strong biophysical gradients. Landscape and Urban Planning 202:103879. https://doi.org/10.1016/j.landurbplan.2020.103879

Spake, R., R. Lasseur, E. Crouzat, J. M. Bullock, S. Lavorel, K. E. Parks, M. Schaafsma, E. M. Bennett, J. Maes, M. Mulligan, M. Mouchet, G. D. Peterson, C. J. E. Schulp, W. Thuiller, M. G. Turner, P. H. Verburg, and F. Eigenbrod. 2017. Unpacking ecosystem service bundles: Towards predictive mapping of synergies and trade-offs between ecosystem services. Global Environmental Change 47:37-50. https://doi.org/10.1016/j.gloenvcha.2017.08.004

Swinton, S., F. Lupi, G. P. Robertson, and S. K. Hamilton. 2007. Ecosystem services and agriculture: Cultivating agricultural ecosystems for diverse benefits. Ecological Economics 64:245-252. https://doi.org/10.1016/j.ecolecon.2007.09.020

Sylla, M., N. Hagemann, and S. Szewrański. 2020. Mapping trade-offs and synergies among peri-urban ecosystem services to address spatial policy. Environmental Science, and Policy 112:79-90. https://doi.org/10.1016/j.envsci.2020.06.002

Thiele, J., C. Von Haaren, and C. Albert. 2019. Are river landscapes outstanding in providing cultural ecosystem services? An indicator-based exploration in Germany. Ecological Indicators 101:31-40. https://doi.org/10.1016/j.ecolind.2019.01.003

Tucker, G., B. Allen, M. Conway, I. Dickie, K. Hart, M. Rayment, C. Schulp, and A. Van Teeffelen. 2013. Policy options for an EU no net loss initiative. Report to the European Commission. Institute for European Environmental Policy, London.

Turner, K. G., M. V. Odgaard, P. K. Bøcher, T. Dalgaard, and J. C. Svenning. 2014. Bundling ecosystem services in Denmark: Trade-offs and synergies in a cultural landscape. Landscape and Urban Planning 125:89-104. https://doi.org/10.1016/j.landurbplan.2014.02.007

Van Berkel, D. B., and P. H. Verburg. 2011. Sensitising rural policy: Assessing spatial variation in rural development options for Europe. Land Use Policy 28:447-459. https://doi.org/10.1016/j.landusepol.2010.09.002

Verhagen, W., E. H. Van Der Zanden, M. Strauch, A. J. A. Van Teeffelen, and P. H. Verburg. 2018. Optimizing the allocation of agri-environment measures to navigate the trade-offs between ecosystem services, biodiversity and agricultural production. Environmental Science, and Policy 84:186-196. https://doi.org/10.1016/j.envsci.2018.03.013

Verkerk, P. J., C. Levers, T. Kuemmerle, M. Lindner, R. Valbuena, P. H. Verburg, and S. Zudin. 2015. Mapping wood production in European forests. Forest Ecology and Management 357:228-238. https://doi.org/10.1016/j.foreco.2015.08.007

Verkerk, P. J., J. B. Fitzgerald, P. Datta, M. Dees, G. M. Hengeveld, M. Lindner, and S. Zudin. 2019. Spatial distribution of the potential forest biomass availability in Europe. Forest Ecosystems 6:5. https://doi.org/10.1186/s40663-019-0163-5

Wood, S. L. R., S. K. Jones, J. A. Johnson, K. A. Brauman, R. Chaplin-Kramer, A. Fremier, E. Girvetz, L. J. Gordon, C. V. Kappel, L. Mandle, M. Mulligan, P. O'farrell, W. K. Smith, L. Willemen, W. Zhang, and F. A. Declerck. 2018. Distilling the role of ecosystem services in the Sustainable Development Goals. Ecosystem Services 29:70-82. https://doi.org/10.1016/j.ecoser.2017.10.010

Wratten, S. D., M. Gillespie, A. Decourtye, E. Mader, and N. Desneux. 2012. Pollinator habitat enhancement: Benefits to other ecosystem services. Agriculture, Ecosystems, and Environment 159:112-122. https://doi.org/10.1016/j.agee.2012.06.020