Landscape Conservation Planning to Sustain Ecosystem Services under Climate Change

Sustainable conservation aims to ensure the sustained conservation of landscape multi-functionality which in turn requires ensuring ecosystem service (ES) and habitat quality (HQ) sustainability with inclusive landscape-scale conservation planning. This study proposes a landscape conservation planning (LCP) framework for landscape-scale ES-HQ conservation and sustainability. Spatially explicit hotspots for five ESs and HQs are identified via InVEST and LISA software. Spatiotemporal changes in ES-HQ hotspots, in terms of stability and resilience, are delineated. The Zonation technique is applied to prioritize areas for conservation based on ES-HQ hotspot stability and resilience maps. High priority conservation areas are identified and are used as reserve area inputs for land use modeling with CLUE-S software to simulate future land use change under climate change scenarios. This study reports that varied rainfall and climate are major driving factors of ES-HQ sustainability disturbance in the study area. Furthermore, our proposed conservation Strategy 2 demonstrates that a larger extent of landscape multi-functionality can be sustained when the existing conservation area includes the total area of identified ES-HQ resilient hotspots. This study effectively identifies the stability and resiliency of ES-HQ hotspot areas affected by disturbances for high priority landscape conservation requirements to ensure ES-HQ sustainability and landscape multi-functionality in the study area.


