Uncertainty of future projections of species distributions in mountainous regions

Multiple factors introduce uncertainty into projections of species distributions under climate change. The uncertainty introduced by the choice of baseline climate information used to calibrate a species distribution model and to downscale global climate model (GCM) simulations to a finer spatial resolution is a particular concern for mountainous regions, as the spatial resolution of climate observing networks is often insufficient to detect the steep climatic gradients in these areas. Using the maximum entropy (MaxEnt) modeling framework together with occurrence data on 21 understory bamboo species distributed across the mountainous geographic range of the Giant Panda, we examined the differences in projected species distributions obtained from two contrasting sources of baseline climate information, one derived from spatial interpolation of coarse-scale station observations and the other derived from fine-spatial resolution satellite measurements. For each bamboo species, the MaxEnt model was calibrated separately for the two datasets and applied to 17 GCM simulations downscaled using the delta method. Greater differences in the projected spatial distributions of the bamboo species were observed for the models calibrated using the different baseline datasets than between the different downscaled GCM simulations for the same calibration. In terms of the projected future climatically-suitable area by species, quantification using a multi-factor analysis of variance suggested that the sum of the variance explained by the baseline climate dataset used for model calibration and the interaction between the baseline climate data and the GCM simulation via downscaling accounted for, on average, 40% of the total variation among the future projections. Our analyses illustrate that the combined use of gridded datasets developed from station observations and satellite measurements can help estimate the uncertainty introduced by the choice of baseline climate information to the projected changes in species distribution.


Introduction
Uncertainty is an important consideration for all climate change assessments.Ignoring or minimizing the importance of uncertainty can negatively affect the usefulness of assessment outcomes for decision making and planning [1].Uncertainty is a particular concern for climate change assessments of future species distributions, as sensitivity analyses have identified multiple sources of uncertainty that can substantially impact assessment findings.These include the availability and quality of species information (e.g., [2][3]), methodologies used to develop species distribution models (e.g., [4][5][6][7][8]), selection of predictor variables that capture environmental influences on species distributions (e.g., [9][10][11]), thresholds used to convert likelihood of occurrence to binary predictions of species presence (e.g., [12][13]), parameterizations and tuning of model settings (e.g., [14]), and choice of future climate simulations (e.g., [15][16][17]).Although few assessments explicitly consider all these uncertainty sources, often due to resource constraints, their importance when interpreting and applying the assessment findings for conservation planning is well established in the literature.
A less well understood uncertainty source is the choice of baseline climate information.As pointed out by Perdinan and Winkler [18], climate observations are the "backbone" of any climate change assessment.Baseline climate information is used in ecological assessments to calibrate species distribution models and, for many assessments, to downscale simulations from climate models to a finer spatial resolution [19].A challenge is that the spatial resolution of climate observing networks, ranging from tens of kilometers in developed regions to hundreds of kilometers in remote areas or at high elevations [20][21], may be insufficient to capture the local and regional climatic gradients that influence the distribution of a particular species [22].
Because of the limitations of climate observing networks, gridded baseline climate layers that have been generated at a fine resolution through spatial interpolation of station observations are frequently used in climate change assessments of future species distributions.Gridded climate datasets differ in terms of their spatial resolution (a few meters to tens of kilometers), temporal resolution (from sub-daily to long-term means), spatial extent (regional to global), and the complexity of the spatial interpolation technique [23].No "best" gridded dataset for ecological assessments exists [24], and users need be aware of the strengths and weaknesses of the available climate datasets.Moreover, a recent sensitivity analysis found that the uncertainty introduced into projected future distributions of African bird species by the differences among several popular gridded climate datasets in their depiction of regional to subcontinental climatic gradients was often larger than the uncertainty introduced by the choice of future climate projection [19].
The potential contribution of the choice of baseline climate information to assessment uncertainty is arguably larger for climate change assessments in mountainous regions, where complex local and regional climate gradients arise from climatic controls such as elevation, aspect, rain shadows, and cold air drainage [25].Datasets developed from coarse-resolution climate networks are unlikely to capture these fine-scale climatic gradients [26][27].For most mountainous regions, especially in more remote locations, high-density climate networks are unavailable, and those climate stations that do exist are primarily located at lower elevations where population densities are higher [28].Furthermore, time-dependent biases in the observational record of a climate station, introduced by changes in location, instrumentation, time of observation, and the surrounding environment [29], can have a large influence on gridded climate datasets when the number of stations is small relative to the scale of important climatic gradients.Consequently, the uncertainty introduced by baseline climate information to assessments of future species distributions, and its relative contribution compared to other uncertainty sources, is largely unknown for mountainous regions.Yet mountainous regions are home to rich biodiversity and account for approximately a quarter of the global land area [30].
We propose that the joint application of conventional and remotely-sensed gridded climate datasets can provide needed insights into the relative magnitude of the uncertainty contributed by the choice of baseline climate information to future projections of species distributions in mountainous regions.Remotely-sensed measurements are estimates of climate parameters obtained from radiative fluxes at various electromagnetic wavelengths [31] or by blending remotely-sensed estimates with observations from climate observing networks [32].Compared to conventional observations, remotely-sensed datasets have a more uniform spatial coverage and finer native spatial resolutions ranging from a few meters to a few kilometers.These datasets have their own limitations, though, such as atmospheric interference (e.g., clouds), accuracy deviations of the on-orbit sensor calibrations, and limitations of the algorithms used to estimate climate parameters from the radiation measurements, among others [33].In addition, remotely-sensed datasets usually have shorter record lengths than climate station observations.Therefore, remotely-sensed datasets are not a replacement for conventional datasets but potentially a complementary source of baseline climate information.
Below we use two sources of gridded climate data to evaluate the relative contribution of the choice of baseline climate information to the uncertainty of future species distributions in a mountainous environment.The first source is the popular WorldClim dataset (hereafter referred to as "WC") [34] derived from spatial interpolation of conventional climate station observations.The second source is a newly compiled dataset of remotely-sensed measurements of temperature and precipitation (hereafter referred to as "RS") [35].Specifically, we investigate: 1) the variance structure of the alternative sources of baseline climate information; 2) the influence of the baseline climate layers on the calibration and interpretation of species distribution models; 3) spatial variations in the projected future probabilities of species presence obtained from the alternative baseline datasets; and 4) the relative magnitude of the uncertainty introduced by the choice of baseline climate information into projections of future climatically-suitable area including the interaction with other sources of uncertainty.
The example study region is the current geographic range of the Giant Panda (Ailuropoda melanoleuca) which covers six mountain regions (i.e., the Qinling, Minshan, Qionglaishan, Daxiangling, Xiaoxiangling and Liangshan Mountains) (Fig 1) in southwest China [36].Twenty-one understory bamboo species consumed by the Giant Panda serve as the focal species.The Giant Panda is a global icon of biodiversity conservation and has been the subject of extensive conservation efforts [37][38].Future changes in the distribution of understory bamboo species due to climate change are anticipated to impact the extent and degree of fragmentation of panda habitat and the connectivity among panda subpopulations [39][40].However, the uncertainty introduced by the choice of baseline climate information on future projections of bamboo or panda distributions has not been considered.Earlier studies employed only a single source of baseline climate information, usually the WC dataset [41][42][43][44][45][46].Therefore, a better understanding of the relative contribution of the choice of baseline climate information to the overall uncertainty of future projections of bamboo distribution is essential for Giant Panda conservation efforts.

