Prioritizing Spatially Aggregated Cost-Effective Sites in Natural Reserves to Mitigate Human-Induced Threats: A Case Study of the Qinghai Plateau, China

Abstract

Anthropogenic activities often lead to the degradation of valuable natural habitats. Many efforts have been taken to counteract this degradation process, including the mitigation of human-induced stressors. However, knowing-doing gaps exist in stakeholder’s decision-making of prioritizing sites to allocate limited resources in these mitigation activities in both spatially aggregated and cost-effective manner. In this study, we present a spatially explicit prioritization framework that integrates basic cost effectiveness analysis (CEA) and spatial clustering statistics. The advantages of the proposed framework lie in its straightforward logic and ease of implementation to assist stakeholders in the identification of threat mitigation actions that are both spatially clumped and cost-effective using innovative prioritization indicators. We compared the utility of three local autocorrelation-based clustering statistics, including local Moran’s I, Getis-Ord Gi*, and AMOEBA, in quantifying the spatial aggregation of identified sites under given budgets. It is our finding that the CEA method produced threat mitigation sites that are more cost-effective but are dispersed in space. Spatial clustering statistics could help identify spatially aggregated management sites with only minor loss in cost effectiveness. We concluded that integrating basic CEA with spatial clustering statistics provides stakeholders with straightforward and reliable information in prioritizing spatially clustered cost-effective actions for habitat threat mitigation.

1. Introduction

Expansion and intensification of anthropogenic events, especially those from land use activities, such as suburbanization [1] and urban sprawl [2], have led to the degradation and fragmentation of wildlife habitats, which have caused extensive loss of biodiversity regionally and globally [3]. Various ecological compensation efforts have been invested to mitigate these human-induced stresses and restore degraded land and habitats [4], including migration of local communities and conversion of associated land use legacies, such as residential land, farmland and local roads to natural lands, such as forest and grassland [5]. Investments on the alleviation of land use intensity and associated threats on valuable habitats are commonly recognized as a key and typically first-step process in systematic conservation planning [6,7]. However, knowing-doing gaps in stakeholder’s practice of prioritizing locations for these mitigation actions have lowered the efficiency of many conservation efforts [8,9]. This leads to a major concern among stakeholders—that is, where to allocate limited resources in an effective manner to actions aiming at eliminating anthropogenic impacts on wildlife habitats under various constraints.

A number of approaches, such as threat mapping, benefit maximizing, and cost minimizing have been traditionally used to address this prioritization concern [10]. The threat mapping approach focuses on habitats that are severely imperiled by stress factors but ignores the returns and costs associated with identified threat mitigation actions, which might result in inefficient and unrealistic actions [11]. Moreover, the efficacy of the benefit maximizing and cost minimizing approaches largely depends upon the intensity of spatial correlation between costs and benefits [12], which may lead to selecting uneconomic sites and omitting cost-effective ones. Comparatively, by using return on investment (ROI) analysis [13,14] to prioritize habitat restoration actions, stakeholders are informed of making transparent and logical decisions in a systematic manner in terms of where to allocate actions that could achieve the largest biodiversity returns with lower costs [15,16]. When using the ROI method, ranking schemes are often applied to compare alternative sites for actions, frequently based on their returns on per unit investment, which is easy to understand and implement [17].

A common implementation of ROI is the cost effectiveness analysis (CEA) [18,19], which identifies alternative, preferred actions, based on the integrated consideration of management benefits and corresponding costs. Critical steps in applying CEA to threat mitigation activities include evaluation of expected benefits, estimation of costs and setting of priorities for actions [14]. Typically, benefits are measured in nonmonetary units in the context of biodiversity or habitat conservation. Estimated cost could be expressed either in monetized or non-monetized units based on data availability, each having its advantages and disadvantages [8,12,20,21]. Priority settings for alternative actions can be realized by directly ranking schemes in order of benefit in per unit cost [17,22], or using more complicated single- or multi-objective optimization algorithms [23]. The utility of CEA in identifying cost-effective conservation and restoration actions has been demonstrated by a number of studies [19,24,25]. Very sophisticated CEA framework could be designed by incorporating additional ecological considerations (such as multiple ecosystem services, landscape connectivity, and policy incentives) into benefits estimation [26,27,28], utilizing complex economic models in cost assessment [12], and integrating sophisticated optimization algorithms in priority setting [29]. Although providing substantial informative guidance, these complicated CEA approaches are not extensively being used by stakeholders due to the perceived difficulties in comprehension, implementation and interpretation [22], which reveals a common gap in ecological conservation and restoration practice [30].

Further, actions identified by CEA could be dispersed in space due to the lack of considering inherent information on spatially correlated structures between prioritized sites [24]. Generally, stakeholders tend to focus their management efforts more in spatially clumped locations rather than dispersed ones [31,32,33], which is a basic principle of systematic conservation planning [7,34]. Because of the spillover or halo effect in biodiversity between adjacent habitat patches [35,36], the summation of spatially dispersed threat management actions is not the same as spatially clumped actions with respect to the restoration of degraded habitats [33,37]. Gathering actions across space contributes more in improving habitat connectivity and compactness [38,39], and thus can solidify energy and material flow among neighboring habitat patches, due to corridor effects on the landscape [35]. Through the removal of barricades in habitat connectivity pathway, spatially aggregated actions of threat mitigation are ensured to be more efficient in repairing the damaged patches and links in habitat networks [40,41,42]. This additional management effectiveness from clustered actions tends to increase along with the intensification of spatial aggregation among identified actions [43]. In addition, clustering threat management efforts in space could help reduce costs [33]. For example, cooperation and interaction between neighboring actions could help reduce labor, management and transportation costs due to economies of scale [43]. Moreover, clumped habitat restoration actions facilitate positive endogenous spatial interactions between neighboring landowners, thus leading to increase in social acceptability and decrease in costs when official agencies seek to purchase conservation easements for these habitat patches [44]. Apparently, spatial autocorrelation that reveals local association or similarity in neighboring sites is an important concern in spatial prioritization of habitat conservation and compensation actions [7]. Spatially correlated conservation and compensation actions are more efficient in promoting the functionalities of habitat networks through the enhancement of compactness, connectivity, and contiguity among habitat patches and corridors, thus improving the elasticity of landscapes and strengthening the interactions across habitat mosaics [45]. Integrating spatial autocorrelation analysis with CEA has the potential of allowing stakeholders to identify cost-effective threat management actions that are also aggregated across space. Although promoting the cost-effectiveness of threat management actions, there is less evidence that spatially correlated structures among actions are being extensively taken into consideration in stakeholder’s practice of prioritizing locations to invest limited resources when using ROI [22,46]. This knowing-doing gap is in part due to the lack of a straightforward and easily applied approach for stakeholders to integrate spatial autocorrelation information among actions into the prioritization indicators [9,24].

