Expansion of Protected Areas under Climate Change: An Example of Mountainous Tree Species in Taiwan

Abstract

Tree species in mountainous areas are expected to shift their distribution upward in elevation in response to climate change, calling for a potential redesign of existing protected areas. This study aims to predict whether or not the distributions of two high-mountain tree species, Abies (Abies kawakamii) and Tsuga (Tsuga chinensis var. formosana), will significantly shift upward due to temperature change, and whether current protected areas will be suitable for conserving these species. Future temperature change was projected for 15 different future scenarios produced from five global climate models. Shifts in Abies and Tsuga distributions were then predicted through the use of species distribution models (SDMs) which included occurrence data of Abies and Tsuga, as well as seasonal temperature, and elevation. The 25 km × 25 km downscaled General Circulation Model (GCMs) data for 2020–2039 produced by the Taiwan Climate Change Projection and Information Platform was adopted in this study. Habitat suitability in the study area was calculated using maximum entropy model under different climatic scenarios. A bootstrap method was applied to assess the parameter uncertainty of the maximum entropy model. In comparison to the baseline projection, we found that there are significant differences in suitable habitat distributions for Abies and Tsuga under seven of the 15 scenarios. The results suggest that mountainous ecosystems will be substantially impacted by climate change. We also found that the uncertainty originating from GCMs and the parameters of the SDM contribute most to the overall level of variability in species distributions. Finally, based on the uncertainty analysis and the shift in habitat suitability, we applied systematic conservation planning approaches to identify suitable areas to add to Taiwan’s protected area network.

1. Introduction

Forest tree species generally have long life spans and this makes them particularly vulnerable to climate change [1]. The ranges of vegetation in boreal, temperate, and tropical forest ecosystems are shifting in latitude and elevation in response to climate change [2,3,4,5]. This is due to the intrinsic physiological tolerances plant species have for various climatic and environmental variables as well as biotic factors, such as interspecies competition and parasite prevalence [6]. In addition, plant reproductive success can be affected by different climate conditions during a variety of developmental stages. Therefore, climate change could have substantial influence on their reproductive cycles [7]. Rapid changes in temperature and other climate parameters, such as changes in precipitation at high elevations, have significant impacts on plant communities [8,9]. For instance, the biological impacts of global warming have already been observed in the Alps; these include: tree line shifts to higher altitudes [10] and the reduced ranges of alpine plant species as their distributions shift to ever higher altitudes [9,11,12,13]. However, not all species tend to shift to higher elevations as expected; in fact, they may respond individualistically to changes in the environment [14,15,16,17]. Hence, by comparing differences in the direction and magnitude of shifting distributions between species, we can come to a better understanding of the exact forces and processes behind species distribution dynamics [18] and the variability of predictions made about future distributions.

Species distribution models (SDMs) are statistical algorithms used to predict the potential distribution of a species using presence-only or presence-absence records of a species and their corresponding independent variables which describe climatic and environmental factors [19,20]. SDMs are a widely used approach to evaluate the impacts of climate change on species and communities [21,22,23,24,25]. SDMs are essential tools for both predicting the impacts of climatic change on biodiversity, and designing conservation plans and policies which best mitigate the effects of climate change [26]. Recent notable examples are by [19,20,27,28,29,30,31,32,33]. However, it is recognized that there are significant limitations in projecting future species distribution using SDMs [34,35].

Most protected area networks have not yet taken the effects of future climate change into account [36]. However, if climate change progresses as projected, current protected areas may be inadequate to ensure the long-term persistence of certain species [36]. It is important to modify our biodiversity protection strategies in order to mitigate the impacts of climate change, especially in mountainous areas which are expected to be greatly impacted by climate change induced range shifts. It is incumbent on policy makers to take projections of future species distributions into account, and systematic conservation planning approaches [37] are increasingly being used to support these decisions [33].

In practice, inherent uncertainty due to the variability of both SDMs and projected future climatic models under different emission scenarios [38,39,40,41] is a particular challenge when considering the reconfiguration of conservation area networks. The variability arising from four different sources, namely: input data, SDMs, General Circulation Models (GCMs) and Greenhouse Gas Emissions Scenarios (GESs) [41,42], is critical in assessing future species distribution uncertainty. According to the study of Buisson et al. (2010), the main contributing factors to variation in projections of future species distributions are SDMs and GCMs [41]. A bootstrapping method, used to estimate the parameter uncertainty, climate uncertainty, and model uncertainty [38], provides an easy-to-implement tool to assess uncertainty without the assumption of a statistical distribution [43]. Selle and Hannah (2010) proposed a bootstrap approach to evaluate parameter uncertainty in dynamic catchment models [43]. To fully incorporate uncertainty in the process, Guilhaumon et al. (2008) used the nonparametric bootstrapping procedure to calculate the confidence interval of the parameters in nonlinear regression models for predicting bird species richness [44]. In this study, the variability in species distributions caused by the three different uncertainty sources, including SDMs, GCMs, and GESs were assessed using the bootstrap resampling technique. The relative proportion of each source of uncertainty that contributed to the overall variability of the final species distribution models was then quantified using a general linear model.

Here we project the impacts of climate change on the distribution of two high-mountain tree species of conservation importance in Taiwan: Abies (Abies kawakamii) and Tsuga (Tsuga chinensis var. formosana). The presence-only dataset of Abies and Tsuga are used to identify correlations to their corresponding background climatic and environmental variables. We use the trained maximum entropy model and future climate data to forecast the species future distributions. We then use the bootstrap resampling technique to analyze the variability from our dataset and from the parameters in the maximum entropy model, GCMs, and GESs. The Wilcoxon nonparametric statistical test is applied to check for significant differences in elevation shifts of the Abies and Tsuga ranges under future climate scenarios and the baseline scenario. We conclude by assessing the effectiveness of current protected areas under different climate change scenarios and apply the Zonation conservation planning software to identify additional conservation areas which provide refugia for the modeled high-mountain species faced with climate change induced range shifts.

2. Materials and Method

2.1. Study Area and Sampling Data

Taiwan is a mountainous island in eastern Asia. It is bisected by the Tropic of Cancer, which divides the island’s climate into tropical monsoon in the south and sub-tropical monsoon in the north. Approximately 60% of the island is comprised of mountainous regions, including the north-south Central Range, the Hsueh Shan Range, and the Yu Shan Range, which contain more than 200 summits over 3000 m in elevation. Due to the complicated geology and mountainous environment, there is a high degree of diversity in vegetation. The vegetation of Taiwan can be divided into six zones [45] including the alpine zone, the Abies zone, the Tsuga-Picea zone, the Quercus zone, the Machilus-Castanopsis zone, and the Ficus-Machilus zone, based on the dominant species and climate at different altitudes. The coniferous forests are divided into three zones [45]:

(1)

Alpine vegetation zone (mainly composed of Juniperus squamata forests and scrubs);

(2)

Subalpine coniferous forests (Abies zone);

(3)

Upper montane coniferous forests (Tsuga-Picea zone).

The high-mountain coniferous forest zones are dominated by species such as Abies kawakamii, Tsuga chinensis var. formosana, Picea morrisonicola and Pinus armandii var. mastersiana. Due to the high prevalence of Abies (Abies kawakamii) and Tsuga (Tsuga chinensis var. formosana), these two species were selected as the target species in this study. The Abies species, which is endemic to Taiwan, has an altitudinal range between 3000 m and 3500 m above sea level [46] with a cold and humid environment [47]. Abies, one of the 40 fir species known worldwide, is also endangered. This is possibly due to its specific survival needs, which restrict it to a narrow habitat. This species has received little conservation attention [46]. The ecotone of Abies and Tsuga exists at the lower elevations around 3000 m. The Tsuga is widely distributed on moist north facing slopes at 2500–3100 m elevation; and it is mixed with Chamaecyparis formosensis and C. obtusa var. formosana in the prevalent cloud forests around 1800–2400 m. In the northern part of the island, which is affected by the northeastern monsoon, the lower limit of Tsuga is approximately 1500 m.

We compiled the raw dataset used in this study from the Third Forest Resources and Land Use Inventory in Taiwan (TFRLUI) published by the Forestry Bureau in the Council of Agriculture [48]. The inventory data, from 1990 to 1995, included land-use types, forest types, timber, animal resources and soil types. The investigation of each forest plot involved systematic sampling every three kilometers from the starting point at 121.5135° E, 25.03534° N. A total of 3996 ground plots were selected from aerial photos of the national forests. For this study, vegetation raster cells, currently classified as sustaining an Abies and/or a Tsuga population, were chosen for modeling. All forests’ survey data (Figure 1c,d) were used for predicting spatial distributions based on various climate scenarios from different GCMs.

In this study, the five GCMs, CSIRO-Mk3.5 (CS), GFDL-CM2.1 (GF), MIROC3.2 (medres) (MI), MPI-ECHAM5 (MP), and MRI-CGCM2.3.2a (MR), and three greenhouse gas emissions scenarios, B1, A1B, and A2, were adopted to represent different temperature sensitivities for the period 2020–2039 [49]. Unfortunately, GCMs output data has a horizontal resolution between 250 and 600 km; since SDMs generally require a grid resolutions on the order of 1–5 km² [50], the direct application of SDMs on GCM output data is rarely possible. The original GCMs output data for 2020–2039 was downscaled to the resolution of 25 km × 25 km by the Taiwan Climate Change Projection and Information Platform (TCCIP), and then was downscaled to the resolution of 3 km ×3 km (Figure 1f) by Area-to-Point (ATP) co-kriging and validated through the use of weather station networks from the Central Weather Bureau (Figure 1b) and a high resolution digital elevation model (Figure 1e) [51,52].

Figure 1. (a) The topography of Taiwan; (b) locations of weather stations; (c) sample sites (presence-only) for Abies; (d) sample sites (presence-only) for Tsuga; (e) elevation data from a digital elevation model; (f) downscaled annual temperature data for 2020–2039.

Figure 1. (a) The topography of Taiwan; (b) locations of weather stations; (c) sample sites (presence-only) for Abies; (d) sample sites (presence-only) for Tsuga; (e) elevation data from a digital elevation model; (f) downscaled annual temperature data for 2020–2039.

2.2. Maximum Entropy

The distributions of the focal species were estimated by maximum entropy based on the site where the species were recorded and the correlating background environmental variables (driving factors). A logistic regression analysis using SPSS (Statistical Package for the Social Sciences, SPSS Inc., Chicago, IL, USA) identified the environmental variables with greatest influence on distribution, which were then used in the SDMs. The candidate variables included: moisture, precipitation, slope, aspect, elevation, and temperature. However only four significant environmental variables were identified by the logistic regression, and subsequently used in this study: average annual temperature, average spring monthly temperature, average summer monthly temperature, and elevation. Though precipitation is often a crucial factor determining tree species distribution, for mountainous regions of Taiwan, where annual precipitation often exceeds 3000 mm, from the macroscopic prospective this may not be a limiting factor in most regions of our study area. Aspect is also an environmental factor which is often crucial to tree species distribution. However, in relatively low latitude regions, such as Taiwan, aspect may not be as important compared to high latitude regions, especially when considering the humid conditions intrinsic to Taiwan. Our field experience substantiates the findings of our logistic regression model because distributions of both Abies kawakamii and Tsuga chinensis var. formosana were observed across large ranges which span both steep slopes and valleys, indicating that their distribution is predominantly determined by macro-climate environmental factors such as temperature.

Following Dudík et al. (2007), an unknown probability distribution π(x) over a space X across the entire study area was estimated, and the presence records were used to determine a set of constraints that were likely to be satisfied by π(x) [53]. The focal species were only present at certain locations in the set of location samples ($x_1,\dots ,x_{\text{m}}$), each of which corresponds to specific environmental variables. The expected value $\pi [z]$ of $z$, selected environment variables, under the distribution $\pi$ is defined as:

$$\pi \left[z\right]={\displaystyle \sum }_{x\in X}\pi \left(x\right)z(x)$$

(1)

where $X$ is the entire study area.

Since the empirical average of z ($\tilde{\pi }\left[z\right]$) and the expectation of z ($p\left[z\right]$) are both very close to its true expectation ($\tilde{\pi }[z]\approx \pi [z]$ and $p[z]\approx \pi [z]$) [54], many distributions $p(x)$ satisfy the constraint $p[z]\approx \pi [z]$ however the $p(x)$ whose entropy is largest should be chosen. The entropy of $p(x)$ is defined as follows [54]:

Maximize

$$H\left(p\right)={\displaystyle \sum }_{x\in X}p\left(x\right)lnp(x)$$

(2)

Subject to

$$\left|p\left[f\right]-\tilde{\pi }\left[f\right]\right|\le \beta$$

(3)

where β is the estimated upper bound (tolerance) of how close $p\left[f\right]$ is to $\tilde{\pi }\left[z\right]$ (the empirical distribution of the focal species). The $p\left(x\right)$ can be written as the following function (Gibbs model):

$$p\left(x\right)=\frac{\text{exp}({{\displaystyle \sum }}_{j=1}^n{\lambda }_j{z}_j(x))}{{{\displaystyle \sum }}_{x\in X}\text{exp}({{\displaystyle \sum }}_{j=1}^n{\lambda }_j{z}_j(x))}$$

(4)

where ${\lambda }_j$ is the coefficient of jth environmental variables (features) $z_j$.

2.3. The Bootstrapping Method

2.3.1. Evaluating Data and Model Performance

The bootstrapping method was implemented by randomly selecting 80% of the presence records from the survey data used to estimate the distribution of focal species by the maximum entropy method. The remaining 20% of the presence data were used as a test dataset and applied to validate the performance by computing the area under the receiver operating characteristic curve (AUC) value. This procedure was done iteratively creating 1000 maximum entropy species distribution models and corresponding validation AUC values.

2.3.2. The Confidence Region of SDM Parameters

The ellipsoidal confidence regions, which are generalizations of univariate symmetrical confidence intervals, for the coefficients of the environmental variables in the maximum entropy model, ${\lambda }_j$, are generated by using the nonparametric bootstrapping procedure [55]. The confidence regions can be written as the following equation [55]:

$$R=\left\{\hat{\theta }-\sqrt{\hat{\Omega }}\times r:r\in R_{\alpha }\right\}$$

(5)

where $R_{\alpha }$ is any region such that $P\left[\sqrt{\hat{\Omega }}\times \left(\theta -\hat{\theta }\right)\in R_{\alpha }\right]=\alpha$ and α was set to 0.95. The shape of the region R depends on the shape of the region $R_{\text{α}}$, which is ellipsoid in this study.

2.4. Variability Analysis

Uncertainties intrinsic to future species distribution models based on future climate models are many and varied [38,39,40,41]. Although it has been demonstrated that SDMs contribute significantly to the total variability of future species distribution projections [22,60], the overall uncertainty in suitable future habitats arising from this variability has rarely been considered [56,57,58,59,60]. In order to enhance the reliability of projections, it is necessary to quantify the uncertainty from several sources. In this way, we can construct an ensemble of relatively reliable forecasts [61].

In this study, the SDM model variability arising from three different sources of uncertainty, including SDM parameters, GCMs and GESs was quantified. Evaluation of relative variability originating in SDM parameters was undertaken by randomly generating 1000 sets of coefficients via a sample presence bootstrapping technique. Coefficients falling within the 95% confidence range were then compared with those coefficients falling outside of the confidence range. The Wilcoxon test was adopted in this study to identify the significant differences in elevation of projected suitable habitats for the two focal species, under 15 future scenarios (5 GCM × 3 GES) and the baseline scenario, based on the parameters from 1000 sets of coefficients. This procedure led to a factorial design which crossed 1000 sets of the SDM’s coefficients, five GCMs and three GESs, thus resulting in 15,000 different projections of future distribution for each species. Based on this design, the proportion of the three different sources of uncertainty, GCM, GES, and SDM model’s parameters, were quantified using a general linear model, which was assessed by calculating the ratio between the deviances explained by each source against the null deviance [41].

Besides the sources of uncertainties, the local and spatial uncertainties were also quantified in this study. The local uncertainty at location k (${P_l}^k)$ can be defined by the following equation [62]:

$${P_l}^k=n_p(g_k)/1000$$

(6)

where 1000 is the number of realizations, and $n_p(g_k)$ is the number from the 1000 realizations in which the simulated species were present at location k.

The joint probability (Pj, global uncertainty) in m occurrence locations $\left(u_1,u_2,u_3,\cdots u_m\right)$ was written as follows [63,64]:

$$P_j=n_p\left(u_1,u_2,u_3,\cdots u_m\right)/1000$$

(7)

where 1000 is the number of realizations, and $n_p\left(u_1,u_2,u_3,\cdots u_m\right)$ is the number from the 1000 realizations in which the simulated species were present at m locations; m is the number of cells, out of 1000 realizations, in which the given critical probability, or cutoff point that could be viewed as the given minimal habitat suitability threshold, is satisfied; n_p is the number of realizations, out of 1000, whose projected focal species distributions are entirely within the given critical probability area.

2.5. Identifying Protected Areas

Zonation is a systematic conservation planning decision support tool. It generates a hierarchy of conservation priorities over all cells in a landscape based on their conservation value [65]. We used Zonation to delineate protected areas for the modeled species under the baseline and the climate change scenarios. The selections were made based on occurrences of conservation targets and pseudo targets, which represent the possible occurrences of target species following distribution predictions. In core-area Zonation, cell removal is done in a manner that minimizes biological loss by picking cell i that has the smallest δ:

$${\delta }_i=\underset{j}{Max}\frac{Q_{ij}(S)w_j}{c_i}$$

(8)

where w_j is the weight (or priority) of species j and ci is the cost of adding cell i to the reserve network [66]. The critical part of the equation is Q_ij(S), which is the proportion of the remaining distribution of species j located in cell i in the remaining set of cells, S. When a part of the distribution of a species is removed, the proportion located in each remaining cell increases.

We then computed an index which integrates both the estimated habitat suitability and conservation status, as indicated by Zonation, to identify cells of greatest suitability for Abies and Tsuga within a protected area network under 15 future scenarios and the baseline scenario. The matching of climate suitability with conservation area for species $r$, in a given grid cell $w$, can be expressed as ($S_w^r$) [67]:

$$S_w^r=H{S}_w^r\times C{A}_w$$

(9)

where $H{S}_w^r$ is the suitability score of grid cell $\text{w}$ for a given species $r$, and $C{A}_w$ is 1, if the grid cell $w$ was selected for protection; otherwise, $C{A}_w$ is 0. The $S_w^r$ ranges from 0, for grid cells which are unsuitable for species or are not selected for conservation, to 1, for grid cells which have both high suitability for the focal species and are also selected for conservation. The habitat suitability within protected areas for species $\text{ r}$ over the study area can be written as the following equation:

$$S^r={\displaystyle \sum }_{w=1}^N{S}_w^r$$

(10)

where $N$ is the number of grid cells in the study area. Moreover, the proportions of habitat suitability covered by protected areas for species $r$, ($C{P}^r$), can be expressed as:

$$C{P}^r=S^r/{\displaystyle \sum }_{w=1}^{N}H{S}_w^r$$

(11)

Large $C{P}^r$ values correspond to large protected area networks which possess a large number of high quality habitat cells.

3. Results and Discussion

3.1. SDM Performance, Focal Species Distribution and Their Shifts in Elevation

The validation AUC values calculated from 1000 iterations of our SDM range from 0.8 to 0.96 for Abies and from 0.85 to 0.95 for Tsuga. The mean and the standard deviation of the 1000 AUC values are 0.88 and 0.021 for Abies and 0.91 and 0.015 for Tsuga. The higher mean AUC value and lower standard deviation correlate with more precise maximum entropy predictions. The high accuracy of the predictions leads to higher consistencies between predicted distributions of focal species and presence data distributions of Abies (Figure 1b and Figure 2a) and Tsuga (Figure 1c and Figure 2b).

Figure 2. The predicted baseline habitat suitability for: (a) Abies; (b) Tsuga.

Figure 2. The predicted baseline habitat suitability for: (a) Abies; (b) Tsuga.

The trees species future distributions under different scenarios are presented in Figure 3 and Figure 4. The distributions are slightly different for the two species; however, the SDMs reveal that the species prefer to live in areas with low temperatures and high altitudes.

Figure 3. Predicted habitat suitability in Central Mountain Range area for Abies in 12 future scenarios (3 GES × 2 GCM × 2 parameter) for 2020–2039.

Figure 3. Predicted habitat suitability in Central Mountain Range area for Abies in 12 future scenarios (3 GES × 2 GCM × 2 parameter) for 2020–2039.

Figure 4. Predicted habitat suitability in Central Mountain Range area for Tsuga in 12 future scenarios (3 GES × 2 GCM × 2 parameter) for 2020–2039.

Figure 4. Predicted habitat suitability in Central Mountain Range area for Tsuga in 12 future scenarios (3 GES × 2 GCM × 2 parameter) for 2020–2039.

3.2. Variability Analysis

The discrepancies in results of species habitat shifts with constant GCMs and GESs were interpreted as the uncertainty introduced by maximum entropy model parameters. The distribution of coefficients of the environmental variables in the maximum entropy model for Abies and Tsuga, which are elevation, average annual temperatures, average spring monthly temperatures and average summer monthly temperatures, based on bootstrap resampling are shown in Figure 5. For Abies, they are in the range of 0.29 to 4.70, −3.42 to −1.19, −2.02 to 1.11, and −4.59 to −0.97, respectively. For Tsuga, they are in the range of −1.26 to 1.61, −3.29 to −1.98, −2.57 to −0.31, and −3.49 to −1.20, respectively. The dashed lines in Figure 5 show the SDM’s parameters fitted by the original data, which approximate the mean parameters of 1000 sets of bootstrap resampling data indicating the low bias of the estimated parameters. For evaluating the variability of parameters, 1000 sets of coefficients which were within the 95% bootstrap confidence interval were used to estimate the future distribution of the modeled species. The nonparametric Wilcoxon test was then used to identify significant differences in habitat elevation between the baseline scenario versus each future scenario (Table 1). The shifts in elevation were based on species current distribution, and the significant environmental variables which these distributions correspond to.

Figure 5. The bootstrap histogram of the coefficients in maximum entropy model for Abies: (a) elevation; (b) average annual temperature; (c) average monthly spring temperature; (d) average monthly summer temperature, and for Tsuga: (e) elevation; (f) average annual temperature; (g) average monthly spring temperature; (h) average monthly summer temperature. The blue line is the fitted kernel density.

Figure 5. The bootstrap histogram of the coefficients in maximum entropy model for Abies: (a) elevation; (b) average annual temperature; (c) average monthly spring temperature; (d) average monthly summer temperature, and for Tsuga: (e) elevation; (f) average annual temperature; (g) average monthly spring temperature; (h) average monthly summer temperature. The blue line is the fitted kernel density.

According to the results of predicted temperature data from different GCMs and GESs (Table 1), greater temperature increases correlate with a higher “count”, which represents the number of statistically significant differences in habitat elevation under baseline and future climate scenarios in 1000 SDM outputs. Moreover, a higher “count” indicates a greater reliability in predicting the upward shift in focal species distribution. The homogeneity of the “count” between the 15 scenarios was tested using the chi-square test, which contrasts the non-homogeneity of the “count” between the 15 scenarios (p value < 0.0001) for each species. The “count” under scenarios 2, 3, 6, 7, 8, 11 and 13 for Abies and scenarios 1, 2, 3, 6, 7, 8, 11 and 13 for Tsuga are significantly higher (+) than that under other scenarios (Table 1). In addition, scenario 8 has the largest “count” for Abies (668) and Tsuga (841). The highest average annual temperatures and average summer monthly temperatures in scenario 8 lead to the most drastic upward shift in distribution for both species. As for other scenarios, the temperature in scenario 2, 3, 6, 7, 11 and 13 all increase leading to a higher “count” for each species. Moreover, as can be seen from Table 1 as the “count” approaches extreme values, which are 0 or 1000, only slight variability of parameters in the maximum entropy model are observed. The scenarios with large variability from parameters in predicting the upward shift of the focal species distributions are scenarios 2, 6 and 13 for Abies and scenarios 1, 2 and 7 for Tsuga. This result indicates that the intermediate increments of temperature change lead to the highest variability in SDM parameters. In other words, we have higher confidence when predicting the upward shift of tree species under scenarios which project larger increments of temperature change.

Table 1. The difference in temperature under each climate change scenario (for 2020–2039) in comparison to the baseline scenario, and the number of statistically significant differences in the optimal elevation of focal species habitat in 1000 different sets of outputs, which are contained in the 95% bootstrap confidence interval.

Scenario	GES	GCM	Change of Average Temperature (°C)			Significant Shift in Elevation of Habitat
			Annual	Spring	Summer	Abies		Tsuga
						Count	Uncertainty	Count	Uncertainty
1	B1	CS	0.77	1.03	0.75	193 (−)	0.16	489 (+)	0.25
2		GF	0.88	0.76	0.93	466 (+)	0.25	571 (+)	0.25
3		MI	1.02	1.22	0.93	576 (+)	0.24	805 (+)	0.16
4		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
5		MR	0.34	0.17	0.59	0 (−)	0.00	2 (−)	0.00
6	A1B	CS	0.88	1.17	1.02	514 (+)	0.25	722 (+)	0.20
7		GF	0.84	0.52	0.82	341 (+)	0.23	442 (+)	0.25
8		MI	1.09	1.08	1.07	668 (+)	0.22	841 (+)	0.13
9		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
10		MR	0.38	0.33	0.67	3 (−)	0.00	29 (−)	0.03
11	A2	CS	0.78	1.13	0.89	355 (+)	0.23	589 (+)	0.24
12		GF	0.67	0.69	0.79	149 (−)	0.13	368 (0)	0.23
13		MI	0.91	1.29	0.91	489 (+)	0.25	757 (+)	0.18
14		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
15		MR	0.35	0.31	0.51	0 (−)	0.00	1 (−)	0.00

Note: The symbols (+), (0) and (−) represent the “count” as being significantly higher than, equal to, and significantly lower than expected based on the chi-square test; GES: Greenhouse Gas Emissions Scenarios; GCM: General Circulation Model; CS:CSIRO-Mk3.5; GF: GFDL-CM2.1; MI : MIROC3.2 (medres); MP : MPI-ECHAM5 (MP); MR : MRI-CGCM2.3.2a.

Table 1. The difference in temperature under each climate change scenario (for 2020–2039) in comparison to the baseline scenario, and the number of statistically significant differences in the optimal elevation of focal species habitat in 1000 different sets of outputs, which are contained in the 95% bootstrap confidence interval.

Scenario	GES	GCM	Change of Average Temperature (°C)			Significant Shift in Elevation of Habitat
			Annual	Spring	Summer	Abies		Tsuga
						Count	Uncertainty	Count	Uncertainty
1	B1	CS	0.77	1.03	0.75	193 (−)	0.16	489 (+)	0.25
2		GF	0.88	0.76	0.93	466 (+)	0.25	571 (+)	0.25
3		MI	1.02	1.22	0.93	576 (+)	0.24	805 (+)	0.16
4		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
5		MR	0.34	0.17	0.59	0 (−)	0.00	2 (−)	0.00
6	A1B	CS	0.88	1.17	1.02	514 (+)	0.25	722 (+)	0.20
7		GF	0.84	0.52	0.82	341 (+)	0.23	442 (+)	0.25
8		MI	1.09	1.08	1.07	668 (+)	0.22	841 (+)	0.13
9		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
10		MR	0.38	0.33	0.67	3 (−)	0.00	29 (−)	0.03
11	A2	CS	0.78	1.13	0.89	355 (+)	0.23	589 (+)	0.24
12		GF	0.67	0.69	0.79	149 (−)	0.13	368 (0)	0.23
13		MI	0.91	1.29	0.91	489 (+)	0.25	757 (+)	0.18
14		MP	−0.03	−0.06	−0.2	0 (−)	0.00	0 (−)	0.00
15		MR	0.35	0.31	0.51	0 (−)	0.00	1 (−)	0.00

Note: The symbols (+), (0) and (−) represent the “count” as being significantly higher than, equal to, and significantly lower than expected based on the chi-square test; GES: Greenhouse Gas Emissions Scenarios; GCM: General Circulation Model; CS:CSIRO-Mk3.5; GF: GFDL-CM2.1; MI : MIROC3.2 (medres); MP : MPI-ECHAM5 (MP); MR : MRI-CGCM2.3.2a.

The variability in projecting future species distributions does not only originate from SDM parameters but also arises from GCMs and GESs [59,61]. The outputs of all GCM × GES combinations cover a wide range of potential future climate conditions; thus, the need to evaluate the variability from both GCMs and GESs is increasing [41]. The variability arising from the three uncertainty sources analyzed in this study, GCMs, GESs, and bootstrapped parameters, are listed in Table 2. The proportions of variability are estimated by calculating the ratio between the deviance explained by one source and the null deviance with the generalized linear model based on 15,000 different future projections for each species [41]. Due to the high variation of increasing temperature between each GCM based on the same GES, the contributions of GCMs, which are 31.3% for Abies and 42.8% for Tsuga, are significantly higher than that of GESs, which are 0.9% for Abies and 0.2% for Tsuga. This result is consistent with that of the “count” in five GCMs based on the same GES which are significantly different from one another in Table 1. This implies that a high level of variability between different GCMs exists, and that this variability affects predictions of future species habitat distributions.

Table 2. Proportion of variability from three different uncertainty sources, including model parameters, GCMs and GESs.

	Source of Uncertainty	Contribution (%)
	parameters	67.9
Abies	GCMs	31.3
	GESs	0.9
	parameters	57
Tsuga	GCMs	42.8
	GESs	0.2

Table 2. Proportion of variability from three different uncertainty sources, including model parameters, GCMs and GESs.

	Source of Uncertainty	Contribution (%)
	parameters	67.9
Abies	GCMs	31.3
	GESs	0.9
	parameters	57
Tsuga	GCMs	42.8
	GESs	0.2

In contrast, the proportions of variability originating from GESs are small because, all else being equal, the predicted temperatures under the three GESs are similar. This result corroborates the studies by both Stott & Kettleborough (2002) and Buisson et al. (2010), who found that the variability from different GESs is small [41,68]. On the other hand, the proportion of variability in predicting upward shifts of tree species originating from SDM parameters is quite large, being 67.9% for Abies and 57.0% for Tsuga. This means that the variability from the model's parameters has a strong effect on predicting the shift of focal species. This result is in accordance with the study by Buisson et al. (2010), who suggested that the high proportions of the total variability in future projections are from SDMs [41]. In summary, the above findings suggest that variability originating from both SDM parameters and GCMs should be of first concern as they contribute most to the total variability in predicting the upward shift of focal species.

3.3. Local and Global Uncertainties

Figure 6 shows the proportion of habitat in 1000 realizations under the baseline scenario for Abies and Tsuga, which provides a measure of local uncertainty associated with a single location (grid). Under each scenario, the areas with high (95%) habitat coverage were located in the center of Taiwan where the central mountain range is located. As for the global uncertainty, given a critical probability of species occurrence, a higher joint probability means greater reliability (lower global uncertainty) of the mapped locations of species occurrences. The joint probabilities increase with the critical probability under each scenario (Table 3). Kerry et al. (2010) and Goovaerts (2009) suggest the critical probabilities 75% and 90% are enough to generate a reliable joint probability. In this study, the joint probabilities under each scenario are larger than 0.5 when the critical probability reaches 0.75 [63,69].

Table 3. The global uncertainty under baseline and 15 future scenarios.

	Cut-off	Abies						Tsuga
Scenario		0.7	0.75	0.8	0.85	0.9	0.95	0.7	0.75	0.8	0.85	0.9	0.95
B		0.43	0.58	0.65	0.7	0.78	0.9	0.43	0.58	0.65	0.7	0.78	0.9
s1		0.49	0.56	0.62	0.72	0.82	0.89	0.49	0.56	0.62	0.72	0.82	0.89
s2		0.45	0.54	0.65	0.73	0.82	0.91	0.45	0.54	0.65	0.73	0.82	0.91
s3		0.47	0.62	0.66	0.73	0.83	0.91	0.47	0.62	0.66	0.73	0.83	0.91
s4		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s5		0.46	0.54	0.66	0.71	0.8	0.92	0.46	0.54	0.66	0.71	0.8	0.92
s6		0.45	0.51	0.67	0.72	0.84	0.91	0.45	0.51	0.67	0.72	0.84	0.91
s7		0.44	0.53	0.59	0.77	0.81	0.89	0.44	0.53	0.59	0.77	0.81	0.89
s8		0.55	0.6	0.67	0.74	0.81	0.89	0.55	0.6	0.67	0.74	0.81	0.89
s9		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s10		0.48	0.54	0.62	0.72	0.81	0.91	0.48	0.54	0.62	0.72	0.81	0.91
s11		0.5	0.53	0.66	0.75	0.82	0.9	0.5	0.53	0.66	0.75	0.82	0.9
s12		0.47	0.53	0.62	0.74	0.86	0.91	0.47	0.53	0.62	0.74	0.86	0.91
s13		0.45	0.52	0.68	0.74	0.83	0.91	0.45	0.52	0.68	0.74	0.83	0.91
s14		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s15		0.46	0.56	0.64	0.7	0.82	0.9	0.46	0.56	0.64	0.7	0.82	0.9

Table 3. The global uncertainty under baseline and 15 future scenarios.

	Cut-off	Abies						Tsuga
Scenario		0.7	0.75	0.8	0.85	0.9	0.95	0.7	0.75	0.8	0.85	0.9	0.95
B		0.43	0.58	0.65	0.7	0.78	0.9	0.43	0.58	0.65	0.7	0.78	0.9
s1		0.49	0.56	0.62	0.72	0.82	0.89	0.49	0.56	0.62	0.72	0.82	0.89
s2		0.45	0.54	0.65	0.73	0.82	0.91	0.45	0.54	0.65	0.73	0.82	0.91
s3		0.47	0.62	0.66	0.73	0.83	0.91	0.47	0.62	0.66	0.73	0.83	0.91
s4		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s5		0.46	0.54	0.66	0.71	0.8	0.92	0.46	0.54	0.66	0.71	0.8	0.92
s6		0.45	0.51	0.67	0.72	0.84	0.91	0.45	0.51	0.67	0.72	0.84	0.91
s7		0.44	0.53	0.59	0.77	0.81	0.89	0.44	0.53	0.59	0.77	0.81	0.89
s8		0.55	0.6	0.67	0.74	0.81	0.89	0.55	0.6	0.67	0.74	0.81	0.89
s9		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s10		0.48	0.54	0.62	0.72	0.81	0.91	0.48	0.54	0.62	0.72	0.81	0.91
s11		0.5	0.53	0.66	0.75	0.82	0.9	0.5	0.53	0.66	0.75	0.82	0.9
s12		0.47	0.53	0.62	0.74	0.86	0.91	0.47	0.53	0.62	0.74	0.86	0.91
s13		0.45	0.52	0.68	0.74	0.83	0.91	0.45	0.52	0.68	0.74	0.83	0.91
s14		0.42	0.52	0.69	0.74	0.79	0.88	0.42	0.52	0.69	0.74	0.79	0.88
s15		0.46	0.56	0.64	0.7	0.82	0.9	0.46	0.56	0.64	0.7	0.82	0.9

Figure 6. The proportion of habitat in 1000 realizations (local uncertainty) under the baseline scenario for: (a) Abies; (b) Tsuga.

Figure 6. The proportion of habitat in 1000 realizations (local uncertainty) under the baseline scenario for: (a) Abies; (b) Tsuga.

3.4. Implications for Protected Areas

Only a limited number of conservation planning studies have taken the effects of climate change into account. However, climate change can induce rapid shifts of species distribution resulting in existing protected areas no longer ensuring the long-term persistence of certain species [36,67]. In this study, for Abies and Tsuga, the climate suitability (S) and the climate suitability within protected areas (Scons), including Taiwan National Parks, Taiwan Wildlife Refuges, and Taiwan Major Wildlife Habitats are likely to decrease under most future scenarios (Table 4). However, the proportions of suitable habitat protected for Abies and Tsuga under the baseline and 15 future climate scenarios are very similar to each other, which are 87.2%–90.9% and 85.3%–88.7% for Abies and Tsuga, respectively. This result is not consistent with previous studies’ concern that climate change may lead to the failure of protected areas to cater to the needs of focal species [36,67,70,71,72]. In our study area, despite the shrinkage of focal species habitat under many of the future scenarios, the existing conservation areas are still adequate. In fact, the proportions of species habitat protected under many of the future climate scenarios increase in comparison to the baseline scenario as reduced total area of suitable habitat force species into conservation areas.

Table 4. Total habitat suitability (S) and habitat suitability within protected areas (Scons) and the proportion of suitable habitat under existing protected areas (CP) under baseline and 15 future climate scenarios (for 2020–2039) for Abies and Tsuga.

Scenario	S	Abies Scons	CP (%)	S	Tsuga Scons	CP (%)
B	199.7	174.2	87.2	295.3	252.8	85.6
145.4	161.6		90	226.2	199.6	88.2
s2	155.6	140.6	90.4	221.9	196.3	88.5
s3	152.8	138.4	90.6	211.7	187.7	88.7
s4	207.0	179.7	86.8	305.5	260.7	85.3
s5	175.5	155.8	88.8	258.7	223.6	86.4
s6	154.3	139.6	90.5	215.4	190.9	88.6
s7	159.5	143.8	90.1	230.2	202.5	88
s8	149.6	136.0	90.9	209.2	185.5	88.7
s9	207.0	179.7	86.8	305.5	260.7	85.3
s10	172.3	153.5	89	252.0	218.3	86.7
s11	158.2	142.7	90.2	220.8	195.4	88.5
s12	163.0	146.4	89.8	233.3	204.8	87.8
s13	155.1	140.2	90.4	214.1	189.8	88.6
s14	207.0	179.7	86.8	305.5	260.7	85.3
s15	176.7	156.7	88.7	257.7	222.8	86.5

Table 4. Total habitat suitability (S) and habitat suitability within protected areas (Scons) and the proportion of suitable habitat under existing protected areas (CP) under baseline and 15 future climate scenarios (for 2020–2039) for Abies and Tsuga.

Scenario	S	Abies Scons	CP (%)	S	Tsuga Scons	CP (%)
B	199.7	174.2	87.2	295.3	252.8	85.6
145.4	161.6		90	226.2	199.6	88.2
s2	155.6	140.6	90.4	221.9	196.3	88.5
s3	152.8	138.4	90.6	211.7	187.7	88.7
s4	207.0	179.7	86.8	305.5	260.7	85.3
s5	175.5	155.8	88.8	258.7	223.6	86.4
s6	154.3	139.6	90.5	215.4	190.9	88.6
s7	159.5	143.8	90.1	230.2	202.5	88
s8	149.6	136.0	90.9	209.2	185.5	88.7
s9	207.0	179.7	86.8	305.5	260.7	85.3
s10	172.3	153.5	89	252.0	218.3	86.7
s11	158.2	142.7	90.2	220.8	195.4	88.5
s12	163.0	146.4	89.8	233.3	204.8	87.8
s13	155.1	140.2	90.4	214.1	189.8	88.6
s14	207.0	179.7	86.8	305.5	260.7	85.3
s15	176.7	156.7	88.7	257.7	222.8	86.5

Nonetheless, new protected areas can play an important role in protecting Abies and Tsuga under climate change. In most of our climate change scenarios the proportion of suitable habitat covered by existing protected areas is higher than the baseline scenario because the total area of future suitable habitat declines. As such, increasing protected areas as determined by Zonation will be necessary to offset the loss of suitable habitat [71] (Figure 7). The additional area to be protected depends on the level of climate change. Therefore the union of all areas suitable under the 15 climate change scenarios and not currently protected would be optimal in establishing the most robust conservation area (Figure 7b). Interestingly, when considering the worst case scenario, that is the one which projects the largest increment of temperature change, a smaller protected area is identified in comparison to that correlating to all 15 climate change scenarios (Figure 7c).

Figure 7. Current protected areas in Taiwan, including Wildlife Refuges, National Parks and Major Wildlife Habitats, and the additional areas determined by Zonation which are currently outside of existing protected areas: (a) under baseline; (b) all future climate scenarios (the extra conserved area under future climate scenario is the union of that under 15 future scenarios); and (c) future climate scenario (the worst scenario with the largest increase of temperature). Note: * The difference in additional conservation area between baseline and future climate scenarios. The period of future climate scenario is for 2020–2039.

Figure 7. Current protected areas in Taiwan, including Wildlife Refuges, National Parks and Major Wildlife Habitats, and the additional areas determined by Zonation which are currently outside of existing protected areas: (a) under baseline; (b) all future climate scenarios (the extra conserved area under future climate scenario is the union of that under 15 future scenarios); and (c) future climate scenario (the worst scenario with the largest increase of temperature). Note: * The difference in additional conservation area between baseline and future climate scenarios. The period of future climate scenario is for 2020–2039.

4. Conclusions

The upward shift of mountainous and alpine habitats is expected under climate change. In Taiwan, an island that is 60% mountainous, the high elevation habitats will be directly influenced by increasing temperatures. In this study, the ensemble forecasting framework used to predict shifts in elevation of Abies and Tsuga habitats allowed for more robust decisions in conservation planning by taking multiple sources of variability into account, namely: SDM parameters, GCMs and GESs. The findings indicate that the projected shift in elevation of each species habitat differs among climate change models and among 1000 SDM realizations based on different sets of parameters. To a lesser degree, variability in predictions originating from different emission scenarios was also found to contribute to overall uncertainty. Our results show that most of the variability originates from the model’s parameters as well as the GCMs, indicating that uncertainty has a strong influence on predicting the shift of the modeled species distribution. Although most of the climate change induced shrinkage of habitats is already covered by existing protected areas, in order to offset the loss of suitable habitat and reduce uncertainties caused by climate change, we suggest additional habitat be added to the network of protected areas.