Effectiveness of conservation areas for protecting biodiversity and ecosystem services: a multi-criteria approach

ABSTRACT In land planning strategies, methods to quantify ecosystem services (ESs) are now used to complement biodiversity assessments. Tension arises when areas important for biodiversity do not spatially co-occur with important areas for ESs. We investigate the effectiveness of protected areas in simultaneously protecting biodiversity and ESs in central Colombia and identify complementary areas. We map, integrate using a multi-criteria technique and correlate five indicators (sensitive species, ecological systems, habitat quality, scenic beauty and water provision). Reflecting the uncertainty in criteria weights, multiple maps were created and overlain with current protected areas to investigate their effectiveness. A consensus and an uncertainty map were calculated to identify multifunctional areas (high value for biodiversity and high provision of ESs and low uncertainty). Protected areas show low to intermediate levels of effectiveness (3–56% percentage overlap with simulated areas), with water provision being the service least effectively protected. Indicators do not show high levels of correlation (all p < 0.57). Sensitive species are negatively correlated with all other indicators. High representativeness levels were found around the city of Bogotá; still, extensive multifunctional areas are not contemplated under any protection status. We advocate the use of our approach to identify multi-purpose areas that are robust to divergent stakeholder opinions. EDITED BY Berta Martin-Lopez


Introduction
Current strategies for environmental conservation have gradually changed their focus from purely biodiversity and landscape-scale management to an ecosystem services (ESs) approach (Egoh et al. 2007;Chan et al. 2011a;Albert et al. 2014). This tendency is reflected by an increase in the number of research projects and publications dealing with ESs. These projects and publications have, among other things, developed tools, methods and models (Bagstad et al. 2013;Mulligan 2013;Peh et al. 2013;Sharp et al. 2014), have promoted ESs as a basis for area prioritization (Chan et al. 2006;Egoh et al. 2007;Luck et al. 2012), and have accounted for ESs in land planning and management (Chan et al. 2011b;Albert et al. 2014Albert et al. , 2016Sitas et al. 2014). International agreements (i.e., The Aichi Targets of the Convention on Biological Diversity -CBD), assessment initiatives (i.e., the Millennium Ecosystem Assessments -MA and The Economics of Biodiversity and Biodiversity -TEEB) and government interest to close the gap between science and policy (i.e., the Intergovernmental Platform on Biodiversity and Ecosystem Services -IPBES) have followed.
As a member of the CBD, Colombia is committed to fulfilling the Aichi Targets which recognize the importance, maintenance and conservation of ESs for the benefit of human well-being. Following the CBD agreement, the incorporation of the ESs approach into the development of new environmental policies, land planning guidelines and decision making tools has been growing rapidly during the last few years. This development is evidenced by the explicit integration of the ESs approach into the National Development Plan (DNP 2015), the recent release of the National Policy for Biodiversity and Ecosystem Services Management (MADS 2013) and the incorporation of ESs analysis in land planning procedures. In these procedures, ESs are meant to contribute to the definition of the 'environmental conservation and protection category' contemplated in the zonation process, for example, see the 'Technical guide for water basin management' (MADS 2014a). This new category is meant to group together current protected areas, ecosystems of special interest, areas with high levels of ESs provision as well as other key biodiversity indicators (e.g., endemic, threatened and migratory species). The consideration of different key biodiversity indicators is based on the evidence that areas of high species richness are not necessarily the most effective for conserving other species with particular characteristics, like endemic (Orme et al. 2005) or migratory species (Faaborg et al. 2010).
Examples of management strategies and decision support in land use planning for ESs are very scarce in comparison to those tackling biodiversity or landscape conservation in the conventional sense (Geneletti 2013;Hodder et al. 2014). In Latin American countries, the current knowledge on ESs, their quantity, distribution and valuation is particularly limited (Balvanera et al. 2012) and Colombia is no exception (Ruiz-Agudelo & Bello 2014). This is partly due to weak governmental and local environmental authorities and partly due to the lack of expertise in applying technical and methodological frameworks already available (Ruiz-Agudelo 2014). The ability of current conservation and protection areas to protect both key biodiversity indicators and ecosystems services in Colombia has never been evaluated.
The new approach begs the question whether conservation and protection areas established by traditional methodologies (i.e., by expert knowledge, based on a particular ecosystem of interest, based on few biodiversity or services indicators like umbrella species or water service, respectively) simultaneously protect areas of high value for biodiversity and provision of ESs, and which areas might be considered in addition.
Research on this topic, at various spatial scales, has shown that the protection of ESs within the existing networks of protected areas is important in some cases but still limited (Cabello et al. 2013;Durán et al. 2013;Hodder et al. 2014;Palomo et al. 2014;Castro et al. 2015). In general, these studies have demonstrated the potentially added benefits of accounting for biodiversity and ESs in conservation planning, and the need to develop and apply prioritization frameworks and methodologies specifically for that purpose (Luck et al. 2012). Closely related are those approaches aiming to quantify the spatial cooccurrence of different ESs (Bryan et al. 2011;Rodríguez et al. 2015) and between ESs and biodiversity (Egoh et al. 2009;Rogers et al. 2010;Dickie et al. 2011;Armenteras et al. 2015). The prevailing findings are that spatial correlations between different ESs and between ESs and biodiversity are not necessarily high, and that there is always a trade-off in which focusing on the protection of one service will compromise the protection of others.
One of the most commonly used methodologies in land planning is spatial multi-criteria analysis -(sMCA) (Ciarleglio et al. 2009;Bottero et al. 2013;Geneletti 2013). Land use types or areas suitable for a specific purpose are defined by the superposition of different information layers, each describing a degree of suitability. The procedure involves the selection, quantification and mapping of different criteria, data standardization and the definition of a weight for each criterion describing its relative importance. Typically, these weights are defined by consensus with stakeholders, sidestepping considerations of bias associated with personal interests, established stakeholder hierarchies, and political agendas and power relations. Hence, assigning criteria weights has been recognized as a major source of uncertainty in model outputs of multi-criteria analysis (Store & Kangas 2001) and is an active field of research (Ligmann-Zielinska & Jankowski 2014).
In this study, we investigated the effectiveness of the current network of protected areas in simultaneously protecting biodiversity and ESs in the Andean high planes of the eastern mountain ranges in Colombia, and identified potential complementary priority areas for conservation. We applied a multicriteria approach with explicit consideration of the uncertainties in criteria weighting. Our case study holds practical lessons for the inclusion of ESs and uncertainties in local land planning.

