Mapping and Assessment of Ecosystems Services under the Proposed MAES European Common Framework: Methodological Challenges and Opportunities

The EU Biodiversity Strategy for 2020 was a driving force behind spatially explicit quantifications of Ecosystem Services (ES) in Europe. In Portugal, the MAES initiative (ptMAES–Mapping and Assessment of Ecosystem and their Services) was conducted in 2014 to address Target 2 (Action 5) of the Strategy, namely mapping and assessing ecosystems, ecosystems’ condition (EC), and ES. In this study covering the NUTS II Alentejo region, EC was assessed and mapped based on four indicators (soil organic matter, plant and bird diversity, and ecological value of plant communities) and five ES were assessed and mapped (soil protection, carbon sequestration, and fiber/crop/livestock production). Assessments were performed under a multi-tiered approach, ranging from spatialization of statistical data to analytical modeling, based on the most detailed land-use/land-cover cartography available. In this paper, we detail the methodological and analytical framework applied in ptMAES and present its main outcomes. Our goal is to (1) discuss the main methodological challenges encountered to inform future MAES initiatives in Portugal and other member states; and (2) further explore the outcomes of ptMAES by looking into spatial relationships between EC and ES supply. We highlight the advantages of the proposed analytical framework and identify constraints that, among others, limited the number of ES and EC indicators analyzed. We also show that MAES can provide useful insights to landscape planning at the regional scale, for instance, red-flagging areas where ES supply may be unsustainable over time.


Introduction
Ecosystem services (ES) are the benefits humans derive from ecosystems [1], which are categorized into provisioning services (e.g., food production), regulating services (e.g., control of erosion rates), and cultural services (e.g., aesthetic enjoyment of nature) by the Common International Classification of Ecosystem Services (CICES) [2], a mainstream classification within the EU. Although the concept of ES is not exempt from criticism [3], it brings into perspective that humans (and human actions by extension) largely depend on well-functioning ecosystems, emphasizing the need for developing sustainable and resilient interactions between human societies and the natural environment [4]. This imposes new challenges when it comes to landscape planning, calling for an understanding of how human actions across different landscapes can best safeguard well-functioning ecosystems and maximize the services they generate [5]. A core concept in understanding and addressing the sustainability of ES is ecosystem condition [6]. Here, we refer to the ecological condition of ecosystems (EC) as the state of ecological systems, which includes their physical, chemical, and biological characteristics and the processes and interactions that connect them [7]. Biodiversity is commonly understood as one possible indicator of the ecological condition of ecosystems [8,9].
Spatially explicit quantifications of ES (i.e., mapping of ES) gained particular prominence in Europe under the EU Biodiversity Strategy for 2020 [10] 1 . Spatial-based approaches facilitate decision making by providing an efficient way of conveying complex information through visual representation [11] and are valuable in systematic conservation planning to ensure the long-term capacity of ecosystems to provide ES [3]. Despite the increasing number of ES mapping and spatial assessment exercises in the literature, only a small amount have explored Action 5, Target 2 of the EU Biodiversity Strategy for 2020 2 from an integrated scientific perspective, namely by pointing out the relationship between the ecological condition of ecosystems and their service supply capacity under a long-term sustainability perspective [12]. A major challenge lies in the complexity of the relationships between EC and ES [13] and in the fact that these relationships greatly vary depending on the scale of analysis [14,15], which may not coincide with the scale of landscape planning.
In Portugal, the MAES initiative (ptMAES-Mapping and Assessment of Ecosystem and their Services) was set out to address Target 2 of the EU Biodiversity Strategy by taking the challenges imposed by its actions 5 and 7a), namely, "Member States, with the assistance of the Commission, will map and assess the state of ecosystems and their services in their national territory by 2014, assess the economic value of such services, and promote the integration of these values into accounting and reporting systems at EU and national level by 2020". The ptMAES initiative was carried out in 2014 as a regional assessment targeting the NUTS II Alentejo region, which covers about one-third of the total area of mainland Portugal. The goal of ptMAES was to map the dominant ecosystems within the NUTS II Alentejo region 3 following the European Nature Information System (EUNIS) habitat classification and to map and assess ecosystem's condition and services provided by the agricultural, forest, and agroforestry ecosystems in the region. The study comprised a multi-tiered approach to spatially quantify EC and ES, following the standardized methodological approaches proposed in the MAES guidelines to support the use of ES in planning [8,16,17]. With this approach, a set of accurate and defendable methods were selected, intended to deliver outcomes that could better assist decision making aimed at landscape planning under nature conservation goals [18]. Although it presented a pioneering exercise in this regard, the study fell short of providing further clues on the spatial relationships between EC and ES, which may potentially assist the actual integration of ES into landscape planning in the region.
In this paper, we detail the methodological and analytical framework developed in ptMAES to map and assess EC and ES supply, explore spatial relationships between EC and ES supply, and discuss the main challenges and opportunities encountered in order to inform the implementation of future MAES initiatives in Portugal and other member states. As such, the goal of this paper is two-fold. We aimed at developing a methodological reference for national policies furthering the implementation of the EU Biodiversity Strategy to 2020 (namely those concerned with Target 2, Action 5) and of the key actions set out in the EU Biodiversity Strategy to 2030. These national policies include, among others, the Portuguese Nature Conservation and Biodiversity Strategy 2015-2020 and 2020-2030, the National Strategy for Forests up to 2030, and the Commitment for Green Growth 2020/2030. We also aimed at exploring the potentialities of the proposed methodology by further looking into spatial relationships between EC and ES.