Data sources and preparation
Species occurrence data.Occurrence data for 21 bamboo species (S1 Table ) consumed by the Giant Panda were obtained from the Fourth National Survey on Giant Pandas [38] for more than 11,900 field survey locations.The Fourth Panda Census covered an area of 43,583 km 2 , and was conducted along 20,513 transects located within Sichuan, Shanxi, and Gansu Provinces [38,47].
Baseline climate data.The WC dataset consists of gridded climate layers for the 1960-1990 period at a resolution of 30 arc-seconds (about 1 km).The gridded layers were produced by applying a thin-plate smoothing spline algorithm to station observations of monthly precipitation and mean, minimum, and maximum temperature using latitude, longitude, and elevation as independent variables [34].The observation stations within the study region that were used in the WC interpolation are shown in Fig 1 in relation to the regional topography.Gridded fields of 19 commonly used bioclimatic variables (bio1-bio19, S2 Table ) were downloaded from http://www.worldclim.org/version1,along with long-term averages of monthly maximum, minimum, and mean temperature and monthly total precipitation.The variables were interpolated and resampled to a 1 km 2 resolution using bilinear interpolation to adjust for sampling bias and to ensure equal area cells [48].
The RS estimates of temperature and precipitation, available at a 0.05 degree (about 6 km) resolution, were compiled by Deblauwe et al. [35] for the period 2001-2013 for temperature and 1981-2013 for precipitation.The temperature estimates were obtained from monthly maximum, minimum, and mean land surface temperatures acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the National Aeronautics and Space Administration Terra satellite (MOD11C3 v. 6.0 product).Precipitation estimates were obtained from the Climate Hazards Group InfraRed Precipitation with Station version 2 (CHIRPS v. 2.0) dataset, which blends satellite imagery with in-situ (i.e., station) precipitation observations.Average monthly values were used to derive the same 19 bioclimatic variables available for WC.The bioclimatic variables were downloaded from https://vdeblauwe.wordpress.com,and re-projected using bilinear interpolation a 1 km 2 spatial resolution.
Future climate projections.WC provides future projections for 2061-2080 of average monthly temperature and precipitation and of the 19 bioclimatic variables.These projections were downscaled by the WorldClim developers from simulations for 17 GCMs in the Coupled Model Intercomparison Project Phase 5 (CMIP5) archive [49] to a 30 arc-second resolution using the popular "delta" method [50] The delta downscaling method was also used to derive similar projections for the RS dataset.We first subtracted the WC layers of average maximum, minimum, and mean temperature for the 1960-1990 observed period (what the WC developers refer to as the "current" climate) from the future projections of average monthly maximum, minimum, and mean temperature for each of the 17 GCMs available from the WC website (http://www.worldclim.org/CMIP5v1).The differences at each grid point are the delta values that the WC developers used to derive the projected future (2061-2080) values of maximum, minimum, and mean temperature from the WC baseline climate data.For precipitation, the relative differences (i.e., ratio) between the WC future mean monthly precipitation projections and the WC baseline mean monthly precipitation fields were calculated.The deltas (i.e., change factors) for the temperature and precipitation variables were then applied to the RS baseline temperature and precipitation layers that had previously been re-projected to the WC grid to obtain future projections for the RS dataset.The 19 bioclimatic variables were generated from the RS future projections of temperature and precipitation using the R package "dismo" [51].(See S1 Appendix for a more detailed explanation of the downscaling procedure.)Thus, the future projections for the RS dataset differ from those for the WC dataset only in terms of the fine resolution baseline climate layers used to downscale the GCM simulations.Pairwise spatial correlation coefficients and principal component analysis (PCA) were used to compare the WC and RS baseline and projected future climate fields.

Habitat suitability modeling
We employed the popular MaxEnt model (version 3.3.3k)[52] for the species distribution modeling.MaxEnt is a presence-only, machine learning algorithm based on maximum entropy theory for determining the niche of a species based on the environmental conditions (e.g., bioclimatic variables) of the areas where the species occurs [52].Although the choice of modeling algorithm has been shown to also introduce uncertainty into species distribution modeling outcomes (e.g., [53][54]), we opted to focus only on the MaxEnt model for this analysis because of its wide use [55].An understanding of the impact of baseline climate data on MaxEnt calibration is helpful for further interpreting the many ecological assessments that have used this model.
The climatic predictor variables used in this study were selected using PCA applied separately to the WC and RS baseline bioclimatic layers.The resulting components were first interpreted in terms of the climate variables with the highest loadings on a particular component (described below), and the variable that best captured each component's interpretation was selected.The use of a single bioclimatic variable per component avoids high correlation among the input variables.An advantage of using PCA to select variables for the MaxEnt analysis, is that, unlike simple correlation between pairs of variables, it explicitly considers the shared and unique variance of the different bioclimatic variables.During the model training process, we treated the grid cells that coincide with at least one of the field survey locations as background and selected 10000 random cells from the background.We thinned the species occurrence data using the default option in MaxEnt that removes duplicate occurrence records in a grid (see S1 Table for the number of presence locations for each species after removing the duplicate locations), and adopted the k-fold cross-validation method.Ten replications were performed for each model calibration and, for the future period, each combination of model calibration and downscaled GCM projection.The mean results of the 10 replications are presented below.Model performance for the baseline climate conditions was evaluated using area under the receiver operating characteristic curve (AUC) [56], the true skill statistic (TSS) [57], and partial AUC [58].

Uncertainty assessment
We initially focused on the logistic outputs from the MaxEnt model, which are typically interpreted as the probability of species presence [55,59].Differences in the probabilities between the future and historical periods were calculated following Guillera-Arroita et al. [59].PCA analyses were performed to explore the similarities and dissimilarities in the spatial patterns of the projections obtained from the different data sources and various GCMs.In addition, future changes in the climatically-suitable area of each bamboo species were assessed.For this analysis, the eleven conversion thresholds provided by MaxEnt (i.e."fixed cumulative value 1", "fixed cumulative value 5", "fixed cumulative value 10", "minimum training presence", "10 percentile training presence", "equal training sensitivity and specificity", "maximum training sensitivity plus specificity", "equal test sensitivity and specificity", "maximum test sensitivity plus specificity", "balance training omission predicted area and threshold", "equate entropy of thresholded and original distributions") were used to convert the projected probabilities for the baseline and future climate conditions to species presence.The inclusion of this uncertainty source in the analysis was motived by the finding of Nenze ´n and Arau ´jo [60] that threshold selection can have a large effect on estimates of future suitable conditions.The relative contribution of the three uncertainty sources (baseline climate dataset, GCM, and conversion threshold) to the change in the climatically-suitable area was analyzed for all 21 species using a three-way analysis of variance (ANOVA).

Visual inspection of the spatial patterns of the bioclimatic variables for the WC (S1 Fig) and RS (S2 Fig) datasets along with plots of the deviations between the two datasets (Fig 2) reveals
substantial differences in climatic gradients.In general, the RS dataset depicts steeper and more complex thermal gradients than does the WC dataset.This is especially evident for mean diurnal (bio2) and annual (bio7) temperature range.For WC, mean diurnal temperature range increases by approximately 1.5˚C from the eastern and southeastern edges of the study region (i.e., the western Chengdu Plain) and the lower-elevation foothills to the mountainous western area (see Fig 1 for a detailed depiction of the topography and location of major mountain ranges).In contrast, a steeper gradient is observed over the eastern foothills of the Daxiangling, Qionglai, and Minshan Mountains for RS, and the largest diurnal temperature ranges of over 2.0˚C are located west of the Minshan Mountains.For WC, the spatial pattern of annual temperature range is generally similar to that for diurnal temperature range, but larger deviations between these two variables are seen for RS.In particular, strong gradients in annual temperature range are evident west of the southern mountain ranges with temperature range decreasing with elevation, and over the northwestern portion of the study region where the largest annual temperature ranges are found west of the Minshan Mountains.Spatial gradients for isothermality (bio3), defined as the ratio of diurnal to annual temperature range, are considerably steeper for RS than WC.Similarly, the spatial distribution of temperature seasonality (bio4) is more complex for RS, whereas a broad north-to-south gradient is evident for WC.During the cold season, mean temperatures (bio6, bio9, bio11) for RS are colder than those for WC in the foothills and lower elevation mountainous regions but warmer in the higher elevation western portion of the study region, again contributing to a steeper thermal gradient for RS.The fine detail structure of annual (bio1), warm season (bio5, bio8, bio10), and cold season (bio6, bio9, bio11) mean temperature for WC closely corresponds with elevation with the coldest temperatures found at high elevations and warmer temperatures at lower elevations.For example, warmer temperatures in mountain valleys are clearly evident in the maps of WC mean temperature.Greater between-season differences are seen for RS.In the warm season, relatively steep temperature gradients are observed west of the Xiaoxiangling Mountains, in the northwestern Qionglai Mountains and along the foothills of the Minshan Mountains.Also, warmer temperatures are observed west of the Minshan Mountains compared to surrounding higher elevation locations.These gradients are either weaker or not present during the cold season.
The two datasets depict a precipitation maximum, evident in both the warm (bio13, bio16, bio18) and cold (bio14, bio17, bio19) seasons, along the western edge of the Chengdu Plain and the foothills of the Xiaoxiangling and Daxiangling Mountains, although this maximum covers a larger area for RS compared to WC. Mean precipitation during the warm, wet season is generally larger, and precipitation gradients weaker, for RS compared to WC.During the cold, dry season, mean precipitation is smaller for RS compared to WC except for a narrow band along the western edge of the Chengdu Plain.Precipitation seasonality (bio15) suggests larger intra-annual variations in precipitation for RS than WC, with the greatest seasonality found in the eastern foothills of the mountain ranges.In sum, the two datasets differ substantially in their fine detail spatial structure and their depiction of local and regional temperature and precipitation gradients.

Variance structure of the WC and RS datasets
To further examine the similarities and differences between the two datasets, Pearson's pairwise correlation and a varimax-rotated PCA were performed.The correlation analysis supports the visual interpretation of the spatial patterns of the bioclimatic variables presented above.The high (>0.8)correlations between WC and RS for most of the bioclimatic variables (Fig 3, main diagonal) indicate that both datasets have similar broad-scale spatial patterns.The weaker correlations for annual temperature range (bio7) and precipitation seasonality (bio15) of 0.6 and 0.7, respectively, indicate greater differences between the two datasets in the spatial patterns of these variables.In general, the within-dataset correlations between the different bioclimatic variables (upper-right portion of Fig 3 for WC, lower-left for RC) are larger (indicated by the larger and darker "circles") for WC compared to RC, although the overall pattern of the correlations is similar for the two datasets (i.e., higher correlations among temperature variables compared to among precipitation variables).In terms of correlation with elevation, temperature-related bioclimatic variables generally have a higher correlation with elevation in WC compared to RS.For both datasets, the correlation between elevation and the precipitation-related bioclimatic variables is weaker than that for temperature.
PCA was also performed since, unlike simple correlation coefficients, it considers the shared variance among variables (Table 1).A varimax rotation was used to facilitate the interpretation of the major dimensions of the two datasets.A scree plot (S3 Fig) suggests that a larger number of components are needed to summarize the variance in RS compared to WC (i.e., five rather than four components with eigenvalues >1).For WC, the first two rotated components explain 78% of the variance in the dataset with the first component broadly representing mean temperature (bio1, bio2, bio5, bio6, bio8, bio9, bio10, bio11) and the second component representing mean precipitation (bio12, bio13, bio14, bio16, bio17, bio18, bio19).In contrast, three components are needed to explain a similar portion of the variance for RS, and the temperature variables, which loaded highly on a single component for WC, are distributed over the second and third components for RS, which loosely can be interpreted as warm season mean temperature and cold season mean temperature.This is in agreement with the larger differences observed above in the spatial patterns of mean temperature between the warm and cold seasons for RS than for WC.On the other hand, the first component for RS represents precipitation, and the variables loading highly on this component are identical to those on the WC precipitation component.For both datasets, precipitation seasonality represents a unique dimension (the third component for WC and fifth component for RS), whereas temperature seasonality is a unique dimension only for WC and temperature range (diurnal and annual) represents a unique dimension only for RS.

Variable selection and interpretation
MaxEnt model calibrations were performed under three sets of climate predictors chosen based on the correlation and PCA analyses.For the first calibration (WC4), one bioclimatic -1 0 0.8

WorldClim
Remote Sensing 0.9 0.8 0.9 0.9 0.9 0.9 0.6 0.9 0.9 0.9 0.9 0.9 0.9 0.8 0.7 0.9 0.8 0.8 0.   variable was selected to represent each of the four rotated principal components of the WC dataset.The variables chosen were: annual mean temperature (bio1); temperature seasonality (bio4); annual precipitation (bio12); and precipitation seasonality (bio15).For the second calibration (RS4), the same four predictors but from the RS dataset were used to compare with the WC4 calibration.A third calibration (RS5) employed five predictors (adding bio7, temperature annual range) from the RS dataset.This resulted in 63 model combinations (21 bamboo species x 3 calibrations) for the baseline climate conditions times 10 replications for a total of 630 simulations.The AUC values ranged from 0.83 to 0.99 for all three calibrations, indicating a high predictive power of the models under the baseline climate conditions, which was confirmed by the values of TSS and partial AUC (S3 Table and S4 Table ).
A comparison of the outputs from the different calibrations highlights that the choice of baseline climate data can influence the interpretation of the relative importance of the environmental variables to the distribution of the species under consideration.Regardless of calibration, only two variables contributed substantially (> 20%) to the MaxEnt model for most of the bamboo species, although some exceptions are evident (Fig 4).For the WC4 and RS4 calibrations, the order of the leading two variables is in agreement for only three of the bamboo species (Fargesia ferax, Fargesia scabrida, Yushania maculate).The estimated contribution is largest for temperature seasonality (bio4), followed by annual precipitation (bio12), for these species.For another eight bamboo species (Bashania faberi, Bashania fargessi, Bashania spanostachya, Fargesia dracocephala, Fargesia qinlingensis, Yushania ailuropodina, Yushania brevipaniculata, Yushania glauca), the WC4 and RS4 calibrations agree on the primary variable in the MaxEnt model, but differ on the variable with the second largest contribution.The leading variable is temperature seasonality for six of the species, but for the other two annual precipitation is the primary contributor.For the remaining 10 species, the WC4 and RC4 calibrations disagree on which of the bioclimatic variables contribute the most to the MaxEnt model, and hence provide substantially different interpretations of the relative influence of the bioclimatic variables on the distribution of these bamboo species.The contribution of the leading two variables is almost identical for RS4 and RS5 with two exceptions.For Bashania spanostachya and Yushania maculata, the contribution of the leading variable for RS5 is smaller than that for RS4, with the other variables contributing more to model calibration.

Future probability of species presence
The species distribution models were employed to project the future (2061-2080) probabilities of species presence using the downscaled bioclimatic variables from 17 GCMs.A total of 10,710 simulations was performed (21 bamboo species x 3 calibrations x 17 GCM projections x 10 replications), and the outputs of the 10 replications were averaged for each bamboo-calibration-GCM combination.Changes in probability were found for each bamboo species by subtracting the probabilities between the future and baseline periods.To help summarize the differences in the spatial patterns, an un-rotated PCA was performed separately for each   .In an unrotated PCA analysis, the first component accounts for the highest variance among all components and represents the direction of the highest agreement among the projections for each species, and the second component, which is orthogonal to the first component, accounts for the largest amount of the remaining variance.On the biplot, each projection is represented by a vector, and the similarity between the different projections is shown by the angle between the vectors (small angle indicates high similarity).
The striking feature of Fig 5 is the large angles between the projections obtained from the WC4 calibration (red arrows) and those obtained from the RS4 (green arrows) and RS5 (blue arrows) calibrations.In contrast, the arrows of the same color are closer together, indicating that the angles between the projections from the different GCMs but same calibration are relatively small.Thus, the choice of baseline climate data source has a larger influence on the spatial patterns of the projected change in species probability than the choice of GCM.Also, for almost all bamboo species, the projections obtained from the RS4 and RS5 calibrations display a high similarity, which suggests that the inclusion of one additional predictor in the RS5 calibration does not have as large an influence on the projections as the initial source of the baseline climate information.Some deviations from these generalizations are evident, however.The angles between the projections obtained from the WC4 calibration and those from the RS4 and RS5 calibrations are smaller, with even some "intermingling" of the vectors, for Bashania spanostachya, Fargesia ferax, and Fargesia robusta.Also, a wider spread among the RS4 and RS5 projections is seen for six of the bamboo species (Chimonobambusa szechuanensis, Fargesia dracocephala, Fargesia qinlingensis, Yushania ailuropodina, Yushania glauca, and Yushania maculata), although the angle between the WC4-and RS-based projections remains large.For the remainder of the bamboo species, though, two distinct vector clusters are evident, one representing the projections obtained from the WC4 calibration and the other cluster composed of the projections obtained from the RS4 and RS5 calibrations.

Projected changes in climatically-suitable area
The projected change in the climatically suitable area was estimated by applying eleven conversion thresholds to the 51 (3 calibrations x 17 GCMs) species presence probabilities for each of the bamboo species and to the probabilities obtained for the baseline climate conditions from the three model calibrations (Fig 7).Substantial differences in future climatically-suitable area are observed between the three calibrations.For three of the bamboo species (Fargesia obliqua, Qiongzhuea opienensis, and Yushania ailuropodine), the direction of the median projected change varies by calibration.For another seven species (Bashania spanostachya, Chimonobambusa szchuanensis, Fargesia dracocephala, Fargesia qinlingensis, Fargesia robusta, Yushania brevipaniculata, and Yushania lineolata), the median value for at least one of the calibrations suggests a marked future change, either an increase or decrease in climatically-suitable area, whereas the median values for one or both of the other calibrations point to a considerably smaller change.On the other hand, the median projected change for all three calibrations agree that climatically suitable area will increase for six species (Bashania fargesii, Fargesia ferax, Phyllostachys nidularia, Qiongzhuea tumidinoda, Yushania glauca, and Yushania maculata) and decrease for another five species (Bashania faberi, Fargesia denudata, Fargesia nitida, Fargesia rufa, and Fargesia scabrida), although for some of these species (Qiongzhuea tumidinoda and Yushania glauca) the magnitude of the projected changes varies substantially.In most cases, the magnitude of the median projected changes is more similar for the RS4 and RS5 calibrations than for the WC4 calibration.
A three-way ANOVA was performed on the WC4 and RS4 projections (374 total projections for each species; 2 calibrations x 17 GCMs x 11 thresholds) in R software to further explore the influence of the different uncertainty sources on the projected changes in climatically-suitable area.The contributions of the uncertainty sources vary substantially by bamboo species (Fig 8 ).For example, the variance among the projections for Fargesia scabrida is primarily explained by the different GCM simulations, whereas for Fargesia qinlingensis and Yushania ailuropodina it is mostly explained by the baseline climate data used in the calibration.The interaction between the baseline climate data and the GCMs is most notable for Fargesia dracocephala.With the exception of Fargesia nitida, contributions of the conversion threshold and its interaction terms are relatively small for the different bamboo species.When summarized across all the bamboo species, the choice of GCM is the single largest contributor to the total variance in the projected future climatically-suitable area, although the sum of the variance explained by the baseline climate data through its effect on species model calibration and the interaction of the baseline climate data and GCM simulation via downscaling is also large.The mean contribution of the GCM to the explained variance is about 38%, and the sum of the mean contribution of the data source and the interaction between data source and GCM accounts for approximately 40% of the total variance.All other factors have a mean contribution below 8%.

Discussion
We evaluated the potential contribution of the choice of baseline climate dataset to the uncertainty of projected future changes in the probability of species presence and climatically-   suitable area, using 21 understory bamboo species occurring within the mountainous geographic range of the Giant Panda in southwest China as examples.We found that variations in local and regional gradients of temperature and precipitation between the WC and RS datasets led to differing interpretations of the relationships among 19 bioclimatic variables commonly used in species distribution modeling.For all but three of the 21 bamboo species, the MaxEnt model ranked the bioclimatic variables differently when calibrated using the two climate datasets, and thus interpretation of the climatic determinants of the distribution of the bamboo species varied by dataset.Moreover, greater differences in the spatial patterns of the projected changes in the probability of species presence were observed between projections obtained from the species distribution models calibrated using the different baseline datasets than between the different downscaled GCMs simulations for the same calibration.Also, a multi-factor ANOVA revealed that a substantial portion of the variance among the future projections of climaticallysuitable area can be explained by the main effect and interactive terms involving the source of baseline climate data.The main effect term represents the influence of the baseline climate dataset on the calibration of the MaxEnt model.We interpret the interaction term as primarily  reflecting the influence of the baseline climate dataset in the downscaling of the GCM simulations.However, the interaction term may also capture the interaction between model calibrations and differences among GCMs in the projected changes of the bioclimatic variables.One simplification of our evaluation is the use of the same change factors (i.e., deltas) for both the WC and RS datasets.Such simplification is perhaps more problematic for the temperature-derived bioclimatic variables, since they constitute measurements of air temperature in WC while estimates of land temperature in RS.We used this approach to focus the analysis on the impact of the different baseline climate datasets.The use of the same delta values for both datasets assumes that air and land temperatures will change by similar amounts in the future, and that the satellite-derived precipitation measurements constitute suitable representations of observed precipitation.Although a similar assumption was made by [61] when developing fineresolution climate projections for Africa, the change factors for future assessments ideally should be calculated from GCM or RCM simulations of variables more similar to those acquired remotely by spaceborne sensors, such as land surface temperature.Also, our study included only future climate projections forced by RCP 8.5 in the uncertainty analysis.However, we anticipate that the relative contribution of the baseline climate information to the overall uncertainty will be larger for the less extreme RCPs, especially as differences in the projected future changes in temperature and precipitation among the 17 GCMs will be smaller reducing the magnitude of this uncertainty source.Thus, the use of RCP 8.5 represents a harsher test of the contribution of the choice of baseline climate data to the overall uncertainty.
The comparison of the spatial patterns for the bioclimatic variables obtained from the WC and RS datasets presented here emphasizes Daly's [25] concern regarding the tendency to equate resolution with realism, and, furthermore, that the fine resolution of gridded datasets may give an appearance of realism that is not consistent with the spatial resolution of the initial observations used to generate them.While the WC dataset has a resolution of ~1 km compared to the coarser ~6 km resolution of the RS dataset, the native resolution of the measurements from which the datasets are derived is much coarser for WC (tens to hundreds of km) compared to RS (30 m to 6 km, depending on the sensor system).Hence, the scale at which the gridded climate layers are provided is closer to the native resolution of RS than of WC, while the finer resolution of the WC dataset is primarily introduced via the resolution of the elevation layer used in the spatial interpolation.
Our analyses also support the findings of Baker et al. [19] that the contribution of the baseline climate dataset to the uncertainty surrounding estimates of future climatic suitability can rival that introduced by future climate projections.Of particular note is that their assessment of climatically-suitable area for avian species in sub-Saharan Africa considered climate datasets interpolated to only 0.44˚(~50 km), compared to the ~1 km focus of our study, and thus emphasized regional to sub-continental climatic gradients rather than the local to regional gradients of the mountainous panda habitat.The similar conclusions of the two studies imply that the choice of baseline dataset is a concern for assessments conducted at a range of spatial scales and in a variety of environments.This interpretation does not appear to depend on the source of future climate projections.In our analysis, we employed downscaled future projections obtained from one realization of 17 different GCMs, focusing on the structural and parameterization differences among GCMs.In contrast, Baker et al. [19] used downscaled projections from five realizations of a single GCM where the initial conditions were modified to evaluate natural variability.In fact, their findings suggest that uncertainty introduced via natural climate variability [62] should be explicitly included in future climate change assessments in addition to that introduced by the choice of GCM.Both studies also found that the relative contributions of the different uncertainty sources varied by species, likely a function of differences in the complexity of the climatic environment in which individual species reside.
An interesting finding of our analysis is that the choice of conversion thresholds to convert probabilities to species presence contributed to only a small portion of variance in the future projections.Nenze ´n and Arau ´jo [60] found that the conversion threshold can induce a 1.7 to 9.9-fold difference in the proportions of species projected to become threatened by climate change.Our results instead suggest that the choice of baseline climate dataset and GCM introduces more uncertainty to the climate change assessment than the choice of conversion threshold, although the differences between the two studies need to be interpreted cautiously as Nenze ´n and Arau ´jo [60] considered an overlapping, but not duplicative, set of potential uncertainty sources.In particular, they used a wider range of threshold-setting methods and, additionally, several approaches to bioclimatic modeling were used rather than the single modeling framework (i.e., MaxEnt) used here.Also, a single baseline climate dataset was employed for all their modeling efforts, compared to the two contrasting datasets included in our analysis.Thus, the relative magnitude of an individual uncertainty source may vary depending on the other uncertainty factors included in the analysis.

Conclusions
The uncertainty associated with climate change assessment findings must be carefully considered.Neglecting this uncertainty can lead to misinformed research agendas, policies, and decisions.Our findings highlight the importance of routinely considering the sensitivity of future projections of species distribution to the choice of baseline climate information, especially in mountainous environments with complex climatic gradients.This uncertainty source is a particular concern as, for many climate change assessments, the baseline climate dataset is used both to calibrate species distribution models and to downscale coarse-scale projections of the future climate to local/regional scales.
Future studies that consider a more complete set of uncertainty factors than those used in the past will advance our understanding of the effects of climate change on species distributions, and the assessment findings will better guide mitigation and adaptation policies and conservation practices for reducing the threats posed by climate change on biodiversity.climate information and four bioclimatic variables; and "RS5" refers to the model calibration using the remotely-sensed baseline climate information and five bioclimatic variables.(DOCX) S4 Table .The partial area under the receiver operating characteristic curve (AUC) test results for 10 replication runs for 21 bamboo species and three calibrations.The calculations were performed using an R package available at https://github.com/vijaybarve/ENMGadgets.The calculations employed a random sample of 50% of the testing points for bootstrapping with 100 replications using thresholds for omission errors greater than 0.1, 0.5, and 0.9, respectively."WC4" indicates the model calibration using the WorldClim baseline climate information and four bioclimatic variables; "RS4" refers to the model calibration using the remotely-sensed baseline climate information and four bioclimatic variables; and "RS5"refers to the model calibration using the remotely-sensed baseline climate information and five bioclimatic variables.The bioclimatic variables were generated based on temperature (units of ˚C Ã 10) and precipitation (units of mm).The bioclimatic variables are long-term averages of annual mean temperature (bio1); mean diurnal range (bio2); isothermality (bio3); temperature seasonality (bio4); maximum temperature of the warmest month (bio5); minimum temperature of the coldest month (bio6); annual temperature range (bio7); mean temperature of the wettest (bio8), driest (bio9), warmest (bio10), and coldest (bio11) quarter; annual precipitation (bio12); precipitation of the wettest (bio13) and driest (bio14) month; precipitation seasonality (bio15); and precipitation of the wettest (bio16), driest (bio17), warmest (bio18) and coldest (bio19) quarter.The bioclimatic variables are generated based on temperature (units of ˚C Ã 10) and precipitation (units of mm).The bioclimatic variables are long-term averages of annual mean temperature (bio1); mean diurnal range (bio2); isothermality (bio3); temperature seasonality (bio4); maximum temperature of the warmest month (bio5); minimum temperature of the coldest month (bio6); annual temperature range (bio7); mean temperature of the wettest (bio8), driest (bio9), warmest (bio10), and coldest (bio11) quarter; annual precipitation (bio12); precipitation of the wettest (bio13) and driest (bio14) month; precipitation seasonality (bio15); and precipitation of the wettest (bio16), driest (bio17), warmest (bio18) and coldest (bio19) quarter.

Fig 3 .
Fig 3. Pearson pairwise correlation matrix between the bioclimatic variables for the WorldClim (WC) dataset (upper-right matrix) and for the remotely-sensed (RS) dataset (lower-left matrix) shown in colored circles, and the correlation between the RS and WC datasets for the same bioclimatic variables (diagonal values from upper-left to lower-right).The numbers from 1 to 19 stand for the 19 bioclimatic variables, while "E" stands for elevation.Blue colors indicate positive correlations and red colors indicate negative correlations as indicated in the figure legend.Both the size and the color of the circles represent the magnitude of the correlation.The bioclimatic variables are long-term averages of annual mean temperature (bio1); mean diurnal range (bio2); isothermality (bio3); temperature seasonality (bio4); maximum temperature of the warmest month (bio5); minimum temperature of the coldest month (bio6); annual temperature range (bio7); mean temperature of the wettest (bio8), driest (bio9), warmest (bio10), and coldest (bio11) quarter; annual precipitation (bio12); precipitation of the wettest (bio13) and driest (bio14) month; precipitation seasonality (bio15); and precipitation of the wettest (bio16), driest (bio17), warmest (bio18) and coldest (bio19) quarter.

Fig 4 .
Fig 4. Contribution of each of the top three bioclimatic variables for the three model calibrations."WC4" indicates the model calibration using the WorldClim baseline climate information and four bioclimatic variables; "RS4" refers to the model calibration using the remotely-sensed baseline climate information and four bioclimatic variables; and "RS5" refers to the model calibration using the remotely-sensed baseline climate information and five bioclimatic variables.The number on the top of each bar refers to the bioclimatic variables: bio1, annual mean temperature; bio4, temperature seasonality; bio7, temperature annual range; bio12, annual precipitation; and bio15, precipitation seasonality.The height of a bar is the percent contribution of a particular variable.https://doi.org/10.1371/journal.pone.0189496.g004

Fig 6 (
and S4 Fig) illustrates the substantial spatial variations between the different calibrations in the projected changes in species presence probability for Fargesia denudata.Such differences were observed for many of the other bamboo species evaluated.Almost all of the projections obtained from the WC4 calibration, regardless of the GCM simulation to which the species distribution model was applied, indicate a decreased probability in the current range (the northeastern Minshan Mountains) of Fargesia denudata but an increase in the northwest portion of the study region.In contrast, almost all of the projections obtained from the RS4 and RS5 calibrations suggest a decreased probability in the current range of Fargesia denudata and for locations northwest of the current range.

Fig 5 .Fig 6 .
Fig 5. Biplots for the 21 bamboo species of the change in the likelihood of occurrence obtained from the three model calibrations expressed as the difference between the future (2061±2080) and baseline climate conditions under RCP 8.5 for 17 downscaled GCMs."WC4" indicates the model calibration using the WorldClim baseline climate information and four bioclimatic variables; "RS4" refers to the model calibration using the remotely-sensed baseline climate information and four bioclimatic variables; and "RS5" refers to the model calibration using the remotely-sensed baseline climate information and five bioclimatic variables.https://doi.org/10.1371/journal.pone.0189496.g005

Fig 7 .
Fig 7. Projected change in the climatically-suitable area for the 21 bamboo species by the three model calibrations developed from the WorldClim (WC) and remotely-sensed (RS) datasets using either four for five bioclimatic variables as predictors (WC4, RS4, and RS5).The projected change is expressed as the ratio of the difference in climatically-suitable area between the future (2061-2080) and baseline climate conditions to the climaticallysuitable area for the baseline conditions (the values can be multiplied by 100 to obtain a percentage change).Each box and whisker plot includes projections obtained from 17 global climate models (GCMs) and 11 conversion thresholds.https://doi.org/10.1371/journal.pone.0189496.g007

Fig 8 .
Fig 8. Proportion of sum-squared error to total squared error for a three-way analysis of variance (ANOVA) of the projected change in suitable area for each of the 21 bamboo species (left), and summarized over all 21 bamboo species (right).The analysis includes two model calibrations developed from the WorldClim (WC) and remotely-sensed (RS) datasets and using four bioclimatic variables as predictors (referred to as WC4 and RS4), 17 future climate projections downscaled from global climate models (GCMs), and 11 conversion thresholds.https://doi.org/10.1371/journal.pone.0189496.g008 (DOCX) S5 Table.Centroids (WGS 1984 UTM Zone 48N projection) of the 1 km x 1 km grid cells with at least one presence point for the 21 bamboo species included in the MaxEnt modeling.(XLSX) S1 Fig. Spatial variations of the 19 bioclimatic variables for the baseline WorldClim (WC) dataset and the elevation of the study region.
(PDF) S2 Fig. Spatial variations of the 19 baseline bioclimatic variables for the baseline remotelysensed (RS) dataset and elevation of the study region.

Table 1 . The varimax-rotated principal component loadings for the WorldClim (WC) and remotely-sensed (RS) datasets
. "RC" refers to rotated principal components with eigenvalues > 1.Four principal components were extracted for the WC dataset and five for the RS dataset.The variables shown in bold were used in the calibration of the species distribution models. https://doi.org/10.1371/journal.pone.0189496.t001