Analysing the Synergies and Trade-Offs between Ecosystem Services to Reorient Land Use Planning in Metropolitan Bilbao (Northern Spain)

In the last decades, some European cities have undergone important changes in search of a more sustainable development. This is the case for the city of Bilbao (Bizkaia, Basque Country), where a Greenbelt has been maintained surrounding the urban areas allowing the periurban areas to deliver ecosystem services (ES) to society. However, the role of the different ecosystems in the provision of ES is not the same, which can lead to conflicts among them. The aim of this study is to analyze the synergies and trade-offs among the eight most important ES in the Bilbao Metropolitan Greenbelt (BMG) to orient their management strategies towards more multifunctional landscapes. We mapped the ES and overlapped them looking for the most relevant areas for the provision of multiple ES and areas that are mostly lacking ES provision. We identify also existing ES trade-offs and synergies between ES using correlations so that managers can prioritize preservation efforts of land use types in the rest of the area. The results show that provisioning ES had trade-offs with regulating and cultural ES and the latter showed synergies between them. The former are mainly delivered by semi-natural ecosystems, while regulating and cultural ES are delivered mainly by natural ecosystems. Moreover, the most relevant areas for the provision of multiple ES were proposed as potential components of a Green Infrastructure (GI). Their identification and ES bundles could help decision-makers to orient their management strategies towards sustainability in metropolitan areas.


Introduction
Multifunctional landscapes have the capacity to simultaneously produce multiple ecosystem services (ES) [1,2]. However, human activities have modified many landscapes throughout the world to supply desired ES, such as food, timber, and fiber, maximizing their production, but resulting in substantial declines in the provision of other ES [3]. Trade-offs in ES arise when management choices entail the optimization of a single or a few ES, leading to the reduction or deterioration of other ES [2] this situation being inevitable when managing multiple ES [4]. However, ES synergies have been defined as 'a situation where the use of one ES directly increases the benefits supplied by another service' [5]. Thus, it is necessary to understand the relationships between ES to improve our ability to sustainably manage landscapes to provide multiple ES [6,7]. Moreover, understanding the synergies and trade-offs resulting from different management objectives may allow stakeholders and the public to conduct informed discussions about future management policies [8,9].
Especially when trade-offs occur, and different interests must be weighed against each other, it is necessary that decisions, e.g., on the enhancement of certain ES, be spatially explicit [10]. The spatial visualization approach involving the use of mapping tools constitutes a powerful tool that can contribute to environmental and landscape decision-making [11][12][13][14], as it helps planners and managers understand the complexity and associated trade-offs and synergies among ES [6] or identify areas of conflicting goals [15]. Moreover, the provision of ES changes over space and time as a result of changing patterns of land use or changes in the composition and structure of different ecosystems [16,17]. Therefore, it is important to assess the spatial heterogeneity of the landscape to provide such services [18], at scales that are relevant for managers and policymakers [19].
Most analyses of ES trade-offs have focused on intact natural vegetation or protected areas [20][21][22][23]. However, in the last years more attention has been paid to the assessment of ES in more complex social-ecological systems such as metropolitan and urban areas [24][25][26][27][28], as the need to protect ES within human-dominated landscapes has been recognized by practitioners [29][30][31]. Urban populations typically obtain most of their ecosystem resources from sources that are distributed over a substantially larger area [32]. Nevertheless, if it is planned and managed appropriately, these areas can support a wide variety of ES, such as air purification, recreation, food production and aesthetic value, playing a key role in the local provision of ES to their occupants and enhancing the well-being of billions of people [33,34]. Therefore, effective management of ES provision within metropolitan areas is an important issue. Thus, analyses of ES synergies and trade-offs, as well as, the identification of multifunctional key areas, e.g., components of a GI, enable policymakers to maximize ES from a preferred policy option and they also facilitate decisions concerning optimized land use in metropolitan areas. Green Infrastructure is a network of natural, semi-natural areas and green spaces that delivers essential ES, which underpin human well-being and quality of life. Maintaining ES through the development of GI is therefore increasingly recognized by policies as a strategy to cope with potentially changing conditions in the future [35,36].
Within this context, in the last decades the environmental changes suffered by the Bilbao Metropolitan Greenbelt (BMG), the periurban green space surrounding the urban areas of metropolitan Bilbao, have reduced its capacity to provide ES as the area has lost multifuncionality. However, the BMG is well valued by its residents and visitors for the different ES that it supplies now and its capacity to supply them in the future, as a transformation of the BMG is expected and required by local citizens [37].
The main purpose of this study is to analyze the synergies and trade-offs among the most important ES provided by the BMG to orient management strategies of this area towards more multifunctional landscapes [38]. The analyzed ES were Provisioning ES: food and timber production; Regulating ES: habitat maintenance, air purification, carbon storage, and water flow regulation; and Cultural ES: recreation, and aesthetic quality [39]. To this end, we specifically: (1) Map these eight ES to explore their spatial distribution.
(2) Identify the most relevant areas for the provision of multiple ES and areas that are mostly lacking ES provision, identifying the corresponding land use type. (3) Identify existing ES trade-offs and synergies so that managers can prioritize preservation efforts of land use types in the rest of the study area. (4) Offer recommendations to reorient land use planning practices to maximize ES delivery in the BMG.