Study Area
The ptMAES initiative was carried out as a regional study in Alentejo (Figure 1), which presents a typically Mediterranean climate regime [19], dominated by rural landscapes with high potential for tourism [20]. The study focused on the rural landscapes of Alentejo, and as such, only forest, agroforestry, agricultural, and shrubland land use/land cover (LULC) were targeted in the analysis (other LULC were classified as "no data", Figure 1).

Mapping Ecosystems' Condition and Service Supply
Following MAES guidelines, the methodological framework chosen to map and assess ecosystems' condition and services consists of a multi-tiered approach, combining methods of varying complexity given data availability. The decision on which EC indicators and ES to assess was based on data availability to test the study under the limited time frame usually associated with the implementation period set out in EU strategies/regulations. As such, the selection here presented was supported by discussions with representatives of the National Authority for Nature Conservation (Instituto da Conservação da Natureza e Florestas-ICNF) to accommodate official requirements while accounting for data availability and approved methods to be used under the proposed time frame for the development of this first MAES initiative in Portugal. Figure 2 presents a schematic representation of our approach, and each method is detailed in the following subsections. All assessments were based on the most detailed LULC cartography publicly available in Portugal, from 2007 (hereinafter COS07), which has a minimum mapping unit of 25 ha. Therefore, the unit of mapping and analysis was the polygons of the COS07 cartography within the study area that belong to the broad LULC categories shown in Figure 1, namely agriculture, agroforestry, forest, and shrubland (n = 149013). We considered the lowest legend level available for COS07 LULC classes, which reports 59 different LULC classes under agriculture, agroforestry, forest, or shrubland categories (see Supplementary Material Table S1).

Ecosystem Condition (EC)
Four indicators were considered to assess ecosystems' condition: Soil Organic Matter, Ecological Value of Plant Communities, Plant Diversity, and Bird Diversity. A summary of the selected indicators is presented in Table 1. The methods are further detailed in the following subsections. Soil Organic Matter content was assessed based primarily on the information presented in the National Greenhouse Gases Inventory Report (NIR), according to its land-use typology (Kyoto Protocol classes) although minor adjustments have been introduced (i.e., changes in organic matter estimates in areas undergoing land-use change). Soil organic matter is indicative of the ecological condition of soils, being essential to maintaining soil ecosystem functions such as stabilization, water infiltration, and conservation of nutrients.

Ecological Value of Plant Communities
Semi-Quantitative Score (1 to 5) The Ecological Value of Plant Communities represents the mean value of five parameters (naturalness, replaceability, threat, rarity, and condition), scored from 1 to 5, was attributed to each of the studied ecosystems (level n). The geobotanical models used, at the geographical scale in which they were implemented, are indicative of the ecological condition of ecosystems by providing integrative information on the structural quality, phytocoenotic integrity, and successional maturity of the present plant communities. This Indicator is thus given by the estimated number of species in grid cells (2 × 2 km) covering the study area, which was reclassified into a 1 (low bird diversity) to 5 (high bird diversity) scale. as can be noted, this differs from the unit of analysis of the other indicators (the LULC polygons from COS07), but this issue has been properly addressed when accounting for spatial relationships. Birds have been widely acknowledged as indicators of the ecological condition of forests and agroecosystems, with bird diversity being one possible good measure of the general ecological condition and overall biodiversity present in an ecosystem.