Study area
The study area is the Corporación Autónoma Regional de Cundinamarca (CAR) jurisdiction located in central Colombia in the east mountain range (cordillera oriental) of the Andes (Figure 1). Elevation varies from 4000 m.a.s.l. on the east side where the páramos ecosystems can be found, to 150 m.a.s.l. on the west side reaching down to the Magdalena River. From the centre to the east side of the study area stretches a large flat terrain ('high plateau') where the capital city of Bogotá is located (2500 m.a.s.l.). The CAR is the regional authority in charge of managing the natural resources of the rural areas and ecosystems surrounding Bogotá, all of them part of the Magdalena river basin. The total population in Bogotá and surrounding cities is about 12 million. The study area covers 19,411 km 2 of which 18% (3679 km 2 ) are currently conservation and protection areas (cCPAs) comprised of national, regional and local protected areas and key ecosystems (i.e., páramos and wetlands). The density of cCPAs is high in the eastern part of the study area around the city of Bogotá and along the high plateau, its surrounding hills and mountains. In contrast, the western part of the study area, characterized by sloping down more than 2000 meters to the Magdalena River, has very few areas with any protection status.

Biodiversity and ecosystem services indicators
Indicators were selected, and in some cases parameterized, during a scoping workshop with technical actors (Conservation International and CAR project members: hereafter referred to as 'project members'). The selection was based on the quality and availability of spatial information and on CAR's short-term environmental management priorities (i.e., water, ecosystems and biodiversity). The five indicators selected were considered a good starting point to evaluate the degree of representativeness when ESs are considered in prioritization approaches. The indicators are described as follows.
2.2.1. Sensitive species 'Sensitive species' is a term that groups together endemic, migratory and endangered species. A database was constructed from museum collections, publications and online databases. The taxonomical nomenclature was carefully checked and corrected following the most up-to-date authorities (Uetz & Etzold 1996;SELVA 2012;Solari et al. 2013;Frost 2014;IUCN 2014;MADS 2014b;Remsen et al. 2015). The indicator used to represent sensitive species was their richness (i.e., the number of sensitive species in a 1 km 2 cell). In total, 75 sensitive species were considered (41 birds, 14 mammals, 14 amphibians, 6 reptiles) (see Supplementary Material I). Species richness of all sensitive species is a well-established criterion of several environmental frameworks in Colombia. For consistency with other national efforts, we grouped all species into one single indicator. However, we acknowledge that the richness patterns of each taxon can be different and, therefore, using them individually might generate a different set of complementary areas and values of representativeness.
The records (i.e., geographical locations) of each species in the database were used to estimate their distribution ranges through species distribution models (SDMs). These models combine the locations of each species with the values of a set of environmental variables in those locations, quantify the relationships and extrapolate an index of habitat suitability or presence probability over the study area (Guisan & Zimmermann 2000). Species with less than five records were removed from the analysis since it has been shown that at least five records are needed to model distribution ranges accurately (Pearson et al. 2007). A database of environmental variables was created consisting of climatic (Hijmans et al. 2005) and terrain variables derived from a digital elevation model (DEM) (USGS 2008). Only variables with a Pearson correlation coefficient with other variables of less than 0.7 were chosen for the analysis. These were elevation, annual temperature range, mean temperature of the driest quarter of the year, annual precipitation and precipitation of the driest month. Support vector machines (SVMs) (Cortes & Vapnik 1995) were used to calculate the distribution ranges of the species. SVMs have been successfully used to model the distribution of species (Guo et al. 2005;Drake et al. 2006;Heubes et al. 2013) and, as a machine learning method, have performed better than other traditional algorithms (Elith et al. 2006). SVMs employ optimization algorithms to locate the boundaries between classes (presence and absence of species) by finding a hyperplane in the feature space that maximally separates these classes (Cristianini & Scholkopf 2002). Absence data were represented by randomly selecting locations outside of the multivariate environmental similarity surface (MESS) (Elith et al. 2010) of the presence data. The mean of a 10fold cross-validation was applied both to measure model performance (mean squared error) and to calculate the ʋ parameter (percentage of outliers in the training data set). The kernel used was radial. The SVMs were implemented in the software environment R (R Core Team 2014) making use of the 'e1071' library (Meyer et al. 2014).
Finally, the richness of sensitive species was calculated by summing the species probability values obtained for each cell. The richness value thus obtained for each cell is the probability-weighted sum of species present. This method was applied because for the analysis, it was more important to have the distribution patterns of richness rather than their absolute numbers. In addition, we avoid the uncertainty of finding a specific threshold for each species to transform the probability values into a binary presence/absence response (Liu et al. 2005; Jimenez Valverde & Lobo 2007).