Heuristic optimization/search algorithm (such as those in the Marxan and Zonation conservation software) are a valuable tool that allows for the exploration of spatial autocorrelation between identified sites when prioritizing conservation and restoration actions [47]. Many pieces of research have demonstrated the utility of these algorithms in promoting compactness, connectivity, and contiguity among locations [48,49,50]. The disadvantage of these algorithms lies in the opaque “black-box” implementation and the complicated algorithms that require a number of inputs and hyperparameters to be set, which usually requires additional expertise and induces additional hindrance in comprehension and interpretation [24,47]. Besides, heuristic optimization algorithms often require multiple runs to get robust outputs, which often calls for intensive computation resources and time in large-scale and fine-grained analysis [51]. In addition, the differences between these multiple outcomes often induce extra difficulties for stakeholders to interpret the results in a defensible way. These disadvantages can lower the acceptability for stakeholders to perform these complicated tools in day-to-day decision-making on habitat restoration actions, and make them remain theoretical in the academic arena. Furthermore, complicated algorithms do not always guarantee better solutions for prioritization of habitat conservation and restoration actions [52].

Spatial statistics, such as local Moran’s I [53] and Getis-Ord Gi* [54] are alternative techniques that are easy-to-use, straightforward, and rigorous, providing strong support for evaluating spatially correlated structures in location-based data. These local clustering approaches have been used to detect spatial autocorrelation characteristics of species, ecosystem services, valuable habitats, thus guiding stakeholders to focus their efforts in space [51]. For example, Getis-Ord Gi* was used to identify social preferences and hot spots of ecosystem services in the Southern Rocky Mountains to promote efficient management of service-based resources [55]. Another example is that clustering types of habitat quality in Italy were analyzed using local Moran’s I statistic to assess the relationship between the distribution of natural habitats and protected areas [56].

In this study, we explore a straightforward and simplified integration of spatial autocorrelation information into CEA using these local clustering statistics. We adopt a basic CEA implementation to integrate the expected benefits and associated costs in threat mitigation actions. The expected contribution of this study is to promote a quick, transparent, and easy-to-use prioritization framework for stakeholder’s decision-making on the systematic mitigation of human-induced threats in degraded land and habitats. Our approach is expected to facilitate the accessibility of ROI for stakeholders to support their decision-making on habitat management actions by showing how just a substitute of the prioritization indicator could help allocate limited resources to actions that are both cost-effective and clustered in space through a straightforward but rigorous way. Therefore, the objective in this study is to investigate the feasibility of integrating a basic CEA and spatial clustering statistics to prioritize spatially aggregated and cost-effective sites for threat mitigation actions. We present a spatially explicit analysis framework to address the following research questions: (1) what is the utility of the CEA method in identifying cost-effective threat management actions, compared with traditional methods that only focus on biodiversity benefits, threats, or costs? (2) how does the integration of spatial clustering statistics with CEA help identify spatially aggregated cost-effective sites? Also, spatial correlation structures in adjacent actions detected by different local clustering statistics could be different [57]. So, it is necessary to investigate how different clustering methods impact the prioritization outcome for threat management actions to shed lights on the selection of appropriate statistical approaches, which will provide more insightful decision-support for stakeholders. We illustrate the prioritization framework using a case study from natural reserves (NRs) in the Qinghai Plateau, China, where many habitats with high-biodiversity value for endemic species are threatened by land use activities from local communities.

2. Materials and Methods

2.1. Study Area

Qinghai Plateau, as a core region of the Tibetan Plateau, is located at the west of China with a total area of 722,300 km², and a population of 5.88 million in 2015. There are three cities and more than 4281 administrative villages distributed in this region. The average elevation is about 3000 m and the highest peak reaches 6860 m. The Yellow River, Yangtze River, and Lancang River originate in Qinghai Province and take about 49.20%, 25%, and 15% of their total amount of water flow respectively, influencing almost 40% of the world population [58]. Qinghai is recognized as one of the world’s highest altitude area with the most abundant biodiversity and unique natural habitats for many endemic species [59]. To protect the unique biodiversity and fragile ecosystems in this region, Qinghai Province has established 11 NRs (Figure 1) with a total area of 221,200 km², which accounts for 30.6% of land in Qinghai and 14.31% of the total area of NRs in China. Each NR can be divided into three function zones: Core zone, buffer zone, and experimental zone [60]. The core zone is designed to completely isolate the most suitable habitats for rare and endangered species from human activities. The experimental zone is characterized by suitable habitats subject to human influence, and the buffer zone serves as a barrier to alleviate impacts exerted by human development in the experimental zone on valuable habitats in the core zone [61]. Some large NRs in the study area comprises of several small protected areas (PAs) that are separated in space. Thus, the total number of PAs is up to 37 and each PA is administrated by an affiliated government agency.