Study Area
The social-ecological system of metropolitan Bilbao is located in the region of Bizkaia in the Basque Country, northern Spain ( Figure 1). It is divided into 29 municipalities and covers an area of 413 km 2 with a population of 893,697 inhabitants, resulting in a high population density (2164 inhabitants per km 2 ). In this region, urban areas are situated in the valley along the estuary of the river Nervión-Ibaizabal, which is delimited by small mountains and the coast to the north. The periurban ecosystems surrounding the urban areas of metropolitan Bilbao form the BMG and occupy almost 75% of the metropolitan area. These ecosystems include beaches, cliffs, rivers, meadows, scrublands, natural forests and forests plantations. Mixed-oak forests, Cantabrian evergreen-oak forests, beech forests and riparian forests represent the potential vegetation of the study area. However, these natural forests occupy only 12.3% of the BMG, while 36.9% of the BMG is occupied by forest plantations, as a result of the reforestation performed during the 20th century. Grasslands and crops cover 27.4% of the area and mainly produce peppers, apples, local white wine (txakoli), milk, cheese and livestock ( Figure 1). periurban ecosystems surrounding the urban areas of metropolitan Bilbao form the BMG and occupy almost 75% of the metropolitan area. These ecosystems include beaches, cliffs, rivers, meadows, scrublands, natural forests and forests plantations. Mixed-oak forests, Cantabrian evergreen-oak forests, beech forests and riparian forests represent the potential vegetation of the study area. However, these natural forests occupy only 12.3% of the BMG, while 36.9% of the BMG is occupied by forest plantations, as a result of the reforestation performed during the 20th century. Grasslands and crops cover 27.4% of the area and mainly produce peppers, apples, local white wine (txakoli), milk, cheese and livestock ( Figure 1).

Selection of Ecosystem Services and Land Use Units
Residents and visitors of the BMG indicated that recreation (71.8%), air purification (26.4%), aesthetic value (9.6%), habitat maintenance (9.0%), food and timber production (1.6%), carbon storage (0.8%), and water flow regulation (0.4%), in that order, were the most important benefits supplied by the BMG [37]. Therefore, these eight ES were selected on the basis of their importance in the area, their relevance to conservation planning of BMG and the availability of data.

Selection of Ecosystem Services and Land Use Units
Residents and visitors of the BMG indicated that recreation (71.8%), air purification (26.4%), aesthetic value (9.6%), habitat maintenance (9.0%), food and timber production (1.6%), carbon storage (0.8%), and water flow regulation (0.4%), in that order, were the most important benefits supplied by the BMG [37]. Therefore, these eight ES were selected on the basis of their importance in the area, their relevance to conservation planning of BMG and the availability of data.
The land use types used to map the ES were defined according to the European Nature Information System (Habitat EUNIS) map, scale 1:10,000 of the Basque Country for year 2009 [40] as it is the only information available for this territory. The 108 habitats found in the study area were aggregated into the 16 land use types most relevant to the area ( Figure 1).

Mapping and Evaluating Selected ES
Although many studies have derived information on ES directly from land use/cover or habitat maps [15,41,42], we used spatially explicit biophysical indicators to analyze the spatial distribution of ES delivery (Table A1). Moreover, once each ES indicator was calculated, they were normalized to a common scale due to their different units of analysis. Therefore, they were scaled and mapped into five ranges (1 = very low or null; 2 = low; 3 = medium; 4 = high; and 5 = very high), related to the capacity of different areas of the BMG to provide each ES. These ranges were defined using Jenks Natural Breaks [23,36,43], except for provisioning services. For geoprocessing we used the software ArcGIS 10.3 [44].

Food Production
Two components were considered in mapping the food production service [36]: agricultural productivity (tonnes ha −1 ) and livestock productivity (tonnes ha −1 ) (Table A1). Thus, the areas considered for providing this service were crops (orchards, vineyard and fruit trees) for agricultural production, and grasslands for livestock production.
The productions and surfaces of the different crops and livestock were obtained from the Statistic Institute of the Basque Country [45]. Data were calculated for all the area of Biscay (Table A2). Based on these data each land use type was assigned the following values: very high or 5 score to orchards, and kiwi crops (≥10.75 tonnes ha −1 ); and high or 4 score to fruit trees, vineyards and grasslands (<10.75 tonnes ha −1 ). Grasslands were given high score based on the livestock productivity of the Basque Country [36]. Finally, food productivity (tonnes ha −1 ) was calculated by aggregating agricultural productivity and livestock productivity maps.