Ecological systems
As part of the development of a new biodiversity offsetting strategy in Colombia (MADS 2012), the government with the support of The Nature Conservancy (TNC) made use of the concept of 'ecological systems' (i.e., landscape units resulting from combining ecosystems and biogeographical units) to rank land based on four different criteria (Saenz et al. 2013): representation (amount of the ecological system within protected areas), rarity (the national-level and local-level rarity of the ecological system), remanence (the percentage of the ecological system that remains relative to its historic distribution) and rate of loss (the rate of loss of the ecological system calculated over the previous six years). This ranking, or offset ratio calculation (Saenz et al. 2013), is used in the present study as the prioritization value of each ecological system where a value of 10 is the maximum for an ecological system with very low representation, that is very rare, has low remanence and a high rate of loss. The ranking ranges from values of 10 down to 5 for 'natural' ecological systems and from 5 to 2.5 for ecological systems with secondary vegetation (i.e., vegetation that has recovered and is no older than 15 years). Ecological systems with no natural or secondary vegetation are assigned a value of zero (see Supplementary Material II).

Habitat quality
Habitat quality is considered here as an indicator of a land cover's 'naturalness'. Thus the different land covers denote the habitats in this analysis. By itself, naturalness is not an ES but it is assumed that the more 'natural' a land cover is the more services and their provision are available (Mittelbach et al. 2003). This indicator is calculated by relating an 'integrity value' assigned to each land cover with the factors threatening its persistence. The integrity value corresponds to the potential of a land cover to serve as habitat for plant and animal species. Values ranged from 0 to 1, and were assigned by discussion and consensus with all project members. Factors threatening the species, the integrity and persistence of a particular land cover are all possible land use types that facilitate the appearance of negative impacts. These threats were considered to be roads and railways, mining, urban areas and settlements, and five different combinations of agricultural and pasture lands. As with the integrity value, each threat is related to each land cover by a value (between 0 and 1) describing how this threat would negatively affect that land cover (Supplementary Material III). For example, a road will have a stronger negative effect if it is built over a forest than over bare soil. At the same time, the maximum impact distance of each threat must be specified as well as the weights (relevance) of the threats in relation to each other. Because there was no reference as to the maximum impact distance of the selected threats, the model was run 10 times, each time with a different distance value, ranging from 1 km (one cell) to 10 km (10 cells). Weights were assigned by discussion and consensus with project members (Supplementary Material III).
The model used to calculate habitat quality is available in the analysis and mapping tool Integrated Valuation of Environmental Services and Trade-offs (InVEST) (Sharp et al. 2014). The CORINE level 2 categories of the land cover map for Colombia (IDEAM 2010) were used to describe the different land cover types in the study area. Official cartography provided by the CAR was used to generate the threat layers. All geoprocessing was done in the statistical software R (R Core Team 2014) and the geographical information system Quantum GIS (QGIS Development Team 2014).