While human activities are legally prohibited in these PAs, many of the PAs report different levels of human-induced habitat degradation and biodiversity loss [62,63]. According to the official land use survey data in 2015, there are more than 7973 spontaneous villages scattered within and surrounding these PAs (see Figure 1). Although the intensity of the human influence in these PAs is reported to be relatively low (compared with PAs in urbanized regions), it keeps rising apparently in recent decades [64]. Additionally, the extremely fragile natural ecosystem in this region amplifies the human-induced impacts on habitats and local biodiversity [65]. This region is heavily engaged with many ecological restoration projects funded by international organizations and Chinese governments to restore degraded natural habitats and ecosystems, as well as native biodiversity. For example, 23.86 million US dollars have been invested in 2013 by the UNDP-GEF (United Nation Development Programme-Global Environmental Finance) [66] and China’s central government to the Sanjiangyuan NR to recover habitat and conserve biodiversity [67]. More efforts are planned to further improve the habitat quality and biodiversity richness in this region, including the establishment of China’s first and largest national park [68]. Reducing or eliminating anthropogenic activities that degrade valuable habitats in these PAs is a key concern in these projects [62]. However, a lack of transparent and easy-to-use tools in supporting the planning and allocation of these projects has led to injudicious utilization of limited resources. Therefore, it is necessary to systematically arrange these actions so as to get better use of the limited resources through spatially prioritizing locations for threat mitigation.

2.2. Materials

