Change of niche in guanaco (Lama guanicoe): the effects of climate change on habitat suitability and lineage conservatism in Chile

Species occurrence data

We built a guanaco occurrence dataset of 3,321 records by complementing previous work by the authors with 359 additional records (González et al., 2013). New records were collected following the same procedures outlined in González et al. (2013), that is, from direct and indirect evidence of guanaco presence collected between the years 2000 and 2016 across several field campaigns. Indirect evidence of guanaco occurrence was assigned to a lineage by genetic and morphological evaluation of biological samples such as feces and dead tissues. Most of new records were collected in the northern section of the country in the Arica, Parinacota (i.e., 17°S latitude) and Coquimbo region (30°S). Each record was assigned to a 1 × 1 km cell defined by the resolution of the environmental datasets employed (see below). This resulted in a total of 298 records for L.g. cacsilensis, 837 for the mixed population, and 2,186 for L.g. guanicoe.

Climate predictors

We limited the selection of environmental predictors to climatic variables (Thuiller et al., 2006; Hu et al., 2015). Similarly to what has been described in the literature (Thuiller et al., 2004), our previous work dismissed the importance of nonclimate predictors for guanaco distribution models in favor of exclusive climatic variables (González et al., 2013). We used all 19 bioclimatic variables from WorldClim (version 1.4) summarizing temperature and precipitation information worldwide (Hijmans et al., 2005). To reduce collinearity, model overfitting, and the number of explanatory variables, we used a paired correlation analysis to inspect pairs of variables and removed variables with a large correlation coefficient (>0.8) (Beaumont, Hughes & Poulsen, 2005).

While the analysis was limited to the administrative bounds of Chile, all WorldClim variables were projected to UTM 19 South, with a one squared-kilometer of spatial resolution, spanning from latitudes 15° to 55°S and longitudes 60° to 80°W and a total area of 5,921,578 km² covering most of the Southern Cone.

Future Climate Projections

The projection of future geographic distribution of niches was performed using the outputs of the Coupled Model Intercomparison Project 5 of the IPCC’s methodology for the Fifth Assessment Report (AR5) (Taylor, Stouffer & Meehl, 2012). The two extreme GHG concentration scenarios, also known as representative concentration pathway (RCP), were used to project future climate niches. RCP2.6, the most optimistic scenario, considers a lower GHG concentration and projects average increases of temperature between 0.3° and 1.6 °C with 0.26–0.55 m increases of sea levels. RCP8.5, the most pessimistic scenario, considers higher GHG concentrations with a 2.6°–4.8 °C projected increase in mean global temperature and a 0.45–0.82 m rise of sea levels (IPCC, 2013). We chose both extreme scenarios to evaluate the minimum and maximum potential impact of climate change in the guanaco’s distribution.

Given the large uncertainties of future climate predictions, the computing power availability and the exploratory nature of such models, we selected five general circulation models (GCM) among the 19 models used to generate the AR5. GCM’s are physical climate models that simulate the interactive biophysical processes between the atmosphere, the ocean and the land (Moss et al., 2010). Selected climate models were: (1) CCSM4 model of the National Center of Atmospheric Research (Gent et al., 2011); (2) GFDL-CM3 model of the Geophysical Fluid Dynamics Laboratory (Donner et al., 2011); (3) GISS-E2-R model of the NASA Goddard Institute for Space Studies (Nazarenko et al., 2015); (4) HadGEM2-AO atmosphere model and (5) HadGEM2-ES earth system model, both of the Met Office Hadley Centre (Collins et al., 2011). Each scenario was evaluated for the short (2050) and medium term (2070).

Niche modeling

Entropy maximization procedures in MaxEnt 3.3.3 k (Phillips, Anderson & Schapire, 2006) were used to model current and future geographic distributions of L. guanicoe and its lineages. MaxEnt uses a machine learning algorithm to generate predictions on the potential distribution of species based on their presence, pseudo-absences and a set of environmental variables. The software analyzes the multivariate distribution of environmental conditions of species occurrences to generate a spatially explicit probability map of lineage occurrence (Franklin, 2009). Such modeling approach has shown to have a good statistical performance compared to other types of modeling techniques (Elith et al., 2006) and is currently one of the most commonly used methods to understand habitat suitability, niche structure, geographical species distribution (Merow, Smith & Silander, 2013) as well as to project environmental niches to future scenarios (Hijmans & Graham, 2006).

We performed 100 cross-validated replicates for each current and projected distribution model with logistic output, that unlike other outputs (i.e., raw and cumulative) assumes that a known observation probability can be assigned to each pixel and has thus been considered as a true approximation of presence (Merow, Smith & Silander, 2013). The “fade by clamping” option was used to avoid predictions beyond the observed geographical range during the training of the future distributions models (Phillips, Anderson & Schapire, 2006). All other parameters were kept at their default configuration (Phillips, Dudík & Schapire, 2004) as they have previously shown good performance in ungulate modeling (Hu & Jiang, 2011; González et al., 2013; Hu et al., 2015; Quevedo et al., 2016) and other taxonomic groups (Phillips & Dudík, 2008; Fourcade et al., 2014).