Soil Organic Matter
To map soil organic matter, we used average organic carbon content present at 0-40 cm deep soil ( Table 2)   The specific land-use classification adopted in the NIR (Kyoto Protocol land-use classes, Table 3) was translated into land-use typologies (COS07) for harmonization of legends and further use in our study. This harmonization (Supplementary Material Table S4) was later evaluated and approved by NIR authors in a meeting. We also considered land-use transitions observed (1990 to 2007) to properly adjust the average soil organic carbon content present in each polygon. This means that the soil organic carbon content in polygons that underwent a land-use transition → is given by the average of the content value of and . All calculations were performed in ArcGIS (v9.3).

Ecological Value of Plant Communities & Plant Diversity
These two EC indicators were assessed as presented in [21]. The indicator Ecological Value of Plant Communities represents the mean value of five parameters (naturalness, replaceability, threat, rarity, and condition), scored from 1 to 5, attributed to each of the studied ecosystems (level N). The geobotanical models used, at the geographical scale in which they were implemented, integrated information on the structural quality, phytocoenotic integrity, and successional maturity of plant communities. Further details on the methods used are provided in [21].
The Plant Diversity indicator was assessed assuming that vegetation series maps provide reliable information on the natural communities occurring at different locations. It was thus possible to consult phytosociological tables of these communities and to know their average or characteristic floristic composition, which reflects species richness and rarity, as well as the presence of endemic or threatened species. Based on 3500 phytosociological inventories, representative of Portuguese natural vegetation, plant diversity was estimated as the weighted average of four different parameters attributed to each plant community (presence of protected species, of other endemic species, of other rare species, and of characteristic species). Further information on the methods used is provided in [21].