Scenic beauty
In general, people prefer to be surrounded by natural (green) landscapes when travelling to rural areas for recreation and relaxation since these landscapes provide a notion of natural beauty (Daniel et al. 2012). Quantifying scenic beauty seeks to capture the quality of a certain place as a source of natural (green) views at the landscape scale that might be inhibited by human disturbances. As the basis for the scenic beauty calculations, a DEM for the study area was extracted from the Global Land Cover Facility (USGS 2008). The human disturbance elements were defined along with the maximum distance in the landscape at which these elements can be seen. The disturbance elements selected were roads and railways, urban and rural settlements, hydroelectric plants, electrical power lines, mining and quarry areas. The default distance of 8 km was used. All data were integrated in the 'Unobstructed Views: Scenic Quality Provision' module of InVEST (Sharp et al. 2014).

Water provision
The water balance was used as the indicator for measuring water provision. It was calculated as rainfall plus fog inputs minus actual evapotranspiration according to climate and remotely sensed vegetation (van Soesbergen & Mulligan 2014) as implemented in the WaterWorld model (Mulligan 2013). Details about model parameterization and algorithms can be found in Mulligan and Burke (2005). This model was used because field data is scarce and the model runs with globally available data.

Standardization
All indicators were rescaled to values between 0 and 10 by dividing the original values by the maximum value of each indicator and multiplying the result with 10. In this way, the difference between indicators with the smallest value at zero (no service at all) and indicators with the smallest value greater than zero remains, compared with rescaling all indicators to a minimum value of zero.

Multi-criteria analysis
A sMCA (Malczewski 2006) was carried out by overlaying all five indicators. In sMCA, weights (i.e., values between 0 and 1) are assigned to each criterion (each indicator in this case). As the ranking of criteria involves value judgements, this step often involves stakeholders. However, the diversity of opinions that exist in a stakeholder group, together with issues of representativeness, group hierarchy and power imbalance can lead to considerable uncertainty around criterion weights. We did not elicit this uncertainty from stakeholders in this study but simulated its effect by trialling possible stakeholder options of assigning weights. The weights for each indicator were calculated randomly from a uniform distribution and rescaled to sum to 1. In total, 1000 simulations were run, including the case where all indicators have the same weight and the five cases, where only one of the indicators is considered (weight = 1), that is, when stakeholders prefer only one indicator over all others. This number of simulations was enough to fill the value space when plotting all criteria on a pairwise basis. For each simulation, the weighted mean of all indicators was calculated.

Representativeness of cCPAs
From each simulation run, simulated Conservation and Protection Areas (sCPAs) were extracted by selecting the grid cells with the highest values until a total extension of 3679 km 2 was reached, the same areal extent as the current conservation and protection areas -cCPAs. The spatial congruence was calculated as the proportion of the area of the cCPAs that spatially intersects with the sCPAs.

Indicator correlation analysis
In order to evaluate the overall pair-wise correlation between all biodiversity and ESs indicators, Spearman's rank correlation coefficient was used.