We used an ensemble forecasting framework to minimize the inherent variability introduced by the various forecast models employed, as proposed by Araújo & New (2007). Therefore, we generated a model from the average of each bioclimatic variable produced by the five GCMs (i.e., Bio1_CC + Bio1_GF Bio1_GS + Bio1_HD + Bio1_HE), and then evaluated an average value for each variable from 100 replicates for both extreme RCP emission scenarios for the years 2050 and 2070. Hence, 400 projected guanaco distribution models were generated (i.e., 2 RCPs × 2 time frames × 100 replicates). The final results are four projected climate models for L. guanicoe, one for each RCP2.6 and RCP8.5 scenarios evaluated for years 2050 and 2070.

Model evaluation, prediction, and spatial projection

Generated niche models were evaluated using a threshold-independent analysis of the area under the curve (AUC) provided by the receiver operator curve (Phillips, Anderson & Schapire, 2006; Acevedo et al., 2010; Anderson & Raza, 2010). These sensitivity tests model accuracy by calculating the proportion of true positives versus false positives. The resulting values range from 0 to 1, where model predictions are considered fair when obtained AUC values are above 0.7 (Swets, 1988; Merow, Smith & Silander, 2013). A 3:1 ratio was used to divide training and testing datasets (Phillips, Anderson & Schapire, 2006). AUC Jackknife analysis allowed to identify the contribution of each variable to final current and future models, and to allow the detection of those variables that significantly improve predictions for the occurrences of each lineage (Phillips, Anderson & Schapire, 2006).

We reclassified predicted habitat using a 0.25 threshold interval to label three habitat suitability classes: low suitability habitat when occurrence probability ranged between 25% and 50%; suitable habitat if occurrence probability was in the 50–75% interval; and high suitability habitat if occurrence probability was over 75%, values below 25% were considered as inappropriate habitat (Hu & Jiang, 2012; González et al., 2013; Shrestha & Bawa, 2014).

Changes in distribution surface and incidence in protected areas

The areal extent for each suitability class predicted by each model, current and projected, were compared to determine habitat loss (or gains) under the various climate change scenarios evaluated. We used the software BioSARN v.1 to calculate the amount of area gained or lost and to estimate differences between models (Heap, 2016). These results were classified into three categories: (a) Areal loss, when future prediction show a decrease of the areal extension compared to current niche models; (b) Areal gain, produced when future prediction add area to current niche models; (c) Unchanged areas, when climate change predictions show no impact on current guanaco’s distribution.

In addition, the fraction of future distribution covered by the system of PA in Chile was estimated. All categories offering some level of protection were considered: national parks, national reserves, biosphere parks, national monuments, national patrimony, and private PA as of 2011. RAMSAR sites (as of 2012) were also included as they constitute the most important feeding grounds for guanaco in the hyper-arid north of Chile (Squeo et al., 2006).

Evaluation of PNC or PND

Phylogenetic niche conservatism and phylogenetic niche divergence among lineages and the mixed population were evaluated through their current and projected niches for the most extreme scenario (i.e., RCP8.5) in 2070 using ENMTools v.1.4.3 (Warren, Glor & Turelli, 2010). Niche overlap between lineages was calculated with the statistical indices “I” (derivative of Hellinger’s distance) and “D” (Shöener’s D) which may take values going from 0 (i.e., no overlap) to 1 (i.e., full overlap between lineages’ niches). We used the Identity test to evaluate the hypothesis of niche equivalence (Graham et al., 2004) and the Background similarity test to evaluate the hypothesis of niche similarity (Peterson, Soberón & Sánchez-Cordero, 1999). The Identity test quantitatively assesses whether the niche space for two compared lineages are equivalent by comparing the actual niche to a null niche model generated from a randomized pool of locations for each lineage. This allows to effectively evaluate whether niche spaces are equal, under the premise that, if they are, they should be able to predict each other (Warren, Glor & Turelli, 2010). Because the Identity test strongly depends on accurate representations of species habitat suitability, it is known to be sensitive to the particular sampling scheme employed, and therefore less suitable to compare allopatric niches (Warren, Glor & Turelli, 2010). The Background similarity test compares the niche difference between allopatric lineages by contrasting the niche of a “focal” lineage to the niche built from the background locations of a second lineage. If there is similarity between these, the null model should predict the niche of the second lineage. We repeated each test 100 times to produce a simulated distribution of I and D values and to evaluate significance using a threshold of 0.1 (two-tailed for background similarity test, and one-tail for identity test) (Warren, Glor & Turelli, 2010; Guisan, Thuiller & Zimmermann, 2017). We considered outcomes as indicative of PNC between lineages when observed I and D values fell within the simulated distribution. On the other hand, when the observed values fell outside of the simulated distribution, they were assumed to be indicative of PND between lineages.

Selection of climate variables and current distribution model