Data used in this study include natural reserve extent, and land use/cover map. The spatial boundaries of PAs, including the extent of function zones, were provided by the Environmental Protection Bureau of Qinghai (http://www.qhepb.gov.cn). A lattice of hexagonal units with a length of 1km is imposed over the study area for fine grained analysis (see Figure S1 in Supplementary Materials). The total number of hexagons in PAs is 89,246, which are regarded as spatial units of the analyses in this study. Land use/cover data of Qinghai Province in 2015 was obtained from the Land Consolidation and Acquisition Center of Land and Resources Department of Qinghai Province (http://www.qhlr.gov.cn/tcnf?vid=2300), which contains the land use/cover type, boundaries of administrative districts, distribution of roads, rural settlements, and mining sites (see Figure S2 in Supplementary Materials). The land use/cover data were officially produced by the county-level department of land use survey in Qinghai province. Visual inspection and field investigation methods were combined to classify different land use/cover types. First, Digital Orthophoto Maps (DOMs) were produced using high-resolution aerial imagery and SPOT 5 satellite imagery. Second, geographic objects in these DOMs were visually identified and vectorized based on the designed interpretation keys in the study area. If the land type of a geographic entity could not be classified through a visual inspection from the DOMs, field investigation will be conducted. Based on the resulting outputs from visual inspection and field investigation, intensive field validation was launched to refine and verify the classified maps to generate the final land use/cover maps. The field validation was implemented by comparing the observed and identified land categories of a certain number (more than 20%) of randomly selected objects in the classified maps. The validation accuracy of the final land use/cover maps should larger than 93%, otherwise visual inspection and field investigation should be implemented again [69,70]. The original land use/cover maps were produced in 2010, and land change survey was launched annually to keep these data up-to-date. All the spatial datasets were projected to the WGS-84 coordinate system and rasterized with a spatial resolution of 30 m by 30 m.

2.3. Methods

In this section, we present the spatial analysis framework that integrates a basic CEA and spatial clustering statistics for the prioritization of spatially aggregated cost-effective sites to mitigate human-induced threats on habitats. Figure 2 shows the overall framework design. First, expected biodiversity benefits from the mitigation of human-induced threats are estimated using a biodiversity model. Moreover, we estimate the associated costs in each site using a non-monetized proxy. Then, cost effectiveness (CE) ratio is calculated through the CEA method. Spatial statistic methods are utilized to measure spatial aggregation of management units. Last, we integrate the CE ratio and spatial clustering statistics to identify sites that are both cost-effective and aggregated in space.

2.3.1. Cost Effectiveness Analysis

• Evaluation of expected biodiversity benefits

We define the biodiversity benefits as the improvement of habitat’s capability to support all levels of biodiversity. Specifically, we calculated the expected biodiversity benefit in a hexagon unit by summing up the increase in habitat quality (HQ) and decrease in habitat degradation (HD) after removing specific human land uses within or surrounding this unit. First, land cover types that are suitable habitats for all levels of biodiversity, and anthropogenic land uses that are threats to those habitats were identified from land use/cover maps. Then, an adjusted land use/cover map was created by removing three human land uses that are farmland, rural residential land, and local roads. Last, the HQ and HD indices for the scenario-based and the original land use/cover map were estimated to quantify the biodiversity benefits in each unit. The HQ and HD were calculated using the biodiversity model (available as a toolbox in the ArcGIS software [71]) provided by the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) modeling framework (https://www.naturalcapitalproject.org/invest/) based on four factors: The relative suitability of each type of habitat to sustain biodiversity, the relative impact intensity of each threat to its surrounding habitats, the relative sensitivity of each habitat type to each threat, and the distance of habitats to the source of threats [56,72,73,74]. The detailed description of how the expected biodiversity benefits were estimated are presented in Appendix A in Supplementary Materials.

• Estimation of threat management costs

To model the distribution of cost related to each unit, we first count the number of cells that have one of the three human land uses in specified buffers of a hexagon unit (using the buffer analysis function in the ArcGIS software). The buffer distances are the same as the maximum influence distance of the specific threat factors. Second, a relative price for the removal of each of the three human land uses is estimated. We assume that the final price for removing a threat cell is only determined by the type of threat factor inside it. That is, the removal price for cells that have the same human land use are identical across space. Based on the former two steps, the estimated cost for removing human land uses related to a hexagon unit could be calculated as follow.

(1)$C_i={\displaystyle \sum }_{j=1}^n{P}_j\times N_j,$

(2)${\displaystyle \sum }_{j=1}^n{P}_j=1,$

where $C_i$ represents the estimated cost of hexagon unit i. $P_j$ is the standardized relative price for removing a cell with threat factor j inside it. $N_j$ is the total number of threat cells that impact habitats in unit i. j stands for the type of human land use that can be removed, including farmland, rural residential land, and local roads.

The relative price of each type of land use indicates the ratio between the final monetized prices for removing the threat cell. As the threat of intense human land use on the surrounding habitats is partially relevant to the associated financial investment [75,76], it is intuitive to prorate the relative management price for each of the three human land uses to their threating intensity W_r. In this case, the relative removal price for farmland, rural residential land, and local roads are set as 6/19, 8/19 and 5/19. It should be noted that the two sets of relative prices differ only if the relative magnitude between prices differs. For example, relative price of 0.2, 0.4, 0.8 is the same as 0.1, 0.2, and 0.4 after standardization.

• Calculation of cost effectiveness ratio

The CE ratio generated through CEA indicates the management benefits in per-unit cost. We calculate the CE ratio for each hexagon unit using the following formula.

(3)$C{E}_i=\frac{B_i}{C_i},$

where $C{E}_i$ is the CE ratio of unit i. B_i and C_i stand for the expected biodiversity benefit and estimated cost for unit i.

2.3.2. Identifying Spatially Clustered Sites for Threat Management

The spatial aggregation feature of CE ratios in units could be analyzed using autocorrelation-based spatial statistical approaches. We chose to use three local spatial clustering approaches that are available in the ArcGIS software: The Getis-Ord Gi* statistic [54], A Multidirectional Optimum Ecotope-Based Algorithm (AMOEBA) [77], and local Moran’s I [53].

The Gi* statistic analyzes the values of each hexagon unit and its neighboring units within a certain distance. It calculates the Gi* statistics as the summation of the product between the spatial weights and CE values in neighboring units divided by the summation of CE values of all units across the study area. The z-score and p value were also reported to test the statistical significance of the Gi* statistics. The z-score can be used to identify spatially clustered units with higher or lower CE ratios. Based on the Getis-Ord Gi* statistic, AMOEBA was applied to identify clustered sites for threat management. While the Getis-Ord Gi* statistic is applicable for detecting regular spatial clusters, AMOEBA is suitable for identifying clusters with irregular spatial shapes. Similar to the Getis-Ord Gi* statistic method, a statistically significant higher positive Gi* value from AMOEBA for a unit indicates inclusion into a spatial cluster with relatively higher CE ratio, and vice versa.

The local Moran’s I statistic identifies spatial similarity and dissimilarity by grouping each unit into a cluster/outlier type. The local Moran’s I statistic is calculated through multiplying the subtraction of the CE value in the focal unit from the average CE value in the study area by the sum of the product between the spatial weights and the subtraction of the CE values in the neighboring units from the average CE value in the study area. Then, the multiplication result was divided by the summation of the subtraction of the CE values into all units from the average CE value in the study area. Similarly, the z-score and p value were computed to test the statistical significance.

The mathematic algorithms of Getis-Ord Gi* and local Moran’s I are similar in general. The major difference between the two methods is how the relative magnitude of value in the focal unit and values in its neighboring units is integrated. When detecting spatial clusters using Getis-Ord Gi* statistic, a unit with relatively low value could be incorporated into a statistically significant hotspot if its neighboring units have higher values. While using local Moran’s I statistic, this unit will be identified as part of a low-high outlier instead of being incorporated into a hotspot (statistically significant high-high cluster). In this case, hotspots detected by Getis-Ord Gi* might include some units with relatively low values, while all units in the local Moran’s I hotspots are guaranteed to have high values. The detailed descriptions of the three autocorrelation-based spatial statistical methods were given in Appendix B in Supplementary Materials.

2.3.3. Priority Setting

Cost effectiveness is widely used as an efficient indicator to select cost-effective management units under given budgets. Stakeholders usually do this by ranking units in order of their CE ratio for the sake of perceptual intuition and execution simplicity [30,78]. To select sites that are both spatially clustered and generate high benefits under given budgets, we rank units in order of their spatial clustering statistics to conduct the prioritization process. Specifically, the prioritization process could be conducted by ranking all the hexagon units from hotpots to outliers, and then to cold spots in descending order of their z-scores when using the Getis-Ord Gi* statistic as the priority indicator [51]. Similarly, the prioritization process based on AMOEBA statistic could be achieved by ranking all the units from high-value clusters to outlier and then low-value clusters in order of the maximized Gi* statistics. As for the local Moran’s I statistic, the prioritization process is launched by ranking all the units in order of high-high hotspots (a high-value unit is surrounded by neighbors with high values), low-high outliers, insignificant outliers, high-low outliers, and low-low cold spots. Units with the same cluster/outlier type are ranked in order of their z-scores.

A key aspect for the calculation of these statistics is the definition of neighborhood for a focal unit i. Contiguity and distance-based neighborhood definitions are two commonly used methods to specify neighboring units. Different neighborhood definition strategy could make differences in the calculation of spatial clustering statistics. We first used the contiguity neighborhood in all the clustering statistics to compare the feasibility of these statistics in identifying spatially aggregated cost-effective units. Then, we investigated how different neighborhood distances impact the accrued benefits and spatial aggregation of selected units under given budgets. Specifically, this investigation starts from a neighborhood distance of 2.7 km. Then, we increase the neighborhood distance with a 1 km increment.

2.3.4. Comparison with Traditional Prioritization Index

To verify the utility of CEA in identifying cost-effective actions, we ranked all the units in descending order of their expected benefits so as to compare the cost effectiveness of the CEA approach with that of selected units when management goal is defined as maximizing benefits without simultaneous consideration of the estimated costs. Similarly, to examine the effectiveness of selected sites when management goal is to minimize costs without any consideration of the biodiversity benefits, we ranked all units in ascending order of their estimated costs. Units were also ranked by their habitat degradation to identify units that have the highest human-induced threats. Moreover, to clarify whether the performance of the CEA method is better than an arbitrary selection situation, we compared the cost effectiveness of randomly selected units with units selected by CEA.

Curves of accrued expected benefits against cumulative estimated costs for all the above prioritization indices are plotted to evaluate how biodiversity benefits are expected to be attained from a given investment budget. A budget is defined as the percentage of investment to the total investment which is the cumulative estimated costs needed to remove all the specific human land uses in PAs. Moreover, to investigate the spatial aggregation features of selected units, the number and mean area of landscape patches that consist of these selected units are calculated as well. Hotspots detected by spatial clustering methods are often regarded as favorable sites for threat mitigation actions. Therefore, we also compared the cost effectiveness and spatial distribution of threat mitigation hotspots identified by different spatial clustering statistics.

3. Results

3.1. Cost Effectiveness of Units Prioritized by CEA and Traditional Index

Compared with traditional prioritization indicators, the accrued expected benefits from threat mitigation actions are the highest under given budgets when units are prioritized by CEA (as shown in Figure 3). In contrast, the accrued benefits in units that are randomly selected present a nearly linear relationship to their accrued estimated costs. The cost effectiveness curves are much similar when units are ranked in order of their expected benefits or degree of habitat degradation. Notably, the accrued benefits are even smaller than those of the random situation when management goal is to minimize costs. The CEA approach can gain the highest benefits at small initial investment, then benefit returns on given investment continue to decrease as more subsequent investments are placed. Comparatively, the biodiversity benefits are extremely low under the same initial budget when using cost as the ranking indicator. For example, when 10% investment is budgeted, the total benefits of units selected by CEA reaches 25.66%, which is 4.43 times higher than that of sites selected by minimizing costs, and 2.63 times larger than that of the random selection situation. Moreover, about half of the biodiversity benefits could be achieved with only 30% of the total investment when CEA method is used.

The distribution of management priorities for units ranked according to CEA and traditional indices varies significantly in space (see Figure 4). The total number of units with non-zero benefits is 21,984. As observed, the spatial distribution of high-priority units identified by CEA are more dispersed in space than units identified by the other three traditional indices. The layout of prioritized units sorted by benefits and threats are coincident, which are more concentrated and contiguous across space (especially in the southern and eastern regions). Apparently, ranking units by their estimated costs will make the distribution of high-priority units almost opposite to the layout of units created by using benefits and threats. In fact, most of the units gain high priorities when ranked by their estimated costs.

3.2. Spatial Aggregation and Cost Effectiveness of Units Prioritized by Spatial Clustering Statistics

Table 1 compares the cost effectiveness of threat mitigation hotspots identified by different spatial clustering statistics. It can be seen that AMOEBA presents the largest area of hotspots, including 2254 units, almost twice of the area of hotspots identified by local Moran’s I. The Getis-Ord Gi* hotspots give the largest benefits (28.50%), as well as the highest cumulative costs (13.72%). However, the local Moran’s I hotspots present the highest CE ratio (248.76%). Comparatively, the CE ratio of units selected by the CEA method under the same cumulative costs with local Moran’s I hotspots (8.96%) is 262.67%, which is moderately higher than the CE ratio of local Moran’s I hotspots.

Figure 5 shows the spatial distribution of threat mitigation hotspots from the three spatial clustering statistics. It reveals that the location of units identified by CEA has a high overall consistency with the layout of AMOEBA hotspots, and is different from hotspots identified by Getis-Ord Gi* and local Moran’s I. Hotspots from the local Moran’s I statistic is completely overlapped by the Getis-Ord Gi* hotspots. Obviously, hotspots detected by Getis-Ord Gi* statistic are more concentrated in space, which can be illustrated by the number and mean area of landscape patches in these hotspots (see Table 1). The total number of landscape patches detected using Getis-Ord Gi* is 336 and the mean area is 16.61km². The number of landscape patches of hotspots identified from AMOEBA is close to that from the CEA method, but is much more than that from Getis-Ord Gi* and local Moran’s I. The mean area of these hotspot patches from AMOEBA is much smaller. Compared with local Moran’s I, the Getis-Ord Gi* hotspots include more units that carry larger biodiversity benefits, and these units are obviously more clustered in space. But the CE ratio of these hotspots is a little smaller than that from local Moran’s I.

The accrued benefits are quite different when units are sorted in order of different spatial clustering statistics (see Figure 6). The total biodiversity benefits identified by these spatial statistics under given budgets are smaller than those of units identified by CEA, but the differences are marginal under small initial budgets. Similar to results from CEA, all the spatial statistics generate their highest benefits with small initial investment (about 10%). After that (higher than 10%), the return on investment starts to diminish. In particular, the local Moran’s I statistic generates almost the same benefits with the CEA method when 10% investment is applied. The biodiversity benefits based on the Getis-Ord Gi* and AMOEBA statistics are both larger than those from the local Moran’s I method. Among the three spatial clustering statistics, accrued benefits from the local Moran’s I statistic show the most rapid decrease, illustrated by its flatter curve for accrued benefits and costs. In fact, the return on investment for the local Moran’s I statistic is more similar to the random selection results when larger investment is budgeted. Comparatively, the shape of cost effectiveness relationship from Getis-Ord Gi* is more similar to that from the CEA method.

As the expected benefits of units are much smaller when units are ranked by the local Moran’s I statistic, we only focus on the spatial distribution of units ranked according to their Getis-Ord Gi* and AMOEBA statistics. Figure 7 indicates that the spatial pattern of the prioritization results from the two spatial clustering statistics are generally consistent and both are similar to those from the CEA method. To further evaluate the spatial aggregation of the selected units, the number and mean area of landscape patches identified from these two clustering methods under given budgets are calculated as well (see Figure 8). Results show that units are more spatially clustered if prioritized by their Getis-Ord Gi* statistic. The number of landscape patches from the Getis-Ord Gi* method is much less than that of the AMOEBA and CEA when cumulative investment is less than 30%. After that (higher than 30%), the gaps in patch number narrows as the cumulative investment grows. The same trend can be observed in the mean area of landscape patches.

Figure 9 shows how accrued benefits, number and mean area of landscape patches change along with neighborhood distance when using Getis-Ord Gi* statistic as the prioritization indicator. As can be seen, with the growth of neighborhood distance, the accumulated benefits decrease substantially at the beginning and then tend to be stable (Figure 9a). The number of landscape patches also shows a similar pattern of change (Figure 9b). The mean area of landscape patches generally keeps growing as an increase in neighborhood distance (Figure 9c). Notably, when neighborhood distance increases from 2700 m to 3700 m, the accrued benefits, the number and mean area of landscape patches show relatively more significant changes. Specifically, when budgeting 30% and using 3700 m as the neighborhood distance in the Getis-Ord Gi* statistic, the accrued benefits of selected units is 41.91% (14.92% less than CEA’s return (56.83%) on the same investment). But the mean area of landscape patches is 19.09 km² larger than that of the CEA results and the number of landscape patches is 943 less than the CEA units.

3.3. Actual Costs of Units Identified by Spatial Clustering Statistics

As mentioned above, the Getis-Ord Gi* method identified the most aggregated units with a small initial investment, and the total benefits of these units are high as well. Therefore we used the Getis-Ord Gi* statistic to select units under a small initial budget of 10% to illustrate how clustered units could help reduce actual costs and where to remove human land use activities. Figure 10 shows the spatial distribution of the human land uses associated with units selected by Getis-Ord Gi* and CEA according to the maximum influence distance of each land use type. It can be seen that the amount of human land uses identified by Getis-Ord Gi* is much less than that identified by CEA method, which means that the actual costs from Getis-Ord Gi* to achieve nearly the same amount of biodiversity benefits is much smaller. Specifically, when using Getis-Ord Gi* statistic, the total area of rural residential land and farmland that should be removed are 26.95 km² and 171.70 km². The length of the identified local road is 2328.73 km, which is much less than that of the CEA method. Furthermore, human land uses identified by Getis-Ord Gi* are more clustered in space.

4. Discussion

Limitations from financial constraints for habitat compensation actions highlight the necessity for stakeholders to allocate limited management resources to the most profitable locations that require human interventions. Due to the spillover effects in biodiversity, human-induced threats that impact biodiversity habitats in a region tend to degrade its adjacent areas. So, habitat management actions that are clustered in space will be more conducive to reap additional biodiversity benefits (such as bolstering landscape connectivity and contiguity) besides relieving human-induced threats. Further, actions located in spatially dispersed units are often more costly than those that are clumped in space, due to the economics of scale [33]. However, the spatial autocorrelation information among management units is often discarded in both the traditional methods and the CEA approach when prioritizing threat mitigation actions. Thus, units identified could be scattered across space, which will hinder stakeholders from obtaining additional benefits and saving costs due to connected spatial structures between actions [38]. This knowing-doing gap is in part due to the lack of a transparent and easy-to-use approach for stakeholders to incorporate the spatially correlated structure among actions into the prioritization indicators. In this study, a basic CEA and straightforward clustering statistics are coupled to achieve the goal of allocating management resources to units that are both cost-effective and clustered in space. The implementation of the promoted prioritization framework requires only basic knowledge on operating commonly used software (typically the ArcGIS and Microsoft Excel software) which would practically encourage the utilization of the ROI approach by stakeholders in decision-making on the mitigation of human-induced threats in degraded habitats.

4.1. Feasibility of the CEA Method

Previous conservation studies defined expected biodiversity benefits as the number of species that benefit from the mitigation of relevant threats. These studies were often conducted at coarse resolutions [20,79]. Species richness data have advantages in the representation of local biodiversity. However, presence and absence data of local species in a fine-grained analysis situation are typically unavailable [80], especially in large spatial extent analysis. Comparatively, fine resolution data on habitat conditions are easier to access and are often regarded as an efficient proxy to demonstrate the capability of a landscape to sustain all levels of biodiversity [25,81]. Researchers have suggested that taking immediate actions with readily available habitat data is sometimes more efficient than delayed actions waiting for high-quality species data [82].

HQ and HD are arguably appropriate indicators for habitat conditions [83]. HQ outlines the “goodness” of a landscape in terms of its capability and tendency of sustaining all levels of biodiversity. Areas with higher HQ have greater possibilities and better natural conditions to accommodate more native communities that anticipate richer biodiversity [84]. HQ has been used as an indicator for local biodiversity in some restoration projects [85,86]. Habitat degradation (HD) caused by human induced threats and other processes leads to the reduction in habitat extent and quality from its original status, resulting in a decrease in biodiversity persistence and resilience [87]. Spatial mapping of these threats, which reveals the spatial distribution of degradation, has been used as feasible indicators for the prioritization of biodiversity and habitat management efforts [88].

HQ and HD could be evaluated from several aspects [89,90]. For example, it could be estimated by comparing field survey data on habitat conditions to ex-ante standardized criteria or by using calibrated models to estimate HQ and HD from the perspective of the spatial relationship among species, threats, and ecosystem types. Furthermore, extracting information on habitat conditions from remote sensing data and comparing them with a historical reference status is also a commonly used approach. While from the viewpoint of supporting conservation management and facilitating stakeholder’s day-to-day decision-making, a measurement that could be used for quick assessment and comprehension is more preferred. Moreover, the calculation procedure should be replicable and transferable at multiple scales. The results are straightforward and can be easily interpreted and applied for temporal and spatial comparison [83]. Thus, the HQ and HD indexes generated by the InVEST biodiversity model are suitable and applicable [72]. This model maps degradation and quality of biodiversity habitats in relation to different types of land covers or ecosystems, as well as human-induced threats posed by surrounding anthropogenic land uses. The model runs based on the assumption that all landscapes are effectively protected and all human-induced threats to habitats are additional. The only dataset required by this model is land use/cover maps, which are more available than species data.

4.2. Spatial Clustering Statistics Identify Spatially Aggregated Cost-Effective Sites

Units ranked in order of benefits are continuous in space (see Figure 4). However, when costs are incorporated, high-priority units identified by CEA could become dispersed across space. This is a side effect that happens when strong positive spatial correlation exists between benefits and costs [91,92]. Our results demonstrate that spatial clustering methods that evaluate local association structure among neighboring units could help identify spatially clumped units with higher cost efficiency. Especially, the Getis-Ord Gi* statistic prioritizes more clustered units with only minor loss in the cost effectiveness when contiguity neighborhood was used. As neighborhood distance increases, the accrued biodiversity benefits drop while the intensity of spatial aggregation is improved. Compared with larger neighborhood distance, the contiguity neighborhood presents higher biodiversity benefits under given investment and the selected units are more clustered in space than the CEA results (see Figure 9). Stakeholders have to trade-off between different neighborhood strategies according to their expectation on cost effectiveness and spatial aggregation.

Different insights into conservation practice may be generated when using various spatial clustering statistics. Therefore, it is key to choose spatial clustering methods to be applied in the spatially explicit prioritization framework. When very limited investments are budgeted, the local Moran’s I statistic is suggested for this study. Hotspots detected by local Moran’s I statistic take the least investment and exhibit the highest cost efficiency, which is the closest to the cost effectiveness of units selected by the CEA approach. Although the spatial aggregation intensity of selected units is moderately smaller than the Getis-Ord Gi* method, it is still much better than that of units identified by CEA and AMOEBA methods. However, if spatial aggregation is of more interest, the Getis-Ord Gi* is a better choice. While the cost effectiveness of units selected by Getis-Ord Gi* is a little lower than CEA, the most aggregated clusters could be detected. When enough investment is budgeted, more units that are not included in statistically significant hotspots will be selected. Thus, ranking units in order of their Getis-Ord Gi* statistics is a more favorable option, because the accrued cost effectiveness of units ranked by the local Moran’s I statistic is much lower than that of units ranked by their Getis-Ord Gi* statistics (see Figure 6). As for the AMOEBA method, though the accrued biodiversity benefits in units selected by this method under the same investment are close to or even higher than those by the Getis-Ord Gi* statistic, units selected by Getis-Ord Gi* statistic are more clustered in space.

The amount of human land uses identified to be mitigated using spatial clustering methods under a small initial budget is less than that of the CEA method. This shows that spatially clustered units have potential in reducing actual costs of threat management actions due to the halo effect of threats in adjacent units. Threats that degrade habitats in a single unit will influence its neighboring units. This halo effect in biodiversity threats is often overlooked by stakeholders when threat intervention actions are allocated [32,93]. Spatial clustering methods based on autocorrelation give key consideration on spatially dependent structure between neighboring units, thus functioning well in identifying these spatially clustered cost-effective units. In addition, if threat management units are clustered, the identified human-induced threats relevant to these units will be more aggregated in space (see Figure 10). Removing these spatially clustered threats contributes to accruing biodiversity benefits in fewer units. For example, the averaged biodiversity benefits of a unit in local Moran’s I hotspots is 0.08%, which is four times higher than that of units selected using CEA method under the same investment (so does the Getis-Ord Gi* hotspots). Therefore, spatially dispersed actions will make more units benefit from the removal of identified threats, but each unit will have fewer improvements in habitat conditions. Comparatively, if the identified threats are relatively more clustered in space, the same amount of biodiversity benefits will be accrued in fewer units, thus each unit will present higher improvements in habitat conditions. In this context, spatially aggregated threat mitigation actions could more efficiently restore the degraded habitat components (patches and corridors) in ecological networks [40].

4.3. Differences between Sites from Alternative Spatial Clustering Statistics

Differences in the cost effectiveness and the spatial aggregation of prioritized units are mainly determined by which spatial clustering statistic is applied. The similarity in algorithms of Getis-Ord Gi* and local Moran’s I explains the spatial consistency between the statistically significant clusters detected by the two methods. The major difference between these two methods shows why the cost effectiveness of the local Moran’s I hotspots is larger than that of the Getis-Ord Gi* and AMOEBA hotspots. Although the AMOEBA method is adapted from the Getis-Ord Gi* statistic, hotspots detected by the two approaches are different. The AMOEBA method attempts to maximize the Getis-Ord Gi* statistic when detecting clusters. A unit with relatively high CE ratio might not be included into an existing cluster if units in this cluster already have even higher values, because the incorporation of this high-value unit does not increase the Getis-Ord Gi* statistic of this existing cluster. As a result, more small hotspots with higher values could be detected by AMOEBA, consistent with finding from other studies [57,94]. This is why the cost effectiveness of the AMOEBA hotspots is close to that of the Getis-Ord Gi* hotspots, or even higher, but these AMOEBA hotspots are less spatially aggregated.

The difference in the cost effectiveness of units selected according to local Moran’s I and Getis-Ord Gi* statistic is mainly caused by how units detected as statistically insignificant outliers are ranked. The statistically insignificant z-scores in outlier units from local Moran’s I and Getis-Ord Gi* could be either positive or negative. When using Getis-Ord Gi* statistic, an insignificant outlier unit with larger positive z-score indicates relatively higher values in its neighboring units. However, when using local Moran’s I statistic, a larger positive z-score for an insignificant outlier unit could indicate either a unit with larger value surrounded by units with relatively higher values or a unit with smaller value surrounded by units with relatively lower values. If most units in the study area have relatively low cost effectiveness (which is often the case in threat management planning), ranking insignificant outlier units by local Moran’s I statistic will make more units with extremely low cost effectiveness identified as preferred action sites. This is why accrued cost effectiveness of units ranked by local Moran’s I statistic is even smaller than the random selection scenario when more units are selected.

5. Conclusions

Our study demonstrates a spatially explicit prioritization framework that allows for identifying spatially aggregated cost-effective sites for the mitigation of human-induced threats in a transparent and systematic approach. The implementation logic of the promoted prioritization framework is straightforward, simplified and accountable, and could be easily applied use common geographic information system and spreadsheets. The spatially explicit framework is based on the integration of a basic CEA with spatial clustering statistics. To guide stakeholders to invest limited resources to the most cost-effective units for threat mitigation, units are ranked by their CE ratio calculated through the CEA method. Furthermore, spatially aggregated cost-effective actions are informed by directly ranking units in order of autocorrelation-based spatial clustering statistics, including Getis-Ord Gi*, local Moran’s I, and AMOEBA.

The CEA method identified the most cost-effective threat management units, compared with traditional methods that only focus on biodiversity benefits, threats, or costs. However, these high-priority units from CEA are apparently dispersed in space, especially when small initial investment is budgeted. Yet, with the same investment, spatial clustering statistics, such as Getis-Ord Gi* and local Moran’s I could identify spatially clustered management units with only minor loss in estimated cost effectiveness, compared with CEA method. Notably, these spatially clustered units contribute to reducing actual costs related to a given amount of biodiversity benefits and these benefits are placed in fewer units, which is more helpful to restore the connectivity and contiguity of wildlife habitats. The cost-effectiveness of selected sites and the identified hotspots from different spatial clustering statistics differ from each other. Comparatively, local Moran’s I statistic is preferred if the investment budget is very limited. When large investment is placed, Getis-Ord Gi* statistic may be preferred to identify threat mitigation actions that are both spatially clustered and cost-effective. Besides, a larger neighborhood in the spatial clustering statistics will identify actions that are more clustered in space, but will also lead to more decrease in the cost-effectiveness. So, stakeholders have to compromise between different neighborhood configurations regarding the spatial and non-spatial aspects of their expectations.

Our work provides quick and transparent support for informing stakeholders on potential sites for actions to mitigate human-induced threats to biodiversity from the perspective of cost efficiency and spatial aggregation [14]. Other considerations, such as social acceptability and government preference, could be quantified and incorporated into our framework in future work by, for example, adding weight schemes during the calculation of CE ratio. Heterogeneous monetized removal price for each grid of the defined threats could be integrated as well in the future when relevant data become available, and the spatial occurrence and absence of species data in NRs.

Supplementary Materials

The following are available online at https://www.mdpi.com/2071-1050/11/5/1346/s1, Figure S1: Hexagonal units in the protected areas, Figure S2: Land use/cover map of the study area, Appendix A: Evaluation of expected biodiversity benefits, Appendix B: Mathematical details for the three autocorrelation-based spatial statistics.

Author Contributions

J.Y., J.G., and W.T. conceived and designed the study and experiments. J.Y. performed the experiments. J.Y. collected and analyzed the data. J.Y., J.G., and W.T. wrote the paper.

Funding

This study was supported by the national social sciences founding project (14BJY057) and the project of Cropland Quality Evaluation in Qinghai Province.

Acknowledgments

The authors thank Jingye Li, Yuling Liu, and Jia Hao for their assistance. The authors would like to thank the anonymous reviewers for their insightful comments and suggestions.

Conflicts of Interest

The authors declare no conflicts of interest.

Figures and Table

Generate image description

Figure 1. Spatial extent and function zones of protected areas on the Qinghai Plateau.

Generate image description

Figure 2. Framework of prioritizing spatially aggregated cost-effective sites to mitigate human-induced threats (LULC, land use and land cover; HQ, habitat quality; HD, habitat degradation).

Generate image description

Figure 3. Cost effectiveness curves for units ranked in order of CE ratio and traditional prioritization indicators (CE, cost effectiveness).

Generate image description

Figure 4. Spatial distribution of priorities for identified sites using different prioritization indices under given budgets (PA, protected area; CEA, cost effectiveness analysis).

Generate image description

Figure 5. Spatial distribution of hotspots identified by Getis-Ord Gi*, local Moran’s I, AMOEBA, and units selected by CEA method when 13.72% of cumulative costs is budgeted.

Generate image description

Figure 6. Cost effectiveness curves for units ranked in order of CE ratio and different spatial clustering statistics (CE, cost effectiveness).

Generate image description

Figure 7. Spatial distribution of priorities for identified sites using Getis-Ord Gi* and AMOEBA (A Multidirectional Optimum Ecotope-Based Algorithm) statistics under given budgets.

Generate image description

Figure 8. Number and mean area of landscape patches of hotspots under given budgets. (a) the total number of patches in the hotspots detected under different budgets; (b) the mean area of patches in the hotspots detected under different budgets.

Generate image description

Figure 9. Accrued benefits, number and mean area of landscape patches in selected units when the neighborhood distance for Getis-Ord Gi* statistic increases. (a) the accrued biodiversity benefits of units identified using the Getis-Ord Gi statistic with different neighborhood distances; (b) the number of patches in units identified using the Getis-Ord Gi statistic with different neighborhood distances; (c) the mean area of patches in units identified using the Getis-Ord Gi* statistic with different neighborhood distances.

Generate image description

Figure 10. Human land uses with the high priority to be mitigated under a small initial budget of 10% investment

Table 1

Features of hotspots detected by different spatial clustering methods compared with CEA under the same budgets.

Method	Number of Units	Number of Patches	Mean Area of Patches (km2)	Cumulative Costs (%)	Accrued Benefits (%)	Cost Effectiveness Ratio
Local Moran’s I	1117	268	10.83	8.96	22.29	248.76%
CEA	1860	956	5.05	8.96	23.55	262.72%
AMOEBA	2254	1029	5.69	12.17	25.22	207.22%
CEA	2639	1211	5.66	12.17	29.67	243.78%
Getis-Ord Gi*	2148	336	16.61	13.72	28.50	207.70%
CEA	2934	1302	5.85	13.72	32.49	236.79%