Timber Production
This ES was calculated based on the annual growth (m 3 ha −1 year −1 ) of the different forest plantations of the Forest Inventory of the Basque Country obtained by LIDAR flights of 2008 and 2011 [46]. Although natural forests have the capacity to supply timber, these ecosystems are not being productively used. Hence, they were considered to have a low timber production value, and only forest plantations were ranged as: low or 2 score (<10 m 3 ha −1 year −1 ); medium or 3 score (10-15 m 3 ha −1 year −1 ); high or 4 score (15-20 m 3 ha −1 year −1 ); and very high or 5 score (>20 m 3 ha −1 year −1 ) to normalize the indicator.

Regulation ES Habitat Maintenance
Habitat maintenance is the ecosystem capacity to supply conditions appropriate for the survival and reproduction of a given organism. The indicator used to evaluate this ES was habitat maintenance index, which was based on three components: native plant species richness, habitat quality and sites of natural interest [23,47,48].
On the other hand, habitat maturity was estimated based on prior literature [23,61,62] as mature habitats present the highest habitat suitability for native species of animal, fungi, and cryptogams [63]. According to these criteria, the assigned values for each land use type were the following: 4 = natural forests and other natural ecosystems, 3 = heathlands and shrubs, 2 = grasslands; 1 = forest plantations and crops, and 0 = urban and artificial areas.
Finally, sites of natural interest were considered because their designation is based on the presence of relevant flora, fauna, and singular habitat. They include those protected areas already designated as part of the Natura 2000 network and by the Nature Conservation Law 16/1994 of the Basque Country (Biotopes). They also include areas that have the habitat listed in Annex I of the European Habitat Directive [64] but that are not protected, as the importance of the habitat is not dependent on the protection status (political decision). The sites of natural interest were valued with 1 and those without a defined natural interest with 0. We have to recognize limitations of this assumption, as many areas that have not been included could contain noticeable environmental values. Finally, habitat maintenance index was calculated by aggregating the values of the three components described above [65] (Table A1).
Subsequently, data of this indicator were ranged as follows: values of 2-3 = 1 or very low or null; values of 4-5 = 2 or low; values of 6-8 = 3 or medium; values of 9-10 = 4 or high; and values of 11-12 = 5 or very high).

Air Purification
The air purification ES focuses on the air pollutant NO 2 , as its abatement is still a pressing challenge in most major urban and periurban areas in Spain [25]. This ES was assessed using the indicator of annual NO 2 removal (g m −2 h −1 ) calculated as the product of NO 2 concentration (g m −3 ) and dry deposition velocity (m s −1 ) [25] (Table A1). The method for the calculation of NO 2 concentration in BMR was based on a combination of a linear regression model and kriging interpolation (geostatistical method) based on the measurements of the maximum daily value of NO 2 (in summer only) of the air quality stations [66] and on the NO 2 dry deposition velocity. The latter was calculated based on literature [67][68][69][70] (Table A3).
Subsequently, data of this indicator were ranged as follows: values of <0.

Carbon Storage
Following [23], we estimated the amount of carbon stored in biomass and soil for the different land use types. For valuation of the carbon stored in soil (tonnes C ha −1 ), we used the "Inventory of organic carbon stored in the first 30 cm of the soil" of the Basque Country [71]. As for carbon stored as biomass, we considered that, in ecosystems other than forests, the amount of carbon stored as biomass was insignificant compared with the carbon stored in forests or in the soil. For example, Spanish shrublands only represent 8.2% of the tree carbon stock in Spanish forests [72] the rest is stored in living trees. Therefore, in this study, we focused on the carbon stored in living trees (aboveground and belowground) calculated multiplying the merchantable volume (m 3 /ha) [46], the biomass expansion factor for the conversion of merchantable volume to aboveground tree biomass to include branches and leaves (dimensionless variable) [73], the root-to-shoot ratio to include belowground tree biomass (dimensionless variable) [73], the basic wood density (tonnes dry matter/m 3 merchantable volume) [74,75] and the carbon fraction of dry matter (tonnes C/tonnes dry matter) [73,76] (Table A3). Finally, carbon storage was calculated by aggregating the values of the carbon stored in living biomass and soil (Table A1).
Subsequently, data of this indicator were ranged as follows: values of <4 tonnes C ha −1 = 1 or very low or null; values of 4-97 tonnes C ha −1 = 2 or low; values of 98-160 tonnes C ha −1 = 3 or medium; values of 161-214 tonnes C ha −1 = 4 or high; and values of >214 tonnes C ha −1 = 5 or very high).