After removing correlated variables, the final subset of independent bioclimatic variables used in this analysis was composed of: annual mean temperature (Bio1), temperature seasonality (Bio4), annual temperature range (Bio7), annual precipitation (Bio12), and precipitation seasonality (Bio15). See correlation analysis in supplemental Fig. S1.

The major contribution to the current distribution of L. guanicoe was given by the annual range of temperature (28.2%, AUC = 0.84), whereas L.g. cacsilensis was dominated by precipitation seasonality (66.7%, AUC = 0.95), the mixed population by annual precipitation (36.7%, AUC = 0.93), and the southernmost L.g. guanicoe was mostly driven by annual mean temperature (35.8%, AUC = 0.91). Detailed Jackknife’s analysis can be seen in supplemental information, Fig. S2.

All current distribution models generated for L. guanicoe and its lineages, presented a good performance with mean AUC values over 0.89. The resulting geographic range for guanaco spanned for about a third of the Chilean continental surface. The geographic areas covered by L.g. cacsilensis and L.g. guanicoe were of 47,148 and 100,539 km², respectively. The mixed population showed a geographic extent of 84,976 km². Interestingly, our models had a 20% difference when comparing areas from the sum of lineages modeled independently and the total area modeled with all the lineages pooled as if they were a single lineage (232,664 vs. 284,499 km², respectively). Full maps and predictions are available in supplemental information Fig. S3.

Projected distribution models

As for current distribution models, the sensitivity analysis yielded a large mean AUC >0.9. After suitability categorization (Fig. 1), our results show that while the geographical distribution pattern of guanaco is conserved, quantitative assessment of the distribution surface reveals a downward trend in both scenarios of climate change for 2050 and 2070 (Table 1).

Marginal decreases of habitat suitability is observed under the RCP2.6 scenario for the years 2050 and 2070 (259,577 vs. 254,979 km², respectively). The small areal reduction under both models is only of 6% and 7.6% of the current area. However, this decrease is more pronounced under the RCP8.5 scenarios in which a 13% and 20.7% reduction is quantified for the years 2050 and 2070, respectively (i.e., 240,505 and 218,841 km²).

Surface change between current and projected distribution models

Projected distribution models under the more optimistic scenario showed an increase of high quality habitat and a net loss of medium and low quality habitat, while the projections under the worst scenario indicate a generalized decrease in habitat suitability. Although no important loss of potential distribution areas are apparent, a large decrease of areas with high quality habitat may occur under the worst climate change scenario evaluated (Table 1).

Surface losses and gains under future climate projections are described in Table 2. Both scenarios showed that a large fraction of the guanaco range will remain unchanged. The more optimistic projection (i.e., RCP2.6) indicated an average loss of 67,042 km² between 2050 and 2070, and a niche displacement (i.e., gain) of 48,225 km² on average between such time periods. A reversed trend was observed under the more pessimistic scenario (i.e., RCP8.5) with the larger change predicted for 2070. Such prediction forecasts a reduction in guanaco’s niche by 37%, with a surface loss of 103,367 km² and a geographic distribution of 172,786 km². Likewise, the smallest niche displacement was observed for this period, with 46,089 km² of areal gain (Fig. 2).

Evaluation of PNC and PND

Niche overlap and equivalence tests showed a large overlap between the contact population and L.g. cacsilensis for niche models under current and worst scenarios (Table 3). While limited niche overlap was reported, the overlap between the northern lineage and the mixed population increased under future climate change scenarios. On the other hand, the overlap between current and projected niches for L.g. guanicoe with the other groups was smaller, particularly with the northern lineage (L.g. cacsilensis). When using the results of this latter analysis as “empirical values” to perform identity and background similarity tests (Tables 3 and 4), we were able to show that statistical differences existed when comparing current niches and projected niches. This indicates that the climatic requirements between the lineages and mixed population are not equivalent.

Background similarity test showed that L.g. cacsilensis share climatic similarities with the mixed population (Table 4). However, such similarity is not reciprocal when comparing the climatic requirements of the mixed population to those of the northern lineage—a possibility mentioned in Warren, Glor & Turelli (2010). A similar situation occurred when comparing projected niches under the worst climate change scenario (RCP8.5), where results suggested that L.g. guanicoe´s niche will closely resemble the future climatic niche of the mixed population, in spite of not sharing any current similarity.

The remaining comparisons between climatic niches, current and future, showed that the climatic requirements for each of the lineages analyzed are significantly different (Table 4).

Projected distribution models and conservation in Chile

The current network of PA in Chile covers a vast area of approximately 256,550 km², according to 2016 data, and shows a limited overlap with modeled guanaco distribution. In fact, our analysis shows that a 9.8% (i.e., 19,402 km²) of the species current distribution overlaps with a PA. When looking at projected distributions for 2070, our results showed that such overlap will decrease to 6.2% (15,772 km²) under the best scenario (RCP2.6). Similarly, under the worst scenario (RCP8.5), the overlap will be of 5.7% (12,434 km²) (Fig. 3).