Introduction
Sustainable conservation, as a concept and strategic aim, includes a number of significant challenges since ensuring the sustained conservation of landscape multi-functionality necessitates ensuring ecosystem services (ES) and habitat quality (HQ) are also sustained, with inclusive landscape-scale conservation planning.While it is well-understood that ES are the set of ecosystem functions that benefit people [1][2][3][4] by providing natural resources, regulating ecological processes, and providing social value [5], the urgency to conserve ES and also to sustain the landscape quality that supports their functionality for our well-being [2,[6][7][8][9] remains.Ecosystem functions generate ES and are supported by biophysical structures and processes [6], while landscape functions are the aggregated ecosystem functions and services at the landscape scale.Though 'ES conservation' is mainly concerned with ES, the need for a broader focus [10] to include landscape sustainability that supports ES functionality has been reported.ES conservation concepts have been applied to sustainable landscape development studies and planning [11], and to understanding how to conserve and sustain those benefits derived from nature [12].Recently, numerous studies have used an ES conservation-based approach to conservation planning.For example, Bhagabati et al. [13] determined conservation scenarios based on ES; Sandifer et al. [14] used ES for biodiversity health enhancement and conservation; Rodriguez-Loinaz et al. [15] proposed multifunctional landscape conservation using an ES conservation approach; Hansen et al. [16] included ES in planning discourses of European and American cities; Mitchell et al. [17] considered ES in their landscape management approach under landscape fragmentation; Darvill and Lindo [18] reported on conservation with priorities for cultural ES values; Zheng et al. [19] considered ES with respect to water conservation policies and management practices; Snall et al. [20] considered ES in green infrastructure design; Mupepele et al. [21] included ES in evidence-based environmental management; Jorgensen et al. [22] reported on local municipal decision making in South Africa that included ES considerations; Maes et al. [4] applied ES conservation concepts in support of the EU Biodiversity Strategy; Verhagen et al. [23] reported on spatial priority patterns for ES conservation; and Lin et al. [10] incorporated ES in social-ecological conservation.Yet, 'landscape sustainability' itself is an important research priority that is a key in shaping the debate among landscape ecologists about the relationships between landscapes, ES, and human well-being [24].
In order to integrate ES conservation and landscape sustainability, research definitions of landscape sustainability are centered on natural capital and ES [25][26][27].Furthermore, landscape sustainability as an end requires maintaining the dynamic relationship between ES and human well-being in landscapes that have changed due to internal and external disturbances and their uncertainties [28].While abiotically-or biotically-caused landscape disturbances create complex heterogeneous patterns of varying severity levels among the affected areas [27], 'environmental resilience' is the ability of an ecological system to absorb and recover from these disturbances in an efficient manner.Landscape sustainability as a measure then, is the capacity of a given landscape to provide constant and long-term, landscape-specific ESs that are essential to our well-being, within a regional context and despite environmental and sociocultural changes [28].Accordingly, sustainable landscapes are those landscapes that can consistently provide long-term, landscape-specific ESs necessary for maintaining and improving human well-being [28].
It is essential to understand the impacts of natural disturbances on landscape functions [29] such as ES, in order to ensure ES within sustainable landscapes.Generally, environmental 'resistance' is characterized as the influence of an ecological community's structure and composition on disturbance, whereas 'resilience' is characterized as the influence of disturbance on an ecological community's subsequent structure and composition [30].Theoretically, lower environmental resilience results in an amplified response to disturbances rendering the area more sensitive to environmental perturbations [31].Additionally, slower responses to disturbances may be evidence of reduced recovery rates in areas approaching critical states [31].Therefore, identifying landscapes with reduced recovery rates, or high ecological sensitivity, is an important step in recognizing regions pending critical ecological change [31].Since ES conservation and landscape sustainability converge at environmental resistance and resilience thresholds (Figure 1).These thresholds then can be considered to be the state of a landscape's functions under disturbances.
Landscape conservation planning (LCP) is a conservation planning approach for efficiently designing conservation areas to comprehensively protect a representative species population and its habitat [32].Due to the intrinsic value of ES [2], LCP-prioritized conservation areas should be identified to ensure the available supply of ES and to meet demands for ES.LCP [33,34], however, has been widely applied to improve existing reserve networks while considering biodiversity patterns [10,[35][36][37][38][39].Though LCP is a well-rounded approach for preserving the natural value of terrestrial [10,[36][37][38]40], marine [41,42], and freshwater [43,44] ecosystems, its main purpose is for locating, designing, and managing protected areas within a landscape [45].Recently, interest in spatially quantitative systematic conservation approaches has been growing, as evidenced by the increased applications of computer-based software modeling tools such as Marxan [46] and Zonation [47] which both implement target-based planning algorithms as their primary planning method [48].Sustainable conservation.Sustainable conservation can be achieved using a systematic conservation planning approach to target ecosystem service resilience, the key to both ecosystem service conservation and landscape sustainability.
Landscape conservation planning (LCP) is a conservation planning approach for efficiently designing conservation areas to comprehensively protect a representative species population and its habitat [32].Due to the intrinsic value of ES [2], LCP-prioritized conservation areas should be identified to ensure the available supply of ES and to meet demands for ES.LCP [33,34], however, has been widely applied to improve existing reserve networks while considering biodiversity patterns [10,[35][36][37][38][39].Though LCP is a well-rounded approach for preserving the natural value of terrestrial [10,[36][37][38]40], marine [41,42], and freshwater [43,44] ecosystems, its main purpose is for locating, designing, and managing protected areas within a landscape [45].Recently, interest in spatially quantitative systematic conservation approaches has been growing, as evidenced by the increased applications of computer-based software modeling tools such as Marxan [46] and Zonation [47] which both implement target-based planning algorithms as their primary planning method [48].
While sustainable landscape planning places physical planning in a broader perspective [49,50], this study uses LCP to prioritize conservation areas for sustainable landscape management in terms of ES.With a focus on sustaining ES after disturbances for landscape sustainability planning, we assess ES change within a landscape that is a spatially heterogeneous, human-environment coupled system.Specifically, the objective of this study is to systematically assess ES as a landscape function within a post-disturbance landscape using stability and resilience analyses to identify changes in ES.

Materials and Methods
For sustainable landscape planning, this study uses the proposed "landscape conservation planning" (LCP) framework to ensure sustained ES, and model land use change to evaluate the effects of future ES changes.This framework was developed in terms of ES stability and resilience and is depicted in Figure 2 as a flow chart.Analyses in this study were carried out as evaluations of ES and HQ at three levels.The Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) software (https://www.naturalcapitalproject.org/invest/) [51] (see Supplementary S1.1) was used to investigate resilience causation and determine resilience scores.As an evaluation tool, InVEST can be used to calculate the extent to which land use and weather changes affect ES and HQ.We further defined 'ES stability' and 'ES resilience' for this study as the ecosystem's ability to return to a predisturbance state within the Chenyulan Watershed, located in the central part of Taiwan (Figure 3).
The Chenyulan River's watershed (449 km 2 ), a typical mountain drainage watershed, is located in Nantou County with a main stream length of 42 km, average stream-bed gradient of 4°, and elevation range of 310-3952 m.The Chenyuland watershed has mean slope of 32°, and relief of 585 m/km.The slates and meta-sandstones are the dominant lithologies in the metamorphic terrains Sustainable conservation.Sustainable conservation can be achieved using a systematic conservation planning approach to target ecosystem service resilience, the key to both ecosystem service conservation and landscape sustainability.
While sustainable landscape planning places physical planning in a broader perspective [49,50], this study uses LCP to prioritize conservation areas for sustainable landscape management in terms of ES.With a focus on sustaining ES after disturbances for landscape sustainability planning, we assess ES change within a landscape that is a spatially heterogeneous, human-environment coupled system.Specifically, the objective of this study is to systematically assess ES as a landscape function within a post-disturbance landscape using stability and resilience analyses to identify changes in ES.

Materials and Methods
For sustainable landscape planning, this study uses the proposed "landscape conservation planning" (LCP) framework to ensure sustained ES, and model land use change to evaluate the effects of future ES changes.This framework was developed in terms of ES stability and resilience and is depicted in Figure 2 as a flow chart.Analyses in this study were carried out as evaluations of ES and HQ at three levels.The Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) software (https://www.naturalcapitalproject.org/invest/) [51] (see Supplementary S1.1) was used to investigate resilience causation and determine resilience scores.As an evaluation tool, InVEST can be used to calculate the extent to which land use and weather changes affect ES and HQ.We further defined 'ES stability' and 'ES resilience' for this study as the ecosystem's ability to return to a pre-disturbance state within the Chenyulan Watershed, located in the central part of Taiwan (Figure 3).
The Chenyulan River's watershed (449 km 2 ), a typical mountain drainage watershed, is located in Nantou County with a main stream length of 42 km, average stream-bed gradient of 4 • , and elevation range of 310-3952 m.The Chenyuland watershed has mean slope of 32 • , and relief of 585 m/km.The slates and meta-sandstones are the dominant lithologies in the metamorphic terrains [52,53].The annual regional rainfall in the watershed is between 2000 and 5000 mm, with an average of approximately 3500 mm.Typhoons which occur three to four times annually contribute to the annual rainfall levels, 80% of which generally occur between May and October during the rainy season.Because of its steep topography, loose soils and short yet intensive periods of rainfall, it is a consistent historical location for extreme rainfall events and related events, such as landslides, denudation in the mountainous area or the debris flow (Jan & Chen, 2005) [54].On September 21, 1999, a 7.3 magnitude on the Richter scale earthquake with a focal depth of 8.0 km occurred.The "Chi-Chi earthquake" was Taiwan's largest earthquake occurrence in 100 years, resulting in numerous landslides in the Chenyulan Watershed, and reoccurrences still [55].
In order to map the distribution of ES and HQ, the Local Indicator of Spatial Association index (LISA) [10,56] (see Supplementary S1.2) was used to identify ES in five time periods from 1999 to 2005, and HQ hotspots.The study period is divided into 'time periods' as opposed to 'years' since the data: 1) includes data from two separate months (March 1999 and October 1999) which are before and after the Chi-Chi earthquake and are treated as stand-alone time periods; and 2) does not include data for the years of 2002 or 2003.LISA was then used to calculate ES-HQ hotspots for each of the five time periods, and to create stability and resilience maps for each ES.Once maps were drawn, the Zonation modeling tool [57] (see Supplementary S1.3) was used to identify which ES-HQ hotspot areas meet conservation targets.Once these areas were identified, five Zonation landscape protection strategies were applied to recalculate resilience maps.This step requires calculating the number of times a map grid cell's ES hotspot status changed between 1999 and 2005, when each cell's number (from 0-5) represents ES stability and resilience.A grid cell value of 0 denotes no change in ES status as a constant hotspot or non-hotspot, while a cell value greater than 1 indicates the number of times ES status has changed.
Cumming [58] recommends that a baseline against which to measure change, and criterion to define success or failure (in terms of sustaining or conserving ES) is needed.Accordingly, this study uses the March 1999 (henceforth "1999a") ES-HQ hotspot data as a benchmark to assess the resultant scenarios.Based on the recalculated resilience maps, we selected five landscape protection strategies for ES conservation, or, ES conservation strategies: (1) Strategy 1 is the control conservation strategy in which there is no change to the original reserve areas (i.e., Yushan National Park); (2) Strategy 2 is a conservation strategy in which the original reserve area is expanded to include the total area of Zonation-identified ES-HQ resilient hotspot areas from the benchmark data; (3) Strategy 3 is a conservation strategy in which the original reserve area is expanded to include the top 10% of the total Zonation-identified ES-HQ hotspot areas from the benchmark data; (4) Strategy 4 is a conservation strategy in which the original reserve area is expanded to include the top 20% of the total Zonation-identified ES-HQ hotspot areas from the benchmark data; and (5) Strategy 5 is a conservation strategy in which the original reserve area is expanded to include the top 30% of the total Zonation-identified ES-HQ hotspot areas from the benchmark data.In addition, in the reserve area, the forest area is forbidden to change to build-up area or agricultural area.Following this, the conservation strategies are then used as land use change model inputs to simulate land use in 2040 under different future climate environments.The Conversion of Land Use and its Effects at Small regional extent model (CLUE-S) [59] (see Supplementary S1.4) was used to simulate 2040 land use as four scenarios (i.e., Scenarios 1-4) based on the five ES conservation strategy inputs (i.e., Strategies 1-5) to consider the impact of climate change on the sustainability of ES.Specifically, the predicted land allocation algorithm in this study's CLUE-S model takes into account location-specific factors, spatial policies, land use conversion trajectories, and the determined calculation method for land use demand for each scenario.The Kappa statistic is also used to evaluate the level of agreement between datasets as a result of the last step.
For climate change scenarios, two representative concentration pathways (RCP), RCP 2.6 (i.e., low emissions) and RCP 8.5 (i.e., high emissions) of five general circulation models (GCMs), i.e., CCSM4, CESM1-CAM5, GISS-E2-R, HadGEM2-AO, and MIROC5, from the Intergovernmental Panel on Climate Change (IPCC) were adopted in this study.The original IPCC GCM output data, including changing ratio of monthly rainfall and increasing monthly temperatures, was downscaled to 25 × 25 km resolution with a statistical downscaling method used by the Taiwan Climate Change Projection and Information Platform (TCCIP).The climate data for the Chenyulan Watershed study area were then further downscaled to a 100 × 100 m resolution using an area-to-point (ATP) co-kriging technique (see Supplementary S1.5).The future rainfall and temperature for each cell was then calculated from historical data and downscaled climate data."No-change" data indicates the historical average trend   For climate change scenarios, two representative concentration pathways (RCP), RCP 2.6 (i.e., low emissions) and RCP 8.5 (i.e., high emissions) of five general circulation models (GCMs), i.e., CCSM4, CESM1-CAM5, GISS-E2-R, HadGEM2-AO, and MIROC5, from the Intergovernmental Panel on Climate Change (IPCC) were adopted in this study.The original IPCC GCM output data, including changing ratio of monthly rainfall and increasing monthly temperatures, was downscaled to 25 × 25 km resolution with a statistical downscaling method used by the Taiwan Climate Change Projection and Information Platform (TCCIP).The climate data for the Chenyulan Watershed study area were then further downscaled to a 100 × 100 m resolution using an area-to-point (ATP) co-kriging technique (see Supplementary S1

Changes in Ecosystem Services Induced by Large Disturbances
Hotspot analysis was further conducted on the ES and HQ outputs generated by the InVEST model (i.e., water yield, nitrogen retention, phosphorus retention, sediment retention, and carbon storage ES) for each period when large disturbances occurred.Specifically, for all ES and HQ (S2a in Supplementary) from 1999 to 2005, we used LISA statistics to identify grid cells as hotspots.Generally, it was found that ES and HQ hotspots (S2b in Supplementary) were highly correlated with rainfall, landscape, or elevation, which tend to be altered by large disturbances.Water yield hotspot locations are highly dependent on rainfall distribution, especially in the southwest of the watershed (Figure S10 in Supplementary).Water yield hotspots range from 27.87% to 34.76% of the watershed during the 1999-2005 period (Table 1).Nitrogen retention hotspots tend to occur in the vicinities of agriculture lands which are one of the main sources of nitrogen emission (Figure S11 in Supplementary).Nitrogen retention hotspots range from 2.39% to 3.66% of the watershed during the 1999-2005 period (Table 1).Unlike nitrogen retention hotspots, agricultural land distribution does not impact the locations of phosphorus retention hotspots.Instead, phosphorus retention hotspot locations tend to be correlated with elevation, and especially in the low-elevation areas (Figure S12 in Supplementary).Phosphorus retention hotspots range from 1.61% to 2.67% of the watershed during the 1999-2005 period (Table 1).
The distribution of sediment retention hotspots during the 1999-2005 period is consistent, ranging from 6.27% to 7.60% of the watershed (Table 1 and Figure S13).Carbon storage hotspots are mostly located in forest lands since forest lands have higher carbon density in aboveground mass, belowground mass, soil, and dead mass than other land use types (Figure S14 in Supplementary).Carbon storage hotspots range from 39.48% to 43.69% of the watershed during the 1999-2005 period (Table 1).Similarly, forest lands provide habitats; thus, habitat quality hotspots are mostly located in forest lands ranging from 51.68% to 52.51% of the watershed (Table 1 and Figure S15).Furthermore, less than 36% (33.31-35.93%) of the watershed provides no ES, while the remaining areas (60.37-63.09%)provide 1-3 ES and/or HQ during the 1999-2005 period (Table 2).Note: Ecosystem service (ES); 1999a refers to the March 1999 stand-alone period; 1999b refers to the October 1999 stand-alone period.

Ecosystem Service Stability and Resilience
This study defined ES stability and resilience as the ability of the ecosystem to return to a pre-disturbance state.First, we calculated the number of times a grid cell's ES hotspot status changed between 1999-2005 and designated a number of 0-5.A grid value of 0 denotes no change in the grid cell's ES status which has remained a constant hotspot or non-hotspot.A grid value greater than 1, however, indicates the number of times the ES hotspot status has changed over time.When identifying which ES hotspot areas had changed ES hotspot status in the following years since 1999, we found that 7.78% of the water yield ES hotspots, 2.76% of the nitrogen retention ES hotspots, 2.24% of the phosphorus retention ES hotspots, 2.17% of the sediment retention ES hotspots, 11.42% of the carbon storage hotspots, and 3.43% of the habitat quality hotspots within the watershed have encountered at least one ES hotspot status change since 1999 (Figure 4 and Table 3).Furthermore, we found that 71.19% of the watershed was not a water yield hotspot in 1999, while 21.02% of the watershed was a constant water yield hotspot from 1999-2005, (Table 3), indicating that 92.22% of the watershed was stable during the study period in terms of water yield provisioning.It was also found that of the water yield hotspots in the southern part of the watershed that have changed hotspot statuses at least once, 5.68% have changed statuses twice (Figure 4a and Table 3).Of the hotspots with statuses that have changed twice, half of those status changes were the result of high precipitation in 2005.More than 93% of the watershed was a nitrogen retention, phosphorus retention, and sediment retention non-hotspot in 1999 (Table 3).Thus, the number of ES hotspot status changes since 1999 are few.Specifically, for the 1999-2005 period, only 0.37-1.23% of nitrogen retention hotspots and 0.28-1.22% of phosphorus retention hotspot areas experienced 1-3 status changes.
Similar to the correlation between agricultural lands and nitrogen retention hotspots observed in each period, stable hotspot (0.33%) and resilient hotspot (2.76%) areas where the nitrogen retention ES hotspot status changed during the study period are mostly located within those areas where More than 93% of the watershed was a nitrogen retention, phosphorus retention, and sediment retention non-hotspot in 1999 (Table 3).Thus, the number of ES hotspot status changes since 1999 are few.Specifically, for the 1999-2005 period, only 0.37-1.23% of nitrogen retention hotspots and 0.28-1.22% of phosphorus retention hotspot areas experienced 1-3 status changes.
Similar to the correlation between agricultural lands and nitrogen retention hotspots observed in each period, stable hotspot (0.33%) and resilient hotspot (2.76%) areas where the nitrogen retention ES hotspot status changed during the study period are mostly located within those areas where agricultural land use remained constant (1.44%), or mostly constant for five periods (2.61%) (Figure 4b, Tables 3 and 4, and Figure S11).Areas where phosphorus retention hotspot status changed during the 1999-2005 study period are located in higher elevation forest lands and downstream of agricultural lands within the watershed (Figure 4c).The areas where sediment retention ES hotspot statuses changed, were located in steeper sloped areas, where landslides often occur (Figure 4d).Landslide areas comprise 0.65-2.00% of the watershed (Table 5), and thus 4.10% of the watershed that is downstream forest remained a sediment retention hotspot during the study period (Table 3).Moreover, 4.10% of the watershed has experienced mostly constant landslides for five periods during the 1999-2005 study period, resulting in a resilient sediment retention ES hotspot (2.18% in Table 3).While 'landslide' is not a land use type per se, it is included for additional insight since areas in Taiwan are designated as high-risk areas for landslide occurrence as a result of recent tragic landslide disasters.Note: 1999a denotes the March 1999 stand-alone period; 1999b denotes the October 1999 stand-alone period; * Landslide is not a land use type per se, but is included for additional insight since areas in Taiwan are designated as high-risk areas for landslide occurrence as a result of recent tragic landslide disasters.
The areas that were constant carbon storage hotspots account for 28.41% of the study area and were mostly located in the forest lands (Table 3 and Figure S14 in Supplementary), while the carbon storage hotspots that have changed hotspot status from 1999-2005 account for 12.02% of the study area (Table 3 and Figure 4e).Resilient carbon storage hotspot areas were located in forest areas that had changed between (grassland or agricultural) land use type.
The areas that were constant HQ hotspots and resilient HQ hotspots account for 48.96% and 3.54% of the watershed, respectively (Table 3).Similar to the resilient carbon storage hotspots, those resilient HQ hotspots that underwent hotspot status change at least once during the study period were located in forest areas that had changed between (grassland or agricultural) land use type (Figure 4e).When we combined the resilient hotspots of all ES with HQ, we found that a total of 31.38% of the watershed was a hotspot area during the benchmark period, and underwent at least one hotspot status change during subsequent periods (Figure 4g).

Impacts of Ecosystem Service Conservation Strategies and Climate Change Scenarios on the Predicted Landscape
Based on the resultant ES-HQ hotspot maps (S2b in Supplementary) and ES stability and resilience maps (Figure 4), five ES conservation strategies to manage the predicted 2040 landscape were proposed which were then used as inputs for four land use change future scenarios.The reserve areas for each of the five ES conservation strategies are as follows: Strategy 1 with 23.38%; Strategy 2 with 40.75%;Strategy 3 with 30.27%;Strategy 4 with 33.74%; and Strategy 5 with 37.62% of the study area designated as reserve areas (Figure 5 and Table 5).Relative to the March 1999 benchmark landscape, forest lands decrease while agricultural lands and urban areas increase in 2040.The predicted land use scenarios under the above strategies are shown in Figures S16-S19 in Supplementary.

Impacts of Ecosystem Service Conservation Strategies and Climate Change Scenarios on Ecosystem Services
The ES-HQ hotspots in 2040 under five different ES conservation strategies and four future land use change scenarios are shown in S4 and Table S2 (in Supplementary).Although the reserve areas vary for different ES conservation strategies, the total area (as percent of watershed) of hotspots of most ES and HQ did not greatly vary among strategies.For the different strategies and scenarios, total area of hotspots ranged between 25.95-29.84%for water yield, 1.88-3.161%for nitrogen retention, 1.67-2.60%for phosphorus retention, 5.62-6.47%for sediment retention, 35.98-45.82%for carbon storage, and 37.87-51.52%for HQ.However, different strategies and scenarios may affect the spatial distribution of hotspot locations as evidenced in the amount of hotspot overlap between the benchmark period data and all five strategies under four land use scenarios: 21.55-26.34% of the watershed overlap for water yield, 0.61-1.89%for nitrogen retention, 0.46-1.68%for phosphorus retention, 4.82-5.74%for sediment retention, 24.64-37.96%for carbon storage, and 32.26-50.48%for HQ (Table S2 in Supplementary).That is, for each of the five conservation strategies simulated under each of the four climate change-affected land use scenarios, the spatial distribution of hotspot locations may vary and can be seen in the distribution similarities of the mapped outputs.When inspecting the overlap all of the mapped outputs, the ratio of overlap indicates the variations in spatial distribution of hotspot locations as a result of the different conservation strategies and scenarios.For example, 100% overlap would indicate no variation in spatial distribution among strategies and scenarios.
Compared to the total area of water yield hotspots in Strategy 1 which only considered Yushan National Park as a reserve area, the total area of water yield hotspots did not increase in Strategy 2 which included the resilient landscape in its reserve area.As the percentage of ES-HQ hotspots increased (Strategies 3-5), so did the total area of water yield hotspots increase.A similar positive relationship between hotspot areas and reserve areas was found between carbon storage and HQ (Table S2 in Supplementary).Carbon storage hotspot areas increased slightly from 36.17-44.62% to 37.73-44.95% in Strategy 3 versus Strategy 5, as did HQ hotspot areas (37.96-51.26% to 39.56-51.32% in Table S2 Supplementary).For water yield, carbon storage, and HQ, over 94.34-99.21% of hotspots in Strategy 3 were overlapped with those in Strategy 1.
Only half of the nitrogen retention or phosphorus retention hotspots in Strategy 1 overlapped with those in Strategy 2, as well as Strategy 3. When comparing Strategies 1 and 2, the difference in the reserve area totals was equivalent to the total area of resilient hotspots in Strategy 2. When comparing Strategies 1 and 3, the difference in the reserve area totals was equivalent to the top 10% of the Zonation-identified ES-HQ hotspots.
(Figure 4(e)).When we combined the resilient hotspots of all ES with HQ, we found that a total of 31.38% of the watershed was a hotspot area during the benchmark period, and underwent at least one hotspot status change during subsequent periods (Figure 4(g)).

Impacts of Ecosystem Service Conservation Strategies and Climate Change Scenarios on the Predicted Landscape
Based on the resultant ES-HQ hotspot maps (S2b in Supplementary) and ES stability and resilience maps (Figure 4), five ES conservation strategies to manage the predicted 2040 landscape were proposed which were then used as inputs for four land use change future scenarios.The reserve areas for each of the five ES conservation strategies are as follows: Strategy 1 with 23.38%; Strategy 2 with 40.75%;Strategy 3 with 30.27%;Strategy 4 with 33.74%; and Strategy 5 with 37.62% of the study area designated as reserve areas (Figure 5 and Table 5).Relative to the March 1999 benchmark landscape, forest lands decrease while agricultural lands and urban areas increase in 2040.The predicted land use scenarios under the above strategies are shown in Figures S16-S19 in Supplementary.

Driving Factors of Changes in Ecosystem Services and Habitat Quality
Numerous factors determine a landscape's succession and recovery rate including species life, community interaction, soil, climate, and the complexities of disturbance legacies [60].For sustainable landscape planning, gaining a better understanding of resilience can quantitatively change disturbance regimes [60].Similar to Thom and Seidl's [29] study, our results indicate that ES such as supporting, provisioning, and regulating services were predominately negatively affected by disturbances.By systematically quantifying factors that focus on ES-HQ resilience under disturbances within the study area, our results indicate that climatic changes may have negative impacts on ES provisioning.Specifically, climactic changes may have severe impacts on 'service-providing units, intermediate-, and final-' ES [12,29].However, short term climactic changes such as spatiotemporal changes in precipitation within the study area, along with land use types can impact water yield ES hotspots [61].
Moreover, our study indicates that varied rainfall and climates induced by large typhoon systems can be considered major driving factors of disturbances.This is in part due to the fact that rainfall on sloped landscape results in flow pathways which are important factors in nitrogen retention and phosphorus retention hotspot resilience [38,61,62].In addition, phosphorus retention hotspots could be affected by land use types such as forested lands and downstream sediment-attached phosphorus deposition (Figure S12 in Supplementary).Topographic features, therefore, are the predominant factors which determine the location of sediment retention hotspots (Figure S13 in Supplementary).In this study, land use type is the most important factor for carbon storage and HQ hotspot resiliency.Furthermore, since forest land provides the most varied habitat among all land use types (Figures S14 and S15 in Supplementary), the impacts of disturbances on carbon storage and HQ will differ [29] though forest is the most significant land use type for either hotspot or carbon storage resilience capacity.

Differences in Ecosystem Service and Habitat Quality Resilience
Ensuring the sustainability of ES provision is a challenging undertaking due to the short-to mid-term changes, and the long-term landscape reorganization, that disturbance regimes affect [60].Similar to the 2005 period phosphorus retention hotspots identified by LISA (Figure S11 in Supplementary), phosphorus retention resilient hotspots were located where cultivated and forested land were in land use type transition, indicating that phosphorus retention hotspot locations did not greatly vary though hotspot status frequently changed during the 1999-2005 study period (Figure 4c).Moreover, the sediment retention hotspot areas that were more resilient to disturbance are more pronounced as depicted in Figure 4d and Figure S13 in Supplementary.Additionally, those areas where carbon storage and habitat quality hotspot status changed more than twice were located on forest lands where the boundaries of agricultural lands and forest changed more frequently (i.e., land use type changed) from 1999-2005.While previous studies have indicated that we can conclude from ES resilience analyses results that landscape stability and resilience could significantly influence ES and HQ [8], this study further demonstrates that landscape stability and resilience does influence ES and habitat quality.
Nitrogen retention resilient hotspots were usually located in agricultural land because agricultural land generates a relatively larger nitrogen load (16,000 kg/ha/year) with minimal retention efficiency (5%) when compared to forest land nitrogen loads (1,600 kg/ha/year) with 80% retention efficiency.Although the watershed during the 1999-2005 study period was comprised of only 10.2-13.66%agricultural land (Table 5), agricultural land retained the most nitrogen when compared to forested land that comprised 73.25-78.6% of the watershed.Phosphorus retention hotspots were usually located in the higher elevation forest lands and in downstream agricultural lands within the watershed.In this study, we assumed that agricultural land generates 500 kg/ha/year of phosphorus with a 5% retention efficiency, while forest land generates 250 kg/ha/year of phosphorus with an 80% retention efficiency.Thus, forest land plays an important role in the watershed's upstream phosphorus retention capacity while agricultural land retains more accumulated phosphorus downstream since, unlike the relationship between agricultural land and nitrogen retention, forest land retained more phosphorus than agricultural land.Our results demonstrate how agricultural land directly impacts nitrogen and phosphorus ES hotspots [7] under disturbances.
Upstream sources such as erosion-prone areas transport sediment downstream by hydrologic flows [63].Sediment retention hotspots then, are greatly affected by precipitation induced soil movement, soil type erodibility, slope, erosion, protection of land use types, and land management practices [13].In this study, sediment retention hotspots were mostly located at the boundaries between tributary streams and forest land.Specifically, at steeper sloped areas where landslides often occur.For the 1999-2005 study period, it was found that in most cases, landslides occurred on grassland and forest land, with an average of 0.14% of the watershed experiencing constant landslides for more than five periods.Because of this, landslides were a relatively stable 'land use type' (see Table 5 note) which rendered the sediment retention hotspots, and their statuses, resilient.Within the study area, constant landslides account for 42% of the grassland and forest land; continuous landslides account for 32% when landslides occurred in grassland and forest land from 1999-2001 and 2005, but not in 2004; continuous landslides account for 15% when landslides occurred in grassland or forest land from October 1999 onward.
Although forest land comprised 73.25-78.60% of the watershed (Table 5), carbon storage hotspots were not necessarily as large as the total area of forest land.During the 1999-2005 study period, 57.16% of the watershed has continuously remained forest land, while 8.68% of the watershed remained forest land for exactly five of the six periods.Transitions between forest-to-grassland and forest-to-agricultural land are mostly observable during the image processing stage.This discrepancy in land classification mainly resulted in the designation of 9.34% of the watershed as a resilient carbon storage hotspot that has undergone a hotspot status change 1-2 times (Table 3).

Impact of Ecosystem Service Conservation Strategies and Climate Change Scenarios on Ecosystem Service-Habitat Quality Hotspots
Maintaining landscape multi-functionality (i.e., multiple ESs within high-quality habitats) can be achieved when landscape functions are 'organized and spatially stacked' [49] so that ES functions at various scales and their physical, chemical, and biological processes [64] are maintained and unimpeded.Wu [28] proposed an integrated landscape sustainability science (LSS) framework for sustainable landscape management that included natural capital, ES, ecosystem valuation, disturbances, ecosystem resilience, ecological restoration, and ecosystem management.In this study, a systematic conservation approach can enhance the proposed ES-HQ conservation framework [38].Our proposed framework, however, provides a systematic method for assessing conservation strategies and statuses specific to, and within a, Taiwan watershed affected by disturbances by identifying landscape responses to disturbances in terms of ES-HQ resilience.Our study indicates that for the study area, increases in agricultural and urban land by 2040 will result in decreased water yield, nitrogen retention, phosphorus retention, sediment retention, and habitat quality hotspot areas (Table S2 in Supplementary).
Prioritization may be needed, along with regulatory changes, to ensure that the intended land use type distribution is made spatially explicit [8].Although forest land is expected to decrease by 2040, with the proposed protection strategies, carbon storage hotspot areas can be retained and even increased to an area greater than the benchmark carbon storage hotspot total area.Water yield hotspot areas increase slightly in Strategy 3 versus Strategy 5 (26.88-28.15%vs. 27.57-29.84%,respectively) as the reserve areas increase (30.27% in Strategy 3 vs.37.62% in Strategy 5), relative to the reserve areas of Strategy 1 (26.43-28.15%).A similar positive relationship between hotspot areas and reserve areas was found in carbon storage and habitat quality (Table S2 in Supplementary).This is due to the fact that forest distribution is the key driving factor for carbon storage and habitat quality hotspots.As the reserve areas increase, more existing forest is preserved within the protected areas.Moreover, when Strategy 1 and Strategy 2 maps are spatially combined to create a single map that delineates resilient hotspots within the conservation area; then when this single map is compared with a Strategy 2 map alone (with delineated resilient hotspots), the ratio that the Strategy 2 map overlaps the combined map (Table S2 in Supplementary) indicates that Strategy 2 protects a larger extent of resilient hotspot areas than Strategy 1.When using the same comparison method for Strategy 1 and Strategy 3, Strategy 3 protects a larger extent than that of Strategy 1. Strategies 2-5 then, can more efficiently protect a larger extent of resilient hotspot areas than Strategy 1 since Strategies 2-5 extend the control conservation area of Scenario 1.
Although agricultural development may impact watershed water quality [19,61], it may also impact sediment retention and water purification so that both capacities increase under agricultural development when compared to the benchmark measurements.The differences in reserve areas resulted in different distributions of nitrogen retention and phosphorus retention hotspots for Strategies 2 and 3 (Figures S14 and S15 in Supplementary).Because Strategies 2 and 3 nitrogen retention hotspots range from 42.92-68.93%,the nitrogen retention hotspot of Strategy 2 is quite similar to that of Strategy 3, though dissimilar to that of Strategy 1.Although the phosphorus retention hotspot distributions of Strategies 1 and 3 are dissimilar, the hotspot total areas are similar.This indicates that by using Strategy 3, additional phosphorus retention hotspots can be sustained.

Conclusions
Landscape management decisions can alter ecosystem functioning and ES provisioning capacities [7].Our analysis results demonstrate that accounting for sustainable ES-HQ hotspots in systematic conservation planning at the watershed scale is possible.Additionally, the proposed framework and methodology is a novel application of systematic conservation planning for sustaining ES and HQ in the study area.Unlike previous studies [4,10,[13][14][15][16][17][18][19][20][21][22][23], our SLCP approach effectively identified both resilient ES-HQ hotspot areas and those affected by disturbances with an interdisciplinary application of scenario planning and resilience concepts [8].Using modeling tools, we then prioritized conservation areas to create various planning scenarios.Furthermore, by systematically modeling land use with conservation strategies and benchmark ES measurements, ES-HQ hotspots can be accurately represented for effective predictive analysis.While our framework may not be suited for assessing conservation strategies at a broader scale, future studies using Zonation conservation prioritization techniques along with CLUE-S modeling may benefit from our approach.The proposed approach does provide, however, a systematic way to manage ES and HQ sustainability with conservation strategies for landscape sustainability planning.Stakeholder engagement can further enhance planners' understanding of the locally lived realities, such as pre-and post-disturbance effects on agricultural activities to determine which conservation strategies are best suited for an area [60,65].Although our results are acceptable for ES-HQ hotspot sustainability management, and though our model used parameters from previous studies [10,61] specific to Taiwan, land modeling uncertainties and stakeholder engagement should be taken into account when determining reserve areas.Validation of all ES calculations should be conducted to enhance ES modeling results.
Maintaining landscape multi-functionality-that is, its ES-can be achieved when landscape functions are identified and protected for sustainable landscape planning.That is, when landscape-essential ES and HQ needed for human well-being are resilient to change over time with appropriate sustainable landscape planning, these ES and HQ functionality can be sustained and thereby a landscape's multi-functionality.In this study, we identified areas in which ES-HQ stability and resilience were maintained at certain levels after multiple disturbances.In our ES-HQ hotspot analysis, the consequent disturbances after earthquakes and large typhoons with varied rainfall and climate, are considered to be major driving factors that disturb the resilience and stability of the spatially explicit ES-HQ in the study area.Moreover, land use type transitions, such as forest land changes after disturbances, are driving factors for carbon storage and HQ hotspot resiliency.Without selecting areas within a landscape for protection based on the stability and resilience of landscape function hotspots, landscape functions in terms of ES-HQ resilience may not effectively be ensured.Unlike previously proposed landscape planning approaches, this study spatiotemporally identifies resilient and stable ES-HQ hotspots for designing ES conservation scenarios under climate change.Furthermore, our SLCP framework is based on spatiotemporal stability and resilience analyses to identify ES-HQ hotspot areas affected by disturbances as well as climate change for sustaining landscape multi-functionality.

Sustainability 2018 , 19 Figure 1 .
Figure1.Sustainable conservation.Sustainable conservation can be achieved using a systematic conservation planning approach to target ecosystem service resilience, the key to both ecosystem service conservation and landscape sustainability.

Figure 1 .
Figure1.Sustainable conservation.Sustainable conservation can be achieved using a systematic conservation planning approach to target ecosystem service resilience, the key to both ecosystem service conservation and landscape sustainability.

19 Figure 2 .
Figure 2. Flow chart of the systematic landscape conservation planning (SLCP) process.Note: 1999a denotes benchmark data from March 1999.

Figure 2 .
Figure 2. Flow chart of the systematic landscape conservation planning (SLCP) process.Note: 1999a denotes benchmark data from March 1999.

Figure 3 .
Figure 3.The Chenyulan Watershed study area in central Taiwan.(a) elevation; (b) land use in March 1999. .
5).The future rainfall and temperature for each cell was then calculated from historical data and downscaled climate data."No-change" data indicates the historical average trend of the climate data (i.e., climate data used is the same as the historical climate data for the study period) in which no considerations of the impact of climate change were taken into account.Besides, the climate data from TCCIP represents the changing of climate from 2021 to 2040.Thus, we use 2040 as the focal year to predict the future land use scenarios.The climate data were also adopted as driving factors, leading to the change of land use, and input into CLUE-S model.For land use demands of each land use type, three scenarios, including calculated with Markov chain, calculated with linear regression and the same as land use demands in 2005, are considered since the RCPs include different socio-economic assumptions related to population, economic growth, and energy consumption.The scenarios considering land use change and climate are as follows: (1)

Figure 3 .
Figure 3.The Chenyulan Watershed study area in central Taiwan.(a) elevation; (b) land use in March 1999.
Note: 1999a refers to the March 1999 stand-alone period.For hotspot status change, the numbers '1' through'5'   indicate that the ecosystem service hotspot status underwent change (i.e., gained hotspot status or lost hotspot status) a total of 1, 2, 3, 4, or 5 times during the 1999-2005 period.Sustainability 2018, 10, x FOR PEER REVIEW 9 of 19

Figure 4 .
Figure 4. Landscape stability and resilience for ecosystem services and habitat quality hotspots.(a) water yield; (b) nitrogen retention; (c) phosphorus retention; (d) sediment retention; (e) carbon storage; (f) habitat quality.Note: red denotes hotspots and Gradient denotes elevation (m).

Figure 4 .
Figure 4. Landscape stability and resilience for ecosystem services and habitat quality hotspots.(a) water yield; (b) nitrogen retention; (c) phosphorus retention; (d) sediment retention; (e) carbon storage; (f) habitat quality.Note: red denotes hotspots and Gradient denotes elevation (m).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2071-1050/11/5/1393/s1.Author Contributions: The scope of this study was developed by Y.-P.L.The first manuscript draft was written by Y.-P.L., W.-Y.L., and J.R.P. was substantially revised by Y.-P.L., W.-Y.L., J.R.P., and L.-C.C.The modeling works have been done by C.-J.C. and W.-H.C. Funding: This research was funded by the Ministry of Science and Technology of the Republic of China, Taiwan, for financially supporting this research under contract numbers: 103-2410-H-002 -161 -MY3 and 105-2621-M-002 -003 -MY3.
Note: 1999a refers to the March 1999 stand-alone period; 1999b refers to the October 1999 stand-alone period.

Table 5 .
Historical and predicted land use during the study period (1999-2005) (in percentage of watershed, %).