Water Flow Regulation
This ES involves the influence of natural systems on the regulation of surface water flows, and it is a function of the storage and retention components of the water flow [77]. We used the fraction of annual water flow stored in the soil to measure the water flow regulation service [23]. Its calculation was based on the hydrological simulation model TETIS developed for the region [78]. The volume of water produced in the area was determined primarily by the rainfall patterns, which depend mainly on abiotic parameters (regional climate and topography). Moreover, ecosystems play a key role in water flow due to the amount of water they retain in the soil and return to the atmosphere via evapotranspiration. Therefore, three components were used to calculate the water flow regulation index: (1) water storage in the soil (mm year −1 ); (2) annual rainfall (mm year −1 ); and (3) corrected annual potential evapotranspiration (mm year −1 ) [23,79] (see in Table A1 the equation used).
The potential evapotranspiration for the different vegetation types was modified using correction factors based on those used by InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs) [47] to obtain a more realistic value for evapotranspiration [23]. Calculation was made for each hydrological unit within the study area (Ibaizabal, Butroe and Barbadun). The water bodies were given the maximum value for water flow regulation.

Recreation
We considered two components for mapping the recreation index: the recreation potential and the recreation opportunity [36,80,81]. The recreation potential is defined as the capacity of ecosystems to provide recreation according to their scenic beauty or specific characteristics, and the recreation opportunity is defined as the capacity of ecosystems to provide recreation according to their accessibility and the presence of infrastructures specific for recreation in their environment [82]. The ES delivery strictly depends on the presence of people in the ecosystem, thus it is necessary for them to get there in order to benefit from ES therein. Moreover, the presence of infrastructures specific to recreation in their environment attracts people to visit the area.
Previous studies have determined that protected areas and ecosystems with high naturalness are more attractive for recreation activities due to their high biodiversity [13,83]. The presence of water bodies, sites of geological interest and mountain summits imply a higher recreational value as well, as people can perform different recreational activities in them (bathing, canoeing, trekking, geological tourism, etc.). Moreover, water body type is also relevant when considering the number of people that can attract. For example, beaches attract much more people than a river in the Basque Country [66]. Therefore, the recreation potential was mapped taking into account five territorial features associated with aesthetic and ecological attractiveness for recreational activities: (1) degree of naturalness; (2) presence of protected areas or sites of naturalistic interest; (3) presence of water bodies; (4) presence of sites of geological interest; and (5) presence of mountain summits. We assessed each feature using specific data of the study area and indicators (see Table A4).
The degree of naturalness is defined as degree of human influence on ecosystems and it comprises two aspects, the transformations caused by humans in ecosystems and how these ecosystems depend on human activity themselves [84]. To estimate the degree of naturalness of the BMG ecosystems we used the index proposed by [84] for the Iberian Peninsula, accommodating it to the land use types defined in the study area (Table A4). Subsequently, in 2007, the presence of protected areas or sites of naturalistic interest, water bodies, sites of geological interest, and mountain summits were assessed with a value of 1, except for beaches and water bodies used for fishing or bathing, which were assessed with a value of 3 and 2, respectively. As mentioned above, these water body types attack more public than other water bodies [66]. Finally, the recreation potential was calculated by aggregating the values of the five features described above. All components were considered equally important, covering complementary aspects of recreation potential; therefore, they were given equal weights, within and among them [85].
In general, good accessibility and good infrastructure networks favor recreational activities [82,83]. Therefore, the recreation opportunity was mapped considering: (1) accessibility of the sites; and (2) natural and constructed infrastructures that were in place to guide or be enjoyed by visitors. We assessed each feature using different data and indicators (see Table A4). BMG is a very populated area where the accessibility is very high as there is a high density of roads and paths [81]. Therefore, we considered as very accessible sites in the BMG those areas that were in a buffer of 200 m around accessible roads to motor vehicles, and accessible sites those areas that were in a buffer of 200 m around limited roads to motor vehicles where the access was on foot [80]. Moreover, we considered that in a buffer of 500 m around natural and constructed infrastructures visitors could use the area for recreation [86]. The recreation opportunity was calculated aggregating the values of the two features described above, and giving equal weights to both of them.
Finally, recreation index was calculated adding the values of the recreation potential and opportunity.
Subsequently, data of this indicator were ranged as follows: values of <4 =1 or very low or null; values of 4 = 2 or low; values of 5 = 3 or medium; values of 6 = 4 or high; and values of >6 = 5 or very high).