Complementary areas
A consensus and an uncertainty map were calculated based on all simulations. The consensus map was determined by calculating the median of the 1000 simulations. The median was preferred over the mean because the distribution of simulation values was skewed and the median in this case fairly represented the central tendency of the data. We selected the central 90% confidence interval from the 1000 simulations as our uncertainty measure. The median and the confidence interval, leading to the consensus and the uncertainty map, respectively, were calculated on a cell-by-cell basis. Thus, a raster cell having a high median value (i.e., high consensus) and a low confidence interval value (i.e., low uncertainty) is a raster cell where all indicators consistently score high (i.e., high biodiversity and high provision of ESs) in the 1000 simulations. If the same cell would have not consistently scored high then it would have obtain a wider confidence interval and therefore high uncertainty. The areas of high consensus between indicators and low uncertainty were deemed the best multifunctional areas. A final ranking was established by making a scatter plot of the consensus and the uncertainty map values on a cell-by-cell basis. Four zones were identified in the scatter plot by selecting a threshold of the consensus values above which all areas represent high multifunctionality and a threshold for the uncertainty values above which the uncertainties are high. These thresholds were defined by the project members as the third quartile for the consensus values and the second quartile for the uncertainty values. The four resultant zones were called: (1) multifunctional areas: high consensus values and low uncertainty, (2) candidate areas: high consensus values and high uncertainty, (3) uncertain bad areas: low consensus values and high uncertainty, and (4) worst case areas: low consensus values and low uncertainty. The four areas were again mapped (Ligmann-Zielinska & Jankowski 2014) and overlain with the current conservation and protection areas to evaluate their degree of overlap, and to identify multifunctional areas without any current protection status.

Representativeness: overlap between current and simulated conservation and protected areas
The representativeness of cCPAs in protecting biodiversity and ESs is intermediate to low, depending on the combination of values assigned to the different criteria to define sCPAs (Figure 2), with values ranging from 3% to a maximum of 56% of overlap between sCPAs and cCPA. This indicates, first, that the selected priority areas for conservation and protection would vary significantly, depending on the relevance (weights) that stakeholders assign to the different criteria. Second, even under the most positive scenario, the overlap with the cCPAs is moderate (56%).
When considering each criterion independently, the overlap between sCPAs and cCPAs varies considerably. Were only ecological systems to be chosen, the priority areas would overlap at 52%. Were only water provision selected, the overlap would be 3%. The median value (48%) represents the case of choosing only habitat quality to select protected areas. Scenic beauty and sensitive species, when considered alone, have a low overlap with the cCPAs (30% and 20%, respectively) ( Figure 2).

Indicator correlation analysis
There is no evident geographical overlap when comparing the distribution patterns of the indicators (Supplementary Material IV). This is reflected in the medium to low global correlation coefficients obtained (Table 1). While high richness of sensitive species is found mostly throughout the high plateau around Bogotá, these areas display low or medium scores for all other indicators. The same situation occurs for areas of high water provision, where all other indicators have low or medium scores.
The highest positive correlation (r = 0.56, p < 0.001) occurred between ecological systems and habitat quality. High values of these indicators coincide mainly at high altitudes overlapping with the páramos ecosystem, while low values coincide throughout the whole study area. In fact, 79% of the study area has zero provision for ecological systems (as these are transformed land covers) (Supplementary Material IV), which corresponds to low values of habitat quality.