Bird Diversity
As birds are usually referred to as an "indicator" group for several environmental parameters [22], including biodiversity and the condition of ecosystems, bird diversity was also chosen to assess EC. We performed a multiple logistic regression with presence/absence records of farmland and forest bird species and a set of explanatory variables (land-use typology, topography, temperature, and rainfall) ( Table 4). See Supplemetary Material Table S1. for aggregation criteria for land-use typologies. Bird data was collected from an extensive and publicly available dataset of observation records (from the PortugalAves database, which now has been merged into eBIRD https://ebird.org/portugal/home). Data processing and validation was required to select only records of true absence/presence and corresponding to a specific temporal scale (2004-2011) close enough to the land use cartography used (2007). This resulted in a final list of 50 species and 15,000 observations (Table 5). Bird observations were recorded in 2 × 2 km cells, and thus we applied a 2 × 2 km grid over the study area and used the cells as the unit of analysis for the Bird Diversity indicator. When comparing with other indicators (see Section 2.5, we integrated the results into the LULC polygons from COS07 by averaging grid cell values within each polygon. We ran a separate multiple logistic regression model for each bird species, which resulted in the probability of observation of each species at any given cell of our study area. Ten out of the fifty species were excluded at this point since the models returned significant errors (either false convergency or null convergency, see marked species in Table 5).
All calculations were performed in R (v. 2.15.13). The multiple logistic regression model was adjusted with the glm (family = binomial logit) function, with step (glw = backwards). AIC (Akaike's information criteria) values were compared before and after stepwise adjustment to ensure the appropriate selection of variables. The Hosmer-Lemeshow goodness-of-fit test [23] was used to compare the log-likelihood of observed and estimated values through the hoslem.test(fitted(backwards)) function (library ResourceSelection). Results for each model are provided in Supplementary Materials (Table S8). We mapped results considering that a species would occur in a given cell if the probability returned by our model for that cell was over 50%. This resulted in a map with the potential distribution of bird diversity in the study area. We reclassified our results on a scale from 1 (low diversity) to 5 (high diversity) ( Table 6).

Ecosystem Services (ES)
Five ES were quantified and mapped: control of erosion rates, climate regulation through carbon sequestration, fiber production, crops, and extensive livestock production. A summary of the selected ES, including its designation under the CICES [2] classification system, is presented in Table 7. The methods are further detailed in the following subsections. Control of Erosion Rates was modeled and mapped based on the Universal Soil Loss Equation (USLE), integrated into a GIS platform, which allowed determining the difference between erosion rates in the current scenario (i.e., erosion rates given actual land cover type) and erosion rates for a worst-case scenario (considering a maximum erosion cover type), as first suggested by [24] 2.4.

Control of Erosion Rates
Control of Erosion Rates was estimated considering the contribution of a given soil cover to reduce soil erosion by comparing it with a worst-case scenario (i.e., the land cover that would generate the highest erosion rate at a given point). Soil erosion rates were estimated using the Universal Soil Loss Equation [25], which consists of a tier 3 approach. The equation estimates average soil loss ( ) such as: As no information was available to estimate factor P accurately, it was kept constant as 1. For our approach, we considered that the Control of Erosion Rates service was given by avoided erosion ( ) when factor =1 (soil without any natural cover), which can be re-written as: Factor K was obtained by combining information from national soil cartography available at http://sniambportal.apambiente.pt and national literature. European estimates for Factor C [26] were refined per Portuguese LULC classes based on [27] (see Supplementary Material Table S2). Factor R was based on national statistics publicly available at http://geo.snirh.pt/AtlasAgua/(required conversion of American to SI units based on [28]. Factors LS (combined) are given by: All calculations were performed in ArcGIS using raster calculation operations. We selected the finest resolution possible for the raster cells (250 × 250 m), but, by definition, when performing raster calculations our results were not mapped using the LULC polygons from the COS07 cartography as unit of analysis (same as explained for Bird Diversity). However, when depicting spatial relationships (see Section 2.1), we averaged raster cell values per LULC polygon to make results comparable.

Climate Regulation through Carbon Sequestration
Carbon Sequestration was estimated as the balance of gains and losses of carbon in both biomass (above and belowground) and soil, considering the land-use transitions that occurred between 1990 and 2007 (and assuming these transitions occurred at a constant rate). IPCC guidelines point to a 20-year transition period to be considered for carbon balances, for the stabilization of carbon fluxes that are slow and/or occurring between two different states [29,30]. Owing to the LULC cartography available (from 1990 and 2007), the analysis was limited to a period of 17 years.
For estimating carbon balances, we first determined carbon gains by calculating mean annual biomass increments as reported in NIR (Tables 8 and 9), converted into tonC.ha −1 .year −1 using the appropriate conversion factors (root-to-shoot, biomass expansion factor, and carbon fraction) (see Supplementary Material Table S3.). Next, we determined carbon losses by using data on timber harvesting, observed land-use transitions (from 1990 to 2007), maps of burned areas, and information on natural mortality to determine: • Mortality by timber harvesting (spatialization of statistical data allowed determining harvesting rate for the given years at each polygon of interest, i.e., polygons with timber harvesting plantations).

•
Mortality by fire (total loss of biomass in polygons that experienced fire events, based on official fire maps for 2007, available at http://www2.icnf.pt/portal/florestas/dfci/inc/mapa).

•
Mortality by transition (total or partial loss of biomass due to land-use transition, from 1990 to 2007, observed in a given polygon).

•
Natural mortality, as determined and reported by the NIR (Table 10).  All information was converted to tonC.ha −1 .year −1 to allow comparisons. The Carbon sequestration service was given by the positive carbon flows estimated (above 0 tonC.ha −1 .year −1 ). Data was managed in a specific geodatabase, and calculations were performed in ArcGIS and Access (Microsoft Office). Emission/sequestration coefficients were thus obtained for each land-use transition, resulting in a carbon emission/sequestration map.

Provisioning ES
Fiber Production was estimated as mean annual increments of forests trees as presented in the NIR (reported for KP forest classes, see carbon sequestration above), deducing biomass losses due to natural mortality. Spatialization of ES supply was possible after harmonization of legends between KP classes and COS07 (see Supplementary Material Table S4.) Mapping of crop production was based on the establishment of a correspondence between the main crops in the study area (as listed in the official statistics) and COS07 classes. We first listed and systematized the main crops, areas (ha) and productivity (ton/ha) in the study area as extracted from the statistics bureau. Once the correspondence with land use classes was made (in certain cases, more than one crop was attributed to the same land-use class; only "pure" land use classes were considered), average productivity was estimated for each land-use class (COS07)-see Supplementary Material Table S5 for established relationships.
Extensive livestock production (or the capacity of ecosystems to support extensive livestock production) was quantified and mapped by determining average livestock densities (for calves, dairy cattle, and sheep) in pasture areas within the study region, using official national statistics at the municipality level. As it was not possible to geographically identify pastures where each unit of livestock production occurs, average livestock density was estimated considering the two main species together (cattle and sheep) for the given pasture area in national statistics for each municipality. Next, average livestock values were spatialized by attributing them to meadow and pasture areas as defined by COS'07. Due to the identified disparities between pasture areas reported in COS'07 and pasture areas reported in the national statistics, the resulting values presented in the map should not be aggregated.

Spatial Relationships and Interactions
The analysis of spatial relationships and interactions between ecosystem condition (EC) and the supply of ecosystem services (ES) is intended to explore the outcomes of the MAES initiative envisioning its potential contribution to landscape planning in the region. This analysis was conducted from two different perspectives: the assessment of statistical correlations and the analysis of the spatial distributions and relationships between EC and ES.
The models used for bird diversity and control of erosion rates required mapping using raster/grid cells of finer resolution as the unit of analysis (as detailed in the corresponding section), and here they were analyzed along with the other variables at the polygon level by averaging cell values per LULC polygon of the COS07 cartography.
For the statistical correlations, we first considered five levels for each EC indicator. We took the semi-quantitative 1 to 5 scale already assessed for bird and plant diversity, as well as the ecological value of plant communities. For the soil organic matter indicator, we categorized the values into five levels using Jenks/Natural breaks [31]: very low (0-34 tC/ha), low , medium (56-69), high (69-82) and very high (82-113). These levels are in accordance with the range of values for soil organic matter previously reported in Alentejo [32]. We then analyzed ES supply for each EC level by applying the Kruskal-Wallis rank test [33] on the EC levels' medians to detect significant differences of ES supply between each level (i.e., to evaluate whether different levels of EC result in different ES supply). We also performed a Jonckheere-Terpstra test [34,35] to identify positive or negative trends in the relationship between the supply of each ES and EC level (i.e., to evaluate whether increasing/decreasing EC levels result in higher/lower ES supply). The statistical work was conducted in RStudio (version 4.0.5).
To show an integrated overview of the distribution of ES supply and EC levels in Alentejo, we normalized all the variables into a 0 to 100 scale to make them comparable. To spatially visualize the relationships between the ES and EC indicators, we overlapped the normalized ES supply and the normalized sum of all EC indicators. We also indicate the relationship per LULC typology, differentiating between agriculture, agroforestry, and forest/shrubland classes. This analysis allowed us to identify how EC condition is related to ES supply in Alentejo.

Results
The results of the assessment are presented in maps showing the distribution of the indicators for ecosystem condition (EC, Figure 3) and the ecosystem services supply (ES, Figure 4) within Alentejo. The breakdown of results per ecosystem type is summarized in Supplementary Material Table S6. Other graphs and maps show the spatial relationships and interactions between EC and ES ( Figures 5-9).

Ecosystem Condition (EC)
Soil organic matter mapping showed average organic matter content in Alentejo at around 50-75 tC/ha (Figure 3a). Only 12% of the study area presented OM content above 85 tC/ha. OM concentrations were overall higher in the north and southwest of the region, mainly in forest and shrubland ecosystems. Results for the EC indicator ecological value of plant communities (Figure 3b) shows most of the region presented lower scores (1)(2). Areas with the highest score (5) can be sparsely found along the west coast (in dune/paleodune and sea cliff ecosystems with a higher preservation value or presence of successional evolution) and to the northeast and southeast of the region (in Quercus suber dominated agroforestry ecosystems with low agricultural intensity). Higher plant diversity was found in the center and southwest of Alentejo (Figure 3c), mainly in agroforestry and forest ecosystems characterized by sandy soils. As for bird diversity (Figure 3d), over 13% of the study area presented the highest score (5). Higher bird diversity scores were found mainly in broadleaved deciduous forests, temperate Mediterranean scrubs, and in inland rock cliff ecosystems.

Ecosystem Services (ES)
In terms of control of erosion rates (Figure 4a), higher service supply (over 50 t/ha/year) was found along mountainous regions to the southwest and northeast of the region. Grasslands and forests (Quercus suber forests) are the most erosion controlling ecosystems. Carbon sequestration was determined using carbon balances, and results show positive rates (sequestration) mainly present throughout the region but predominantly in the southwest (Figure 4b). Negative rates (carbon emission) dominate the eastern part of Alentejo. Higher service supply (i.e., carbon sequestration rates) are associated with transitions to forest ecosystems, namely areas converted to broadleaf forests and Pinus pinea, while the highest emissions rates were found in areas that were converted from forest to non-forest uses over the 17 years analyzed (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007).
Three provisioning ES were assessed based on the mean annual increment of forest species of interest (fiber production) and spatialization of statistic data reports (crop and extensive livestock production). Spatialization of data was limited to specific land-use typologies (e.g., fiber production is not accounted for in grasslands or croplands), which explains the larger areas where ES supply is absent ("0") in these maps (Figure 4c-e, respectively). Fiber production in Alentejo is mostly present at lower rates in large extents of agroforestry ecosystems and Quercus suber and Quercus rotundifolia forests, mainly in central Alentejo (Figure 4c). Nevertheless, other ecosystems dominated by fast-growing species (Eucalyptus sp. and Pinus pinaster) can be found in the eastern part of the region, where much higher rates of service supply were found (up to 9.4 m 3 /ha). As for crop production, a higher service supply can be found along the Tagus Valley (northwest) and sparsely throughout the region where cropland is present. Finally, results for livestock production shows the level of service supply ranging from 0.2 to over 1LU/ha. Levels higher than the 0.6-0.8 ha category are found mainly in grassland ecosystems, whereas lower levels are present in agroforestry ecosystems (extensive grazing).  (Figures 8 and 5, respectively). Contrastingly, the supply of CE and CS seems to decrease with increasing plant diversity and ecological value of plant communities (Figures 7 and 6, respectively).

Spatial Relationships and Interactions
As for provisioning ES, supply levels mainly decrease with increasing levels of EC indicators, except for fiber production (which increases with increasing organic matter, Figure 5) and extensive livestock production (which increases with increasing ecological value of plant communities, Figure 6). These results were expected since (1) higher levels of organic matter were found in production forests and (2) the higher supply of extensive livestock refers mainly to medium range values (around 0.5 LU/ha), which occur in agroforestry systems that scored medium-high in the ecological value EC indicator.
Additionally, we identified overlaps between ES supply and EC indicators. Figure  9a-c show the distribution of normalized ES supply for fiber, crop, and livestock production, respectively, in relation to the normalized sum of all EC indicators (ecological value of plant communities, plant and bird diversity, and soil organic matter). Overlap results for the two regulating ES (control of erosion rates and carbon sequestration) are shown in Supplementary Material Table S7. We also present the breakdown of the overlap between ES and EC per land-use/land cover (LULC), indicating which LULC class is more predominant under different combinations of ES supply and EC level. For instance, in Figure 9c, we see that intermediate levels of EC condition and ES supply are more predominant in agroforestry areas (yellow dots) than in pastures/grasslands (green dots). Areas of high ES supply and below-average EC level (green shades on the top left of the legend) can be found throughout Alentejo for the three provisioning ES, particularly to the north and northeast of the region (Figure 9a,b) but also sparsely in the central part of Alentejo (Figure 9b,c). This indicates a high supply of ES in areas of poor ecological condition.

Discussion
Our discussion is framed by the objectives earlier stated: first by considering the proposed analytical framework and outcomes of the ptMAES study in the context of future MAES initiatives, and second by discussing the assessed relationships between ecosystem services (ES) and ecosystems' condition (EC).

Proposed Analytical Framework
We identify several advantages and challenges of the proposed analytical framework for implementing national policies targeting the EU Biodiversity Strategy to 2020 and 2030, as well as assisting the launch of MAES initiatives in Portugal.
The methods proposed for assessing the ES control of erosion rates bring novelty to the soil arena, being a first attempt to detail a C-factor for the national land-use cartography classes (COS07) and a GIS-based application of the USLE at a regional scale in a context of ES assessment. This resulted in a highly detailed map (Figure 4a), capturing spatial variability with the potential to assist landscape planning at multiple scales. Our estimates proved to be more detailed and reliable than soil erosion maps produced at the European scale (JRC PESERA) based on remote sensing/LULC data, although our range of estimates falls within the reported values for Portugal. However, the proposed approach fails to incorporate the temporal variability of ES supply and the impacts of different land management options, which can impact observed erosion rates [23].
The approach for mapping and assessing the ES carbon sequestration (Figure 4b) presented several advantages: it uses information already produced for national reporting (UN Framework Convention on Climate Change and Kyoto Protocol), innovatively applied on a spatial-based approach, and it accounts for LULC changes (i.e., temporal and spatial dimension), allowing the ES to be genuinely quantified as flows [36]. Accounting for LULC change is critical as its potential impact on ES supply at different scales has already been called out in the literature [37][38][39], and LULC changes in Alentejo are significant and mainly driven by agricultural policies [40].
Though the use of regional process-based models could significantly improve outcomes [41], they provide additional challenges in terms of data requirements that might not be readily available for the scale of interest.
Lack of disaggregated data (i.e., at the farm or local scale) publicly available was also a significant limitation observed in the assessment of provisioning services (Figure 4c,d,e). Other authors have reported this to be a common limitation of tier-1 approaches for ES mapping [42,43], which is usually the case when mapping provisioning services. In the case of fiber production, this limitation may be solved by using more spatialized annual increment coefficients [44]. This refinement can also be possible based on literature and/or expert opinion while keeping the cost-efficient nature of the proposed framework [45]. Better results can be achieved by integrating Land Parcel Identification System information, which allows for a more detailed spatialization of crop and livestock production [46]. However, using this information in research agendas requires disclosure of personal data and raises confidentiality and data protection issues.
Concerning EC indicators, the advantages of the proposed approach to mapping soil organic matter (Figure 3a) are the same as reported for carbon sequestration, in the sense that it relies on information that is available and collected for national reporting on GHG emissions (UNFCC and Kyoto Protocol), but innovatively spatialized based on LULC cartography while accounting for LULC change. The estimates efficiently account for spatial variability and have the potential to assist landscape planning at multiple scales while falling within the range of values previously assessed based on regional models using sampling data [32]. We also highlight the potential of the proposed biodiversity indicators, namely plant diversity and ecological value of plant communities (Figure 3c and Figure  3b, respectively), in assisting decision making regarding landscape planning (e.g., protected areas management), as they focus on ecological and functional measures complemented with expert judgment rather than on the presence or absence of species alone. We argue that such integrated solutions conveying expert knowledge to mapping could better improve the assessment of biodiversity-related indicators with increased potential to assist decision-making at multiple scales [47,48]. As for the bird diversity indicator (Figure  3d), although the method proposed is suitable for upscaling, rare species with conservation interest were not included in the analysis due to their low observation frequency in the study area. This limitation should be addressed in future applications if the purpose of the assessment is to support nature conservation policies targeting protected species [49]. Moreover, it has been argued that high biodiversity levels in managed landscapes are more likely to be maintained for reasons of intrinsic value (e.g., traditional land management practices), cultural values (e.g., "bequest" ES) or use values ("direct use") than for its functional values or its role in maintaining good ecosystem condition [50]. Our findings in this paper support that it is meaningful to include biodiversity assessments both as ES supply (biodiversity with conservation interest) as well as EC indicators (functional biodiversity) in MAES initiatives aimed at the integration of ES into planning/policy design in human-dominated landscapes. Finally, when assessing EC, we argue that priority should be given to indicators that could be related to the ES flow measurement (e.g., soil organic matter for crop production, infiltration rate for water availability, etc.), such as demonstrated by recent findings [51,52].
Overall, the proposed framework can be easily applied at different scales subject to specific refinements dictated by territory spatial variability. This is clearly a major advantage of performing a multi-tiered spatial assessment based on LULC cartography as proposed in the ptMAES study, which has been also pointed out by others [45], However, we acknowledge that in the absence of data limitations (in terms of scale, availability, and coverage), a more comprehensive assessment would be possible, including a larger number of ES and a wider array of EC indicators. Cultural ES, in particular, should not be overlooked as they play an important role in motivating public support and assisting decision-making [53,54]. To make the best out of the available information, certain ES and EC indicators were based on tier-1 approaches (i.e., spatialization of statistical information) (Figure 2), which limited the portrayal of spatial variability and consequent explanatory power of the assessment to better assist decision-making [52]. These issues have been discussed by [55], and our analysis supports the review and directions proposed by the authors. We highlight that the ptMAES also produced a thematic EUNIS habitat cartography based on interpreting LULC units as territorial units of ecological succession. This was an innovative approach in Portugal, and the outcomes are presented in [22]. Results of EC and ES supply per EUNIS habitat are presented in Supplementary Material Table S6.

Spatial Relationships and Interactions
In this paper, we analyzed the spatial relationships and interactions between EC and ES supply in Alentejo. Our goal was to explore the outcomes of the analytical framework proposed and discuss how they could potentially support landscape planning.
Our analysis indicates overall changes in ES supply rate for varying levels of EC (Figures 5-8). Generally speaking, we expected regulating ES supply to increase with increasing EC (as ecosystems under better ecological conditions potentially sequester more carbon and provide more soil protection) and provisioning ES to decrease with increasing EC (as higher production rates are usually found in monocultural crops or intensive pastures under poorer ecological condition) [56]. Despite this, we found that carbon sequestration rates decrease with increasing levels of plant-related EC indicators (Figures 6 and 7). Two underlying factors explain this result: (1) the great representation of agroforestry systems in good EC (scoring 3 and higher, see Figures 5b and 6b); and (2) the lower carbon sequestration rates attributed to Q.suber/Q.rotundifolia agroforestry systems, which are given by the lower mean annual increment considered for these species in the estimated biomass balance, when compared to other forest species. This means that carbon sequestration is potentially underestimated in these ecosystems and is not true a reflection of their ecological condition [43]. In addition, we found crop and fiber production to increase with increasing EC levels (particularly bird diversity, Figure 8a). We believe this is caused by: (1) the consideration of common farmland and arable land bird species in the distribution models, which estimate high bird presence in cropland and harvested forests alike; and (2) fiber and crop production levels are not real estimates collected in situ but rather result from the specialization of statistical data. This result should be interpreted with caution since bird presence is actually expected to decrease in intensive monocultures [57].
Despite the considerations presented above, when looking at spatial overlaps between the supply of provisioning ES and the normalized sum of all EC indicators ( Figure  9), we identify high ES supply rates in areas under below-average EC (green shades in the maps). Since good ecological conditions ultimately underpin the capacity of an ecosystem to supply ES [58,59], these results red-flag areas where the sustainability of ecosystems (and its capacity to supply ES) is threatened in the long term. The breakdown of these results per LULC typology shows that agroforestry land-uses promote higher levels of ES supply (namely crop and extensive livestock production) under best (above-average) ecological conditions (see, for instance, the representation of agroforestry in Figure 9b and Figure  9c, respectively). This is less evident for fiber production as higher production levels are limited to fast-growing species (Eucalyptus sp. and Pinus pinaster), overshadowing fiber production in Q.suber-/Q.rotundifolia-dominated agroforestry systems. However, the bundle of ES supplied by the agroforestry systems present in Alentejo has been widely documented in the literature [60][61][62][63].
Our outcomes show examples of how MAES could potentially assist landscape planning. For instance, our results point to planning interventions that could: (i) target and help improve management practices in croplands, timber forests, and grasslands where ES supply is high, but EC is not optimal, and (ii) target agroforestry areas in good EC and which have the potential to supply many ES in the central part of Alentejo. Notwithstanding, we point out that our achievements are specific to this case study and should not be readily transposed to other research or policy contexts.

Conclusions
In this paper, we present an analytical framework based on LULC cartography to assess ecosystem condition (EC) and ecosystem services (ES) as tested in the ptMAES study. Innovative, multi-tiered information-processing approaches were used, which proved to be effective and suitable for upscaling. Our findings leverage new approaches to the potential integration of ES into landscape planning at the regional and national scales. We highlight the usefulness of integrated approaches (e.g., conveying expert knowledge to spatial assessments), as it can better depict spatial variability in the absence of data available at the desirable scale. However, we identify relevant caveats in the proposed analytical framework that should be addressed in the future for a more comprehensive MAES assessment: (1) data availability (in terms of aggregation, scale, and coverage) limited the inclusion of process-based modeling; (2) refinement of results could be achieved with the use of information collected by the public administration if data protection issues are overcome; (3) a wider range of EC indicators and ES should be considered (particularly cultural services); and (5) the selection of EC indicators should better reflect their relationship to ES.

Acknowledgments:
The authors would like to thank the extended ptMAES working team who greatly contributed to the production of the original data used in this study: Sandra Mesquita, Jorge Capelo, Ivo Gama, Miguel Alves, Vânia Proença, Paulo Canaveira and especially Marco Reis, who, although no longer with us, continues to inspire by his example and dedication to the people he worked with over the course of his fruitful career.

Conflicts of Interest:
The authors declare that they have no conflict of interest. 1 The 2020 EU Biodiversity Strategy (COM 2011) was built around six mutually supportive and inter-dependent targets that addressed the main drivers of biodiversity loss. They aimed to reduce key pressures on nature and ecosystem services in the EU by setting up efforts to fully implement existing EU nature legislation, anchoring biodiversity objectives into key sectoral policies, and closing important policy gaps. Each target was accompanied by a set of focused, time-bound actions to ensure these ambitions are fully realized. 2 The goal of Target 2 of the EU Biodiversity Strategy 2020 is to "maintain and restore ecosystems and their services", with Action 5 set out to "improve knowledge of ecosystems and their services in the EU". 3 NUTS II refers to the second level of the Nomenclature of Territorial Units for Statistics (NUTS) that is used in Portugal.