Aesthetic Value
The aesthetic enjoyment of landscapes that the environment provides to society depends on both the aesthetic perception of the society and the intrinsic characteristics of the landscapes. Therefore, the aesthetic enjoyment of landscape was mapped taking into account six features associated with aesthetic attractiveness: (1) the social perception; (2) the type of relief; (3) the diversity of landscape; (4) the presence of water bodies; (5) the influence of landmarks; and (6) the influence of negative elements in the landscapes. We valued each feature using different data and indicators (see Table A5). In general, the presence of water bodies and/or landmarks in landscape [13,83,87], a higher diversity of landscapes, and a higher difference in relief [88,89] were related to a higher aesthetic value [90]. However, the influence of negative elements in the landscapes as wind farms, active quarries, landfills, roads, and railroads were related to a lower aesthetic value [91,92].
The social perception was calculated using a photo-questionnaire online to which 629 persons of the Basque Country completed [81]. Photo-questionnaire consisted of a battery of photos of the different land use types (two for each land use type that appeared randomly), where respondents were asked to rate each photo on a scale from 1 to 6 according to their scenic beauty with 6 being the highest value and 1 the lowest. The mean values obtained for each land use type were used as a proxy for mapping the social perception.
For assessing the type of relief and the diversity of landscape, we used viewshed as a quantification unit. Viewshed can be considered a unit for describing the territory based on visibility criteria, because it consists of the set of intervisible points [93]. Therefore, we considered that viewsheds with index of relief higher or equal to 32 and index of landscape diversity higher or equal to 1.70 in the Basque Country corresponded to areas with high slopes and to areas with a diverse landscape, respectively [93].
Moreover, the presence of landscape landmarks offers an added aesthetic value to the landscape that owns it. However, you can only enjoy them when they are visible from the point where you are. Therefore, according to the Catalog of Singular and Outstanding Landscapes of the Basque Country [93] we considered that the area of influence of these landmarks is all the area located at a distance of less than two kilometers from which these landmarks are visible.
In the case of the influence of negative elements in the landscapes to determine the area of influence or buffer, we relied also on the previous Catalog [93].
Finally, the aesthetic beauty was calculated by aggregating the values of the first five features described above and subtracting the value of the influence of negative elements in the landscapes, giving equal weights to all of them.
Subsequently, data of this indicator were ranged as follows: values of <4 =1 or very low or null; values of 4 = 2 or low; values of 5 = 3 or medium; values of 6 = 4 or high; and values of >6 = 5 or very high).

Data Analysis
All ES maps were overlapped using the Raster Calculator of the software ArcGIS 10.3 [44] to spatially visualize the most relevant areas for the provision of multiple ES and those areas that are mostly lacking ES provision. Those areas with high or very high values (values = 4 or 5) for at least 5 ES (habitat maintenance + 4 ES) were proposed as principal components of the GI. Those areas with high, very high or medium values (values = 5-3) for at least 5 ES (habitat maintenance + 4 ES) were proposed as secondary components of the GI. In BMG the habitat maintenance is considered as a priority ES due to scarcity of natural areas and fragmentation of the territory. Finally, those areas with low or very low or null values (values = 2 or 1) for all ES were considered as areas that are mostly lacking ES provision ( Figure 2). Subsequently, this map was overlapped with the map of land use types to analyze how they contributed to the different identified areas.
In addition, to identify existing ES synergies and trade-offs so that managers can prioritize preservation efforts of land use types in not proposed as components of the GI areas correlations among ES were performed using Spearman's rank correlation [65] except between carbon storage and timber production, as both rely on similar data. Moreover, using the correlation matrix, a principal component analysis (PCA) was applied to distinguish spatial synergies, trade-offs and gradients between ES and land use types. This information would be useful for managers to design the GI, and to improve the multifunctionality in the landscape. Statistical analysis was performed in the R software environment (Version 3.2.2; R Development Core Team). Finally, the aesthetic beauty was calculated by aggregating the values of the first five features described above and subtracting the value of the influence of negative elements in the landscapes, giving equal weights to all of them.
Subsequently, data of this indicator were ranged as follows: values of <4 =1 or very low or null; values of 4 = 2 or low; values of 5 = 3 or medium; values of 6 = 4 or high; and values of >6 = 5 or very high).

Data Analysis
All ES maps were overlapped using the Raster Calculator of the software ArcGIS 10.3 [44] to spatially visualize the most relevant areas for the provision of multiple ES and those areas that are mostly lacking ES provision. Those areas with high or very high values (values = 4 or 5) for at least 5 ES (habitat maintenance + 4 ES) were proposed as principal components of the GI. Those areas with high, very high or medium values (values = 5-3) for at least 5 ES (habitat maintenance + 4 ES) were proposed as secondary components of the GI. In BMG the habitat maintenance is considered as a priority ES due to scarcity of natural areas and fragmentation of the territory. Finally, those areas with low or very low or null values (values = 2 or 1) for all ES were considered as areas that are mostly lacking ES provision ( Figure 2). Subsequently, this map was overlapped with the map of land use types to analyze how they contributed to the different identified areas.
In addition, to identify existing ES synergies and trade-offs so that managers can prioritize preservation efforts of land use types in not proposed as components of the GI areas correlations among ES were performed using Spearman's rank correlation [65] except between carbon storage and timber production, as both rely on similar data. Moreover, using the correlation matrix, a principal component analysis (PCA) was applied to distinguish spatial synergies, trade-offs and gradients between ES and land use types. This information would be useful for managers to design the GI, and to improve the multifunctionality in the landscape. Statistical analysis was performed in the R software environment (Version 3.2.2; R Development Core Team).

Figure 2.
Workflow followed to produce the map where the most relevant areas for the provision of multiple ecosystem services (ES) and areas that are mostly lacking ES provision were identified.

Mapping and Evaluating Selected ES
The maps indicate high spatial variation of the examined services across the BMG. In the case of provisioning services 20% of the BMG showed high or very high values for both ES. On the contrary, regulating services, air purification and carbon storage, were the ES that were best supplied by the study area (38% and 37% of the BMG, respectively, with high or very high values). However, high or very high values for habitat maintenance and water flow regulation were only found in 15% and 11% of the study area, respectively. For cultural services, namely aesthetic value and recreation, 23% and 21% of the BMG, respectively, showed high or very high values (Figure 3).  The most relevant areas of the BMG for the provision of multiple ES were proposed as potential components of the GI in BMG (principal (3%) and secondary components (22%)) ( Figure 4). The land use types that corresponded to principal components of the GI were mainly natural forests (81% mixed-oak forests, 11% riparian forests, and 3% evergreen oak forests), and those that corresponded to secondary components were grasslands (48%), natural forests (21%) and shrubs and heaths (19%). On the other hand, 20% of the BMG do not provide ES or provide them with low or very low values. These areas corresponded mainly to urban areas (89%), mines and quarries (8%) and urban parks (3%).
As it is very difficult to change the land use type in those areas some measures are proposed later in the discussion section that will help increase the ES provision on a more local scale. In the rest of the BMG, managers need information to prioritize some land use types over others depending on the ES that they want to prioritize. Therefore, the following analysis of synergies and trade-offs between the ES can offer them such information.

Analysis of Synergies and Trade-Offs
The significant positive correlations found between ES mean that a high capacity for carbon storage is related with a high capacity for habitat maintenance, air purification and water flow regulation, while a high capacity for timber production is related with a high capacity for water flow regulation; a high capacity for habitat maintenance with high capacity for recreation; and finally, a high capacity for recreation is related with a high aesthetic value (p < 0.0001) ( Table 1). On the contrary, weak negative correlations found between ES means that a high capacity of land use for food production is related with a low capacity for timber production, habitat maintenance, carbon storage, water flow regulation and recreation, while a low aesthetic value is related with a high capacity for timber production and carbon storage (p < 0.05) ( Table 1).

Figure 4.
Spatial distribution of the most relevant areas for the provision of multiple ES, which will be the potential components of the Green Infrastructure (GI) (principal and secondary components of the GI), areas that are mostly lacking ES provision, and percentage of surface they occupy in BMG. Areas with high or very high values (values = 4 or 5) for at least 5 ES (habitat maintenance + 4 ES) were proposed as principal components of the GI. Areas with high, very high or medium values (values = 5-3) for at least 5 ES (habitat maintenance + 4 ES) were proposed as secondary components of the GI. Areas with low or very low or null values (values = 2 or 1) for all ES were considered as areas that are mostly lacking ES provision.

Analysis of Synergies and Trade-Offs
The significant positive correlations found between ES mean that a high capacity for carbon storage is related with a high capacity for habitat maintenance, air purification and water flow regulation, while a high capacity for timber production is related with a high capacity for water flow regulation; a high capacity for habitat maintenance with high capacity for recreation; and finally, a high capacity for recreation is related with a high aesthetic value (p < 0.0001) ( Table 1). On the contrary, weak negative correlations found between ES means that a high capacity of land use for food production is related with a low capacity for timber production, habitat maintenance, carbon storage, water flow regulation and recreation, while a low aesthetic value is related with a high capacity for timber production and carbon storage (p < 0.05) ( Table 1). The relationships between the ES provided by the BMG and the land use types considered are shown in Figure 5, where the first two components of the PCA are presented. Factor 2 (explaining 44.4% of the variance) and Factor 1 (explaining 23.5% of the variance) represent a trade-off between food production and the rest of the ES analyzed and between ES provided mainly by forest and shrubs ecosystems and those supplied by other land use types, respectively (see Appendix A: Table A6). The upper left quadrant represents food production; the lower left quadrant represents ES provided mainly by forest and shrubs ecosystems (timber production and regulation services except for habitat maintenance); and the lower right quadrant represents those ES supplied by natural ecosystems (cultural services and habitat maintenance). Therefore, land use types subjected to greater human influence, such as forest plantations, mines and quarries, orchards and crops, grasslands and urban parks were located from the diagonal upwards, whereas natural ecosystems (coastal habitats, rocky vegetation, aquatic ecosystems, and natural forests) were located from the diagonal downwards. Therefore, if managers want to prioritize, for example, food production ES in some areas, they have to be aware that regulating and cultural ES provision will be reduced. On the other hand, if they prioritize land use types as natural forest or shrubs, multiple ES provision will increase.  The relationships between the ES provided by the BMG and the land use types considered are shown in Figure 5, where the first two components of the PCA are presented. Factor 2 (explaining 44.4% of the variance) and Factor 1 (explaining 23.5% of the variance) represent a trade-off between food production and the rest of the ES analyzed and between ES provided mainly by forest and shrubs ecosystems and those supplied by other land use types, respectively (see Appendix: Table  A6). The upper left quadrant represents food production; the lower left quadrant represents ES provided mainly by forest and shrubs ecosystems (timber production and regulation services except for habitat maintenance); and the lower right quadrant represents those ES supplied by natural ecosystems (cultural services and habitat maintenance). Therefore, land use types subjected to greater human influence, such as forest plantations, mines and quarries, orchards and crops, grasslands and urban parks were located from the diagonal upwards, whereas natural ecosystems (coastal habitats, rocky vegetation, aquatic ecosystems, and natural forests) were located from the diagonal downwards. Therefore, if managers want to prioritize, for example, food production ES in some areas, they have to be aware that regulating and cultural ES provision will be reduced. On the other hand, if they prioritize land use types as natural forest or shrubs, multiple ES provision will increase.

Provision of Multiple ES
In recent years, efforts are being made to transfer the ES framework to land use planning and policy-making activities [94], as it has much potential to help policy-makers and practitioners to identify, protect and prioritize areas for biodiversity conservation in human-dominated landscapes [95,96]. Moreover, the assessment of multiple ES can help to identify on one hand, the most relevant areas for the provision of multiple ES and on the other, areas that are mostly lacking ES provision.
In this study, the spatial distribution of the provision of the ES assessed in the BMG was highly heterogeneous, as might be expected in a district characterized by intensive human activities and a highly fragmented landscape. Nevertheless, some multifunctional areas were identified, which should be considered as potential components of the GI in the BMG. The determination of the most relevant areas for the provision of multiple ES in a periurban area, such as BMG, could help decision-makers to develop a GI, as not all green areas qualify as GI [97]. In BMG, the land use types that fell within the category of principal components of GI were mainly natural forests that should be conserved. The secondary components of GI corresponded to grasslands, some natural forests and shrubs and heaths that should be considered for their potential contribution in the connectivity of the principal components of GI, as connectivity is an essential feature of GI. Maintaining ES as proposed under target 2 of the EU biodiversity strategy requires expanding the network of GI [35], therefore this type of studies are necessary to expand or create the GI in the periurban areas.
On the other hand, areas that are mostly lacking ES provision corresponded to urban areas, mines and quarries and urban parks. It is true that, on an urban scale, urban parks supply important ES [25], and they can be highly significant for the flow of ES [98]. However, at periurban scale less modified ecosystems, such as natural forests, showed the capacity to supply a greater variety and overall higher level of the assessed ES than these highly modified ecosystems. Therefore, in highly modified ecosystems, where land use type cannot be modified, different strategies at local scale must be developed to increase ES provision, such as increasing the number and size of green areas in the urban areas, naturalization of abandoned mines and quarries, or renaturalization of canalized rivers (nature-based solutions).

Trade-Offs, Synergies and Bundles of ES
The identification of synergies, trade-offs and the land use types that contribute to the ES provision allows policymakers to better understand the hidden consequences of preferring one ES to another [99], helping on the prioritization of management strategies [3,100]. Moreover, it helps also to make better-informed decisions to create landscapes that balance conservation with production and promote ecological, social and economic resilience [6,100,101]. In the case of BMG, synergies were mainly found among the regulating ES (habitat maintenance, air purification, and water flow regulation). Synergies were also found between the cultural ES, and between habitat maintenance and recreation. These synergies have also been found in other studies [22,102,103]. All regulation ES could be considered as a bundle of ES, which are supplied mainly by forest ecosystems; therefore, their conservation is really relevant. However, forest plantations (25% coniferous plantation, most in private land) despite their high contribution to some regulating ES should be carefully considered as their management techniques are quite aggressive, producing disservices associated with soil erosion, water quality [104] and regional diversity [105]. In addition, the current conifer timber production is not such a profitable activity in the area due to globalization of the timber market and other factors, and its subsistence depends heavily on public subsidies [106]. Therefore, some of them could be transformed in other land use types, such as natural forests, grasslands or into agroforestry systems to increase the production of multiple ES in the area and to get a more multifunctional landscape [65]. Moreover, more sustainable forest management practices should be applied as Belgium, France or Germany has already done. For example, plantations could include different species and/or different age classes disposed unevenly across the landscape to mimic natural forests while being harvested here and there periodically [107,108]. This allows them to get the better of two worlds: maintaining ES provision and economic returns in this area. This could be made possible by means of incentives which have the ability to motivate changes in land use and management, especially among the most conflicting land use types [17]. Another ES bundle should come from the synergies between habitat maintenance and cultural ES supplied mainly by natural ecosystems (aquatic systems, natural forests, coastal habitats . . . ).
On the other hand, trade-offs were also observed in BMG mainly between aesthetic value and timber production [21,30]. This can be due to the fact that Basque population considers the presence of forest plantations (coniferous and eucalyptus) to be detrimental to the landscape as they cause homogeneity in it. Other trade-offs were also observed between food production and some regulating and cultural ES [42,109]. Although increasing provisioning ES inevitably results in losses of regulating and cultural ES [35], it is clear that they are also relevant for society, not only from the point of view of production but also for their aesthetic value. Some agroecosystems (grasslands, orchards or vineyards) are considered cultural landscapes, being attractive to people [37,81]. Nevertheless, this land use type should be located in the best areas for it [65] and it must be managed sustainably. The identification of these synergies and trade-offs is important for conservation planning [110], especially in areas subject to intensive human activities [111].

Application of the ES Approach to the Management of Metropolitan Areas. Planning for the Delivery of Multiple ES
Concomitant with the increase in the proportion of the urban population, the awareness of the importance of ES delivery to them has grown [32]. Metropolitan areas are of significance in this regard. These locations provide major opportunities to align planning, landscape ecology and ecosystem management to enhance the well-being of billions of people [7]. Therefore, investments in conservation within human-dominated landscapes are needed, such as those made in wilderness areas, to provide both key ES and biodiversity benefits [30] and thus, reducing the regional and global footprint of cities and towns [32,38,112].
Regional land management has favored the delivery of some ES in metropolitan areas while reducing the capacity of ecosystems to provide relevant ES. However, this is not necessarily irreversible, as changing the actual land use configuration would result in changes in the ES provision [3,16]. If the goal is to maximize the provision of ES, it should be considered substituting, reducing or diversifying the areas occupied by the least multifunctional land use types in favor of more multifunctional landscapes. A useful approach in this regard could be to implement "land sharing" to strengthen the delivery of multiple ES, instead of dedicating large areas to obtaining only provisioning services (land sparing) [113,114]. Nevertheless, prior to deciding which ES we want to favor or maximize and determining the changes that will be necessary for that purpose, it is essential to take into account the interests and demands of the society, as a specific location will need to be managed differently considering alternative situations [115][116][117]. Nevertheless, some measures should be implemented in relation to land use management if social demands are to be met [37]. In this regard, natural forests should be considered high-priority areas because, currently, these forests are fragmented and sparse. Thus, management interventions should be focused on maintaining these areas, connecting the existing patches and increasing their size. An increase in rural activities could potentially be implemented in areas currently abandoned improving the aesthetic value of the area, as rural landscapes have a higher aesthetic value than those where such activities have been abandoned [37,81].

Conclusions
In conclusion, the location of the multifunctional areas based on ES is a good method to highlight the land use types providers of the ES needed or demanded in the area and those land use types not providing them. Thus, this information would be necessary to find the best sites to be integrated into the GI and the potentiality of other sites to be integrated in the GI through management interventions. This would help to manage the area prioritizing actions to preserve the ES provider sites and, also to design plans to convert the low ES provision sites on better providers. A limitation of this assessment is that we assume that all ES will be valued the same, however, unfortunately, land management decisions are mostly based on economic trade-offs rather than ES provision ones, which totally could change the final picture.
However, the design of these plans needs to take into account the interests and needs of the society, not forgetting the land ownership and, thus promoting the transformation of more multifunctional sites through positive incentives from local/regional governments. The design of this landscape based on the ES framework should help to improve the well-being for the inhabitants of the area as well as fulfilling their demands. Acknowledgments: We deeply appreciate the suggestions of the anonymous referees, which greatly improved the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.
Appendix A Table A1. Indicators for mapping ES [39] and the variables and equations used to calculate them.

Ecosystem Service Indicators Variables Used to Calculate Indicators Equations
Provisioning Food production Food productivity (F) (tonnes ha −1 ) A p = Agricultural productivity (tonnes ha −1 )