Introduction
Conservation planning and forecasting rely on detailed knowledge of the ecological and geographical distribution of species. Species distribution models provide detailed predictions about the potential distribution of species by relating presence records to relevant environmental factors (Phillips et al., Reference Phillips, Anderson and Schapire2006; Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011; Jetz et al., Reference Jetz, McPherson and Guralnick2012).
If implemented accurately species distribution models are a powerful and repeatable means of mapping the potential distribution of species (Wintle et al., Reference Wintle, Elith and Potts2005). Models based on bioclimatic variables at macro scales have proven successful in predicting known distributions, and refined algorithms perform well with presence-only data and a limited number of localities (Elith & Leathwick, Reference Elith and Leathwick2009). A major goal of species distribution models is to predict which areas within a region meet the characteristics of a species’ ecological niche, which is part of the species’ potential distribution (Anderson & Martínez-Meyer, Reference Anderson and Martínez-Meyer2004). In this context, MaxEnt, a machine-learning algorithm based on the maximum entropy theory, has been shown to outperform alternative presence-only models when sample sizes are small (Austin, Reference Austin2007; Baldwin, Reference Baldwin2009) and organisms have a restricted range (Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011). MaxEnt is particularly flexible in fitting complex responses: it estimates the target probability distribution of maximum entropy subject to a set of constraints, which imply the expected value for each feature matching their empirical means represented by real-value variables (Phillips et al., Reference Phillips, Anderson and Schapire2006).
We used the MaxEnt algorithm to build a bioclimatic distribution model for the Endangered huemul Hippocamelus bisulcus (Jiménez et al., Reference Jiménez, Guineo, Corti, Smith, Flueck and Vila2008), which is endemic to the southern Andes of Chile and Argentina, and the most threatened of South America's 15 deer species (IUCN, 2013). This medium-sized deer is the only large herbivore inhabiting Patagonian forests. Although it inhabits a variety of environments, from mountain ranges > 2, 500 m (Povilitis, Reference Povilitis1986) to periglacial valleys at sea level in the Patagonian fjords (Frid, Reference Frid2001), its primary habitat is montane forest dominated by Nothofagus spp. on Andean slopes. The huemul population has declined by 99%, and its distribution by 50%, since the arrival of European settlers during the 19th century (Redford & Eisenberg, Reference Redford and Eisenberg1992; Vila et al., Reference Vila, López, Pastore, Faúndez and Serret2006).
Habitat loss through conversion of native forest into farmland, disease transmission (Corti et al., Reference Corti, Saucedo and Herrera2013), and competition with domestic livestock (Frid, Reference Frid2001; Vila et al., Reference Vila, Borrelli and Martínez2009), poaching, the introduction of exotic species, and predation by roaming dogs (Corti et al., Reference Corti, Wittmer and Festa-Bianchet2010) are assumed to be the main causes of the huemul's population decline. Additionally, Wittmer et al. (Reference Wittmer, Elbroch and Marshall2013) suggested that natural predators, such as pumas Puma concolor and culpeo foxes Lycalopex culpaeus, in areas with abundant alternative prey could also, at least locally, be an important cause of huemul decline. Furthermore, the lack of connectivity among the remaining small populations is considered one of the main threats to the species, with the potential risk of inbreeding and possible potential local extinctions (Corti et al., Reference Corti, Shafer, Coltman and Festa-Bianchet2011).
Although there has been research and conservation directed at the huemul since the 1970s, reliable information on the species’ population status and trends is still limited. Chile and Argentina have national conservation plans for the species (Ministerio de Medio Ambiente y Desarrollo Sustenable, 2002; CONAF, SAG, CONAMA, 2009), which reflect both countries’ efforts to coordinate conservation measures, and suggest areas for further research. Nevertheless, population-monitoring efforts have often had shortcomings in terms of funding, and lack of staff and methodology. Intrinsic biological characteristics of the species and difficult access to its habitat make field surveys a complex enterprise, and thus monitoring programmes often fail (Wittmer et al., Reference Wittmer, Corti, Saucedo and Galaz2010).
Post-glacial recolonization of the huemul's northern distribution range probably occurred from multiple refugia in the north-east Patagonian region, in what is now Los Alerces National Park in Argentina, and surrounding areas (Marín et al., Reference Marín, Varas, Vila, López, Orozco-terWengel and Corti2013). The remnant population in Nevados de Chillán, the northernmost extant huemul population, may have been derived from those geographically closest refugia in the North Patagonia region. In this area, huemul inhabit the Valdivian Ecoregion (Vila et al., Reference Vila, Saucedo, Aldridge, Ramilo, Corti, Barbanti Duarte and González2010), a Global 200 Ecoregion (Olson, Reference Olson1998; Lara et al., Reference Lara, Little, Urrutia, McPhee, Álvarez-Garretón and Oyarzún2009). Argentina's National Park Administration has systematically recorded direct and indirect signs of huemul presence within and around protected areas since 1960 (H. Pastore, pers. obs.). However, for neighbouring districts in Chile there is little information about huemul presence, although the species is likely to occupy similar habitat types to those in Argentina.
To build a species distribution model for the Valdivian Ecoregion we used data from areas where huemul are known to occur there (Marín et al., Reference Marín, Varas, Vila, López, Orozco-terWengel and Corti2013), and projected the model predictions northwards to areas where huemul have not been recorded in recent decades. We were thus able to identify potential distribution areas between the known northern Patagonian populations and a small remnant population in the northernmost distribution range, 400 km apart. The results facilitate identification of existing and potential biological corridors between protected areas and areas prioritized for landscape conservation, through which huemul dispersal could occur.
Study area
The study was conducted in the central part of the current huemul distribution (Fig. 1; 37–44°S), where the species’ historical habitat is characterized by the Valdivian Rainforest. The Valdivian Ecoregion covers c. 155,000 km2, overlapping part of the southern Andes of Argentina and Chile (WWF, 2012). Predominant forest types are evergreen, southern beech and cypress forests. The ecotones are dominated by grasses and shrubs, such as Embothrium sp., Maytenus sp., Chiliotrichum sp., Pernettya sp., Berberis spp., Escallonia sp. and Empetrum sp. The area is characterized by steep slopes, with elevations of 20–2,500 m.
Fig. 1 Projected potential distribution of huemul Hippocamelus bisulcus in the Valdivian Ecoregion of Chile and Argentina, predicted by species distribution modelling. Each cell represents 1 km2.
Methods
Presence data and predictor variables
The initial presence data set included 661 records. The National Park Administration in Argentina provided data from their database of huemul presence, which contained 600 records collected during 1960–2012 in four protected areas in Argentina: Los Alerces and Lago Puelo National Parks and neighbouring provincial protected areas (Fig. 1). To these data we added 32 presence records from outside the protected areas, also from the National Park Administration (Fig. 1). Most of the presence records in the National Park Administration's database were from Los Alerces National Park (442 records), where the Park Administration has conducted a huemul monitoring programme since 2001. This programme is based on six established transects that are surveyed annually during the southern spring and autumn. The database includes both direct sightings and indirect signs, such as footprints and faeces, identified by experienced park rangers. Additionally, during the summer (January–March) of 2012 we made 29 new observations in the same National Parks in Argentina (Los Alerces, n = 8; Lago Puelo, n = 12), and in the Futaleufú National Reserve in Chile (n = 9; Fig. 1), mainly footprints and faeces but also observation of five huemul in Futaleufú.
We used ArcGIS v. 9.3 (ESRI, Redlands, USA) to visualize the spatial distribution of huemul observations and their relation to landscape features (lakes, rivers, roads) and topography (elevation, slope). Presence data are more likely to be recorded near roads because these areas are more accessible to observers, and therefore we deleted from the database records that were linearly clustered along roads, to minimize sample bias in presence data and to maintain spatial independence among observation points. Thus our dataset was reduced from 661 to 396 records. To prevent spatial autocorrelation, MaxEnt retains only one presence record per 1 km2 cell (Tognelli et al., Reference Tognelli, Abba, Bender and Seitz2011). This further reduced the number of records to 258, of which 194 (75%) were used to build and train the models and 64 to validate them.
We consider potential predictor variables to be all environmental variables that could potentially influence the species and therefore can be used to predict its potential distribution. The modelling is based on the assumption that the probability of huemul presence is related to the combination of several of these variables in a specific condition at a given location, and attempts to identify the most relevant predictor variables and generate the best possible predictive model.
We used global datasets from open-access sources of bioclimatic, vegetation and topographical variables: (1) eighteen variables from the WorldClim database of interpolated climate surfaces at 1 km resolution, derived from monthly temperature and rainfall values (this database uses only climate databases with at least 10 years of data, and provides mean values for 1960–1990; Hijmans et al., Reference Hijmans, Cameron, Parra, Jones and Jarvis2005); (2) the vegetation variable from the Global Land Cover 2000 Project (GLC2000) at 1 km resolution, providing information about croplands, vegetation types, water bodies, artificial surfaces, bare areas and permanent snow and ice, using satellite images from 2009 (Bartholomé & Belward, Reference Bartholomé and Belward2005); and (3) slope and elevation variables from the Shuttle Radar Topography Mission at 90 m resolution (Farr et al., Reference Farr, Rosen, Caro, Crippen, Duren and Hensley2007). All layers were set to a common scale using ArcGIS, to ensure congruence of all input data (Hu & Jiang, Reference Hu and Jiang2010; Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011).
Using ArcGIS we created a grid of 1 km2 cells (Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011). Using the packages psych and GPArotation in R v. 2.15 (R Development Core Team, 2012) we identified those variables most strongly associated with the first and second components in a principal component analysis that included the 21 predictor variables for the 258 grid cells in which huemul were recorded. These two components together accounted for 83% of the variation in the 21 predictor variables (Table 1).
Table 1 Results of principal component analysis applied to 21 predictor variables for the 258 1 km2 grid cells in which huemul Hippocamelus bisulcus were recorded (see text for details), with the loading of each variable associated with the first three components, the proportion of the variance explained by each component, and the cumulative proportion. The highest loading scores (in bold) indicate the strongest relationship between variable and component (i.e. most of the temperature variables are strongly related to the first component, PC1, and the precipitation variables to the second component, PC2).
Modelling with MaxEnt
Of the bioclimatic variables that were highly associated with huemul presence (Table 1), only those pairs that did not exceed a pairwise Pearson correlation value of 0.70 (results not shown; Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011) were retained. Four models were run, with five temperature variables, two precipitation variables, vegetation, elevation and slope.
The four models were fitted with MaxEnt v. 3.3.3 (Phillips et al., Reference Phillips, Anderson and Schapire2006) using the following settings: regularization multiplier = 1; maximum number of iterations = 1,000; convergence threshold = 10−5; maximum number of background points = 10,000; and adjusted sample radius = −6. The models were run using auto-features, by which MaxEnt computes the default mathematical functions of the environmental features (i.e. linear, quadratic, product, threshold and hinge features). To construct a binary map from the MaxEnt outputs (i.e. differentiating only suitable from unsuitable cells) we used an equal training sensitivity and specificity threshold of P > 0.40 (Jiménez-Valverde & Lobo, Reference Jiménez-Valverde and Lobo2007; Freeman & Moisen, Reference Freeman and Moisen2008). This is a conservative approach that increases the specificity of the model, making its predictions more precise. To facilitate interpretation of the model we selected the logistic output format, which provides a proxy of probability of presence (Peterson et al., Reference Peterson, Soberón, Pearson, Anderson, Martínez-Meyer, Nakamura and Araújo2011). Jackknife analyses were carried out on the regularized gain of the training data to examine the relative importance of each predictor variable for model performance (Hu & Jiang, Reference Hu and Jiang2010). The remaining model training parameters were left at their default settings. We used the logistic outputs, given as probabilities representing degrees of habitat suitability, from 0 = unsuitable to 1 = most suitable habitat (Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011).
The models were trained in the geographical area where the presence records were recorded, and then using the projection tool in MaxEnt their predictions were projected into a larger area: approximately 37–44°S and 70–74°W. Test and difference areas under the curve were used as performance indicators to select the best-fit model (Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011; Warren & Seifert, Reference Warren and Seifert2011). We then used ArcGIS to calculate the potential distribution area of this model, and how much of it currently falls within the protected area system in each country. The shape files for protected areas were obtained from Protected Planet (Bertzky et al., Reference Bertzky, Corrigan, Kemsey, Kenney, Ravilious, Besançon and Burgess2012).
Results
Annual temperature range and precipitation were the predictor variables that contributed the most information about huemul distribution to the four models (Table 2). Summer precipitation was second in importance for model 4 (M4) and mean diurnal range for model 2 (M2). Other temperature variables (i.e. bio1, bio8 and bio10; Table 2) had only limited contributions in all models. We will focus on the predictions of M4 to describe the predicted potential distribution of huemul (Fig. 1) because it had the best performance according to the test and various values of area under the curve (Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011; Warren & Seifert, Reference Warren and Seifert2011). Of the potential distribution area in Argentina, c. 40% (2,759 km2) is within the protected area system, and in Chile < 25% (1,214 km2) of the potential distribution area is in the protected area system.
Table 2 Percentage contribution of predictor variables to models M1–M4, generated by the MaxEnt algorithm, and model performance. Variables that contributed more information, or with information not present in other variables, are highlighted in grey.
The response curves describe the variation of the potential distribution prediction based on its dependence on one predictive variable and on the dependencies induced by correlations between one variable and the other variables. The predicted potential distribution shows a peak at an annual temperature range of c. 22° (bio7; Table 2), decreasing rapidly for larger and smaller temperature ranges (Fig. 2a); the empirical mean of the presence data for this variable was 22.19 ± SD 0.7°C (n = 258). In the case of precipitation during the summer season (bio18; Table 2) there is a pronounced peak at c. 120 mm (empirical mean = 128.96 ± SD 23.17 mm; Fig. 2b). The model attributes little importance to vegetation cover (gcover: contribution 3.3%; Table 2). However, the model predicted zero probability of presence for almost all categories related to wetlands, artificial surfaces and most of the forest categories, 35% of the empirical data were located in the vegetation category ‘closed to open shrubland’, and the second most frequent category (c. 28%) was the vegetation category ‘closed to open broadleaved evergreen or semi-deciduous forest’. Other variables had lesser importance for model performance.
Fig. 2 Variation in the predicted probability of presence of huemul with (a) annual temperature range and (b) summer precipitation.
Discussion
Historically, huemul occurred in the Andean Cordillera from the Cachapoal River (c. 34°S) southwards to the northern shore of the Magellan Strait (54°S; Cabrera & Yepes, Reference Cabrera and Yepes1960). The current known distribution is in two areas separated by c. 400 km. The northern population is at Nevados de Chillán (Chile) and the southern is distributed among several fragmented areas extending southwards from Nahuel Huapi National Park (Argentina). Marín et al. (Reference Marín, Varas, Vila, López, Orozco-terWengel and Corti2013) suggested that the huemul at Nevados de Chillán are genetically derived from the Eastern Andes Refugium, in northern Patagonia, from where we collected our data. Our results show a wide range of potential distribution within the projection area, where there have been no records of huemul presence since the 1980s. The isolation of the northern population may imply that other factors are preventing the species from occupying the identified potential areas. These factors could include natural and artificial dispersal barriers, the distribution of native and exotic predators, and resource competition and interaction with domestic and exotic species (e.g. red deer Cervus elaphus; Corti et al., Reference Corti, Wittmer and Festa-Bianchet2010). Frid (Reference Frid2001) and Briceño et al. (Reference Briceño, Knapp, Silva, Paredes, Avendaño and Vargas2013) found evidence that huemul modify their habitat use when sharing habitat with exotic ungulates. Red deer are known to be present in the study area, and probably influence huemul presence (Jaksic, Reference Jaksic1998). Nonetheless, we limited our analysis to abiotic predictors because of the complexity of interpreting the relative importance of biotic and abiotic factors simultaneously (Guisan & Thuiller, Reference Guisan and Thuiller2005). Additionally, georeferenced information about such biotic factors is scarce for the area we studied.
Excluding the presence records linearly clustered along roads improved the predictive power of the model. When the whole dataset was used, part of the predicted area of potential occupancy comprised water bodies. This is because of the location of the main transportation routes, which increases the likelihood of observing huemul, thus biasing the analysis.
The area of potential distribution (P > 0.40) predicted by the best-fitting model (M4) is c. 12,360 km2, of which 7,000 km2 (57%) is in Argentina and 5,360 km2 (43%) is in Chile. Based on our modelling results, potential areas of distribution were found to be restricted to the Andes range, as expected, markedly decreasing towards the west, probably because of an increase in annual precipitation and humidity, which promotes the growth of dense evergreen forest. Towards the east, precipitation decreases gradually, resulting in the arid conditions of the Argentinian pampas. In Argentina suitable conditions for huemul exist in the southern part of Neuquén province and southwards. In Chile there is the opposite situation: there is little potential area for huemul distribution in the southernmost district but it increases towards the north as the administrative border between Chile and Argentina turns eastwards (Fig. 1).
The temperature suitability range predicted by our model seems narrow compared to the temperature gradients occupied by huemul across their whole range, but this is related to the elevation range used by the species in this part of its distribution. This particular combination of climatic variables has an indirect effect on huemul distribution, mediated through the type of vegetation that exists under these conditions. The dense forest type of the Valdivian Ecoregion permits only limited movement of medium to large ungulates, such as the huemul; it is the optimal habitat for pudu Pudu puda, a dwarf deer of southern South America (Jiménez, Reference Jiménez, Barbanti Duarte & and González2010).
Some presence records are from areas categorized by the model as only moderately suitable (i.e. some records from Futaleufú National Reserve in Chile, where PQ recorded several signs of huemul and directly observed five huemul during this study). This apparent incongruence between model predictions and field observations may imply that some individuals have dispersed to less suitable areas from more suitable areas nearby. However, further research on the dispersal dynamics of huemul populations in these areas is needed to confirm this, and would help us to understand the role and importance of protecting such moderately suitable areas for huemul conservation. It would also help to identify natural corridors between existing protected areas (Gilbert-Norton et al., Reference Gilbert-Norton, Wilson, Stevens and Beard2010).
Temporal differences between field records and environmental variables (especially land cover) used in this type of modelling should be as small as possible. However, because environmental conditions generally do not change rapidly, several years of difference may be acceptable, and historical presence records from museum and species collections may be used for species distribution modelling (Elith & Leathwick, Reference Elith and Leathwick2009). We used the complete huemul dataset provided by Argentina's National Park Administration because most of the presence records (92.9%) were post-1995 and we assumed that the environmental conditions had not changed substantially since then.
The low predictive performance of vegetation and topographical variables in the model may be attributable to the 1 km2 scale used; however, climatic and vegetation variables are usually highly correlated (Wintle et al., Reference Wintle, Elith and Potts2005). In addition, the maps of potential distribution predict suitable bioclimatic conditions but they cannot be used to infer causal effects of environmental variables on species’ distribution, for which other kinds of data, such as telemetry and presence–absence data, and other analytical tools would be needed.
Parts of the study area located on the watersheds between Chile and Argentina include protected areas and potential ecological corridors but are threatened by hydropower development projects (Urrutia et al., Reference Urrutia, Lara and Villalba2005). We recommend that future surveys in these areas use the predicted potential distribution presented here to determine huemul presence; our projections will enhance the information available to decision makers when evaluating the conservation risks associated with the implementation of energy, infrastructure and tourism projects. Conservation planning at a landscape scale should pay special attention to potential natural corridors, which connect areas of potential distribution with areas where huemul populations are known to exist. A network of public and private protected areas interconnected by biological corridors would provide protection and enhance huemul dispersal (Corti et al., Reference Corti, Shafer, Coltman and Festa-Bianchet2011), although huemul conservation on private properties also faces significant challenges (Wittmer et al., Reference Wittmer, Elbroch and Marshall2014).
Our best model suggests the potential presence of huemul in areas for which information on huemul presence is completely lacking since the 1980s. We recommend investigating potentially unrecorded populations and prioritizing survey efforts in these projected areas, which would provide an empirical validation of the model (Lobo et al., Reference Lobo, Jiménez-Valverde and Hortal2010). We need to be cautious, however, when interpreting the model's predictions. This is because of the potential for over-prediction, possibly associated with the choice of the projection area and with the use of a small part of the huemul's habitat range for the model calibration (Marino et al., Reference Marino, Bennett, Cossios, Iriarte, Lucherini and Pliscoff2011). It is also important to train models under various climate change scenarios to predict how the potential huemul distribution will be affected by increasing temperatures (Rebelo et al., Reference Rebelo, Tarroso and Jones2010).
Vila et al. (Reference Vila, López, Pastore, Faúndez and Serret2006) suggested a greater degree of fragmentation occurred in the Chilean portion of the huemul distribution in the Valdivian Ecoregion as a result of human disturbance. Accordingly, our results indicate that the proportion of the area potentially inhabited by huemul is larger in Argentina for this part of the original huemul distribution range. This suggests the need for a binational conservation strategy to coordinate the various stakeholders and to ensure that conservation measures will be implemented in both public and private protected areas throughout the entire range of the huemul's distribution. It is imperative that huemul monitoring continues in national protected areas, to verify the distribution predicted by the model and to check habitat conditions in these areas.
Acknowledgements
We thank Martin Wilmking for his comments and support during the development of this research, H. Wittmer and an anonymous reviewer for constructive comments, and the Rosa Luxemburg Foundation for funding a scholarship to P. Quevedo. We are grateful to Argentina's National Park Administration (APN) and Chile's National Forestry Corporation (CONAF) for collaborating in this research, providing historical data, and logistical and technical support during the field surveys conducted in 2013, and for their ongoing commitment to huemul conservation. We are particularly grateful to César Bastias, Néstor Diocaretz, Neftalí Rozas and Carlos Salvadores from CONAF, and Félix Vidoz, Mauricio Berardi and Javier Montbrun from APN, for assistance during the fieldwork. We also thank Evaristo Araneda for sharing his field knowledge, Marcela Morales and Johannes Horstmann for their assistance during fieldwork, and Ignacio Díaz and Álvaro González for their advice on GIS. Achaz von Hardenberg was financially supported by a Research Education Grant from the European Social Fund of the European Union, the Autonomous Region Aosta Valley and the Italian Ministry of Work and Social Providence. P. Corti was funded by FONDECYT 3110187 and CONICYT 7912010016.
Biographical sketches
Paloma Quevedo is a veterinarian with an interest in habitat modelling for species with conservation problems. She uses geographical information system technologies to identify innovative options for wildlife management. Achaz von Hardenberg is a biologist with an interest in the application of statistical modelling tools in conservation ecology. His research mostly focuses on the population biology and behavioural ecology of mountain wildlife. Hernán Pastore is a biologist interested in management and conservation of biodiversity, and related policies. José Álvarez is a forestry engineer working on the conservation of threatened species and landscapes. Paulo Corti is a veterinarian focused on the behaviour, genetics and demographic variables affecting evolutionary processes in wild populations of mammals and birds, and the application of this information in conservation and management plans; he is a member of the IUCN Deer Specialist Group.