Complementary areas
Despite the weak global correlations between indicators, some areas of high value for all indicators can be recognized from the consensus map. These areas are largely distributed over the páramos ecosystems in the eastern sector of the study area with some patches in the north-central and north-western areas (Figure 3(a)). The lowest scores in the consensus map tend to be concentrated on the high plateau where the density of urban settlements is high. Some areas of high consensus scores show high uncertainty (Figure 3). By assuming a value of 5.5 (third quartile) to separate consensus scores and a value of 3.0 (second quartile) to separate highly uncertain areas, 26.3% of the selected high consensus score areas are highly uncertain (candidate areas: Figure 3(c)). Hence, 73.7% of the high consensus score areas can be considered effective multifunctional areas (3573 km 2 ).
Overlaying these multifunctional areas with the cCPAs indicates that 47.8% of the latter can be The map values can be plotted on a cell-by-cell basis (c) to define four areas based on thresholds assigned to the consensus (third quartile) and the uncertainty (second quartile) scores. These four areas represent regions with high consensus scores and low uncertainty (multifunctional areas), high consensus scores and high uncertainty (candidate areas), low consensus scores and high uncertainty (uncertain bad areas) and low consensus scores and low uncertainty (worst case areas). These four categories can be mapped again (d) and the areas overlapping with the cCPA or not can be identified. Figure 3(c) serves as a legend to Figure 3(d).
considered multifunctional areas. The remaining area of the cCPAs consists of 8.0% candidate areas, 28.3% uncertain bad areas and 15.9% worst case areas.
Of the multifunctional areas, 49.2% corresponds to cCPAs, the remaining 50.8% (1815 km 2 ) can be considered complementary areas. Some of these areas form small isolated patches spread throughout the study area, which are less easily connected to cCPAs. Other areas are contiguous to cCPAs, for example, in the Sumapaz area or next to the DMI Cuchilla de San Antonio in the Northwest (Figure 3  (d)). In the latter a continuous multifunctional area could be recognized as a potential complementary area to protect biodiversity and ESs simultaneously (Figure 3(d)).

Representativeness of protected areas
The current network of protected areas is only partly representing important areas for the protection of biodiversity and ESs. According to the simulations, representativeness of cCPAs for the conservation and protection of biodiversity and ESs will not exceed 56%, and it could be as low as 3% in the worst case. Focussing on areas with high consensus score and low uncertainty would lead to multifunctional areas that overlap with cCPAs at 49%. The sensitivity of outcomes to the choice of weights has also been found in similar studies applying multi-criteria approaches for land suitability analysis (Store & Kangas 2001;Chen et al. 2010;Ligmann-Zielinska & Jankowski 2014), highlighting once more the uncertainty in outcomes that would be generated by divergent stakeholders' preferences. However, through the consensus analysis and explicit integration of uncertainty in the criteria weights, we were able to identify those areas that, regardless of criteria weighting by stakeholders, are recurrently selected for showing high value and provision of biodiversity and ESs. This approach is rarely used within the multi-criteria context (e.g., Ligmann-Zielinska & Jankowski 2014), though Hengl et al. (2009), for example, applied regression-kriging to explicitly account for uncertainty in the prediction of species occurrence. Similar consensus approaches can be found in biodiversity prediction and climate modelling (Grenouillet et al. 2011;Heubes et al. 2013).
When analysing the geographical pattern of the limited overlap between cCPAs and multifunctional areas (Figure 3(d)), overlap is evident in the east of the study area where there is a high density of cCPAs. Although 18% of the study area has a protection status, there is a clear spatial bias towards the remaining fragments of natural vegetation in and around the high plateau area. This area corresponds mainly to the páramos ecosystem and thin stretches of Andean forests surrounding the high density of urban settlements on the high plateau. These ecosystems have been the focus of conservation efforts given the presence of priority species and ecosystems (i.e., páramos), and their assigned importance as buffer zones for areas of water supply and storage for Bogotá.
The indicator most congruent with cCPAs is 'ecological systems'. This was expected since areas of natural vegetation have been the main criterion in selecting protected areas. Noteworthy are the cCPAs that do not overlap with multifunctional areas. Again, these are mostly located in the east of the study area and they represent heavily transformed areas, albeit located within the páramos or wetland ecosystems. For example, one of the protected areas, the Reserva Forestal Protectora Cuenca Alta del Río Bogotá (Figure 1), recently lost more than half of its original area following the argument that most of its territory was currently used for cattle ranching and agricultural purposes. Another factor contributing to the amount non-multifunctional areas within cCPAs is the scale of analysis. The land cover in the eastern part of the study area is highly fragmented. Many of the fragments are part of a protected area, but when aggregating the original land cover data to a raster resolution of 1 km 2 many of the small fragments could be lost leading to low scores of service provision and hence low overlap with cCPAs.
Importantly, our results demonstrate that the most important areas for conserving the water provision service in the CAR jurisdiction, when considered as the most relevant indicator, lie outside the current network of protected areas. They may not have been considered so far as important areas for water protection because there is currently no significant demand placed on them. Indeed, conservation efforts to protect water supply and storage areas for Bogotá and other smaller settlements within the CAR jurisdiction have so far focused on the páramos and cloudy Andean forests that lie outside our study area.
Our results are comparable to studies of other biological groups, ESs and biogeographical zones. These studies all show that current networks of protected areas fail to protect extensive areas important for biodiversity and ESs (Cabello et al. 2013;Durán et al. 2013;Castro et al. 2015). Among the reasons for this trend is, as in this study, the lack of congruence between different ES and biodiversity indicators. Although the cCPAs have been established to fulfil certain conservation objectives, integrative approaches where multiple criteria are considered highlight those areas deserving more attention and protection efforts.

Congruence between indicators
We found no evident geographical congruence between biodiversity and ESs, with correlations ranging from moderate to low (positive and negative). The highest positive correlation between ecological systems and habitat quality (ρ = 0.56) highlights that the natural (non-disturbed) ecosystems in the study area are under high pressure from all the threats considered in this study. Other studies have shown that the amount of spatial co-occurrence between different ESs and between ESs and biodiversity is highly variable, and dependent on the region and scale of analysis, while in some cases spatial co-occurrence can be positive and significant for some biodiversity indicators (i.e., particular taxonomical groups) (Egoh et al. 2009). Schneiders et al. (2012), for example, found high congruence between ESs and biodiversity in an area 10 times smaller than our study area, while Turner et al. (2007) found an enormous variation at the global scale. At the national level, recent studies have shown a low overlap (< 1% of the country) of several ESs . The lack of spatial congruence found in our study highlights the danger of utilizing one or few indicators for delineating conservation areas, since protecting few will not secure protection for others (Zhang et al. 2015). Yet, as the number of indicators increases so does the value judgement of criteria weights, which is prone to stakeholder disagreement and power play. Hence, criteria weights are best deliberated openly and transparently, and remaining disagreement propagated as uncertainty through the analysis chain as in the present study.
Noteworthy is the negative correlation between sensitive species and all other indicators. In general, high richness values of sensitive species are distributed over the high plateau, that is, over highly disturbed areas. Orogenic phenomena, accompanied by processes of isolation and speciation have given rise to a high number of endemic species in this area (Fjeldså et al. 2012;Graham et al. 2014). In addition, the fact that the high plateau was a former lake means there are now several wetlands, which are used as temporary stops for many migratory species (SELVA 2012). However, the high plateau being an area with a high density of human settlements, many species have received a threatened status. Urban areas, in fact, constrain the provision of all ESs considered in this study. In contrast to our findings, other global and regional studies have shown that the spatial co-occurrence of ESs and biodiversity is positive (Mittelbach et al. 2003). However, studies at national and subnational levels have shown that the relation at those scales becomes weaker (Anderson et al. 2009;Armenteras et al. 2015). We have shown in this paper that the relation can become even negative at local scales.
The negative correlation between sites of high water provision and all other services is in line with other studies (Chan et al. 2006;Egoh et al. 2009). Areas of high water provision are located in the northwest of the study area. Although these areas are partially covered by remnants of tropical montane forest (TMF), whose important role in providing hydrological services is well known (Bruijnzeel 2001;Sáenz & Mulligan 2013), we might have overestimated this service given the data sources at hand. The data for the water balance equation used here comes from global datasets which might misrepresent local conditions, leading to over-or under-estimation of both the amount of water produced and the area providing this service.

Complementary areas
Despite the limited congruence between indicators, multifunctional areas, that is, areas of high biodiversity and provision of ESs were identified. Almost half of them intersect spatially with the current network of protected areas but the other half has no protection status and can thus be considered complementary. A similar situation was found in the dry region of Andalusia in the Iberian Peninsula where 59% of the total supply of regulating ESs occurred within protected areas (Castro et al. 2015). Some of the complementary areas act as buffer zones of key ecosystems and protected areas, for example, surrounding the Sumapaz region. Another important complementary area is located adjacent to the Distrito de Manejo Integrado (DMI) Cuchilla de San Antonio in the northwest of the study area (Figure 3 (d)). In general, the existence and distribution of complementary areas leads to the question whether protected areas could be expanded or new ones established. However, given the constraints that such a process would involve, we advocate, as in similar investigations (Palomo et al. 2014),for more practical actions that are feasible to implement and would account for the local communities living in these territories, like voluntary conservation agreements or payments for ESs schemes. Further analysis, incorporating socio-economic factors, such as land tenure, would have to be carried out in order to focus conservation efforts on places where significant benefits for human well-being can be achieved with least resource investment (Luck et al. 2012).

Methodological approach
The multi-criteria approach is one of the most common modelling techniques used in Colombia and abroad for land planning studies. We adapted the methodological workflow proposed by Ligmann-Zielinska and Jankowski (2014) that explicitly integrates spatial uncertainty with the consensus map (the 'average suitability score' map in their case). In the present study, we were interested in the identification of multifunctional areas ('robust areas' in their study), both to evaluate representativeness and define complementary areas. Ligmann-Zielinska and Jankowski (2014), in contrast, focused on candidate areas, that is, regions of high consensus score and high uncertainty in order to identify potential restoration areas for a key species.
Other variants of the multi-criteria approach have been used to evaluate spatially explicit trade-offs between ESs in order to identify priority areas for their conservation (Schröter & Remme 2016). For example, Zhang et al. (2015) and Liu et al. (2013) used ordered weighted averaging (OWA) multi-criteria evaluation, which evaluates different scenarios by iteratively re-ordering weights. Area-selection optimization algorithms, which are common in systematic biodiversity conservation, were used by Chan et al. (2006) to define priority areas for different ESs. They also evaluated the congruence between the selected areas. One of their challenges was to select the units and amount of representativeness of each service for the algorithm to be able to find priority areas. New optimization techniques have been developed to avoid this decision in the modelling process (Lehtomäki & Moilanen 2013).
The methodological approach used in this study and the results thereof can be considered a practical application of what the Colombian legislation labels as the Main Ecological Structure (MES) (Estructura Ecológica Principal) (Andrade et al. 2008). The MES is used to designate the network of areas for the longterm conservation of biodiversity, its functionality and the provision of ESs to secure the well-being of society (IDEAM 2011). Although practical examples at local scales in Colombia are still scarce, we encourage other scientist and decision makers to use our study as an example case to be replicated in land planning efforts. This could be promoted by current initiatives towards adapting regulatory frameworks to include the ESs approach (Albert et al. 2014), such as those going on in Colombia.
However, challenges remain in relation to data availability and quality and the selection of robust tools to model ecological processes. For example, all indicators in this study were assumed to accurately characterize the aspect they represent; uncertainties related to each of the models and their outputs were not explicitly considered. SDMs, for example, are sensitive to data quality (García Márquez et al. 2012), modelling technique (Dormann et al. 2008) and spatial autocorrelation (Dormann et al. 2007), among other things. Water provision, in turn, was estimated from a global database that might be inaccurate at local scales and habitat quality and scenic beauty were subject to expert opinion with potentially different outcomes when different stakeholders parameterize the model (Hamilton et al. 2015).

Conclusion
Current conservation and protection areas in the CAR jurisdiction have a moderate degree of biodiversity and ESs representativeness. A dense network of conservation areas exists on the high plateau and the surrounding hills. Attention should now be directed to the western region of the study area, not only because the number and extent of protected areas is low in this region, but also because important areas for biodiversity and the provision of ESs were found there. Management actions and land planning strategies for these areas require their consideration and acknowledgement as complex social-ecological systems and the factors that would make these areas resilient in the long term (Cumming Forthcoming 2016), as well as the harmonized work of environmental authorities and local communities, for example, through the implementation of conservation agreements or payments for ESs schemes, and not necessarily through the establishment of new protected areas.
This study and the indicators selected, considered as priority management criteria by the local environmental authority, served as a pilot to determine which areas are important for environmental conservation that currently do not have any protection status. New areas could emerge when adding new indicators to the analysis. We encourage the selection of those indicators in a joint effort with local stakeholders. Future research should concentrate on data collection and modelling of different ESs at the local scale. Our simulations highlight the sensitivity of protected area representativeness to stakeholder preferences of indicators. We thus encourage local land planning to consider participatory approaches in which stakeholder preferences are quantified, possibly harmonized (consensus) and any remaining disagreement propagated as uncertainty through the analysis chain; because these decisions will finally configure the location of areas for the long-term conservation of biodiversity and ESs.
Galina Churkina, Jane Coates, Leticia Lima and Marisa Beck for their insightful feedback on earlier versions of this manuscript. We would like to thank two anonymous reviewers for their comments.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by the Corporación Autónoma Regional de Cundinamarca -CAR; Conservation International -Colombia; German Excellence Initiative.