Introduction
Species that are faced with an imminent risk of extinction are conservation priorities (Vane-Wright Reference Vane-Wright, Kellert and Speth2009, IUCN 2012). Such species are typically rare and low in numbers (Gaston Reference Gaston1994), which makes them difficult to detect and as a result, it is difficult to assess their status with confidence. Consequently, assessments are often reliant on subjective expert opinion (van der Ree and McCarthy Reference Van der Ree and McCarthy2005). The monophyletic Himalayan Quail Ophyrisa superciliosa is ‘Critically Endangered’ (IUCN 2012) due to a small putative population, a narrow geographic distribution and intensive habitat modifications (IUCN 2012). Quantitative assessment of the species’ status and potential distribution is essential for the targeting of efforts to rediscover the species and hence its conservation.
Few records of the Himalayan Quail exist: the last specimens date from 137 years ago and there has been a lack of confirmed records since. Re-sighting probability calculations offer a more objective, probabilistic insight into extinction assessments as we have little evidence for, or confidence in, declarations of extinction based on the raw data alone (Butchart et al. Reference Butchart, Stattersfield and Brooks2006). A recent estimate for Himalayan Quail suggests it went extinct in the late 1890s (Collen et al. Reference Collen, Purvis and Mace2010), only c.20 years after the last of a small number of sightings. Recently, further reports have been collated and all records for the species have been subjected to a critical re-evaluation (Boakes et al. Reference Boakes, McGowan, Fuller, Ding, Clark, O’Connor and Mace2010). Consequently, a revised estimate of extinction is needed. If the species is extant, search effort needs to be targeted at the most suitable areas.
Specimen records for the Himalayan Quail come from Mussoorie and Nainital, in the state of Uttarachand, Northern India. Surveys undertaken to date have failed to detect the Himalayan Quail around these areas and nor have further surveys in western areas of Nepal (R. J. Cuthbert unpubl. data). The lack of modern records suggests the need for enhanced techniques in these areas and/or searches in additional areas where the species may still occur. However, the distribution of the species is uncertain and its geographic range is all but unknown or inaccurately mapped (e.g. BirdLife International 2014a). An updated and quantitative assessment of its range made using the best available data is needed if surveys are to be appropriately targeted.
For most rare species, surveys are undertaken in areas of the most suitable habitat to maximise the likelihood of detection. One approach to identifying such habitat is to generate species distribution models (e.g. Guisan and Zimmermann Reference Guisan and Zimmermann2000, Boitani et al. Reference Boitani, Maiorano, Baisero, Falcucci, Visconti and Rondinini2011). However, while historical records may still reflect the climate space of the species they are unlikely to reflect its habitat preferences because reliable records come from historical areas that have been extensively modified, creating a temporal mismatch between specimen records and habitat covariates. Other approaches are therefore needed for this species. The use of proxy species is well documented in conservation biology (Caro et al. Reference Caro, Eadie and Sih2005) but the use of proxy species to identify suitable habitat is novel. Here we use more abundant and widespread species that might share similar habitat requirements as a guide for directing survey efforts. Cheer Pheasant Catreus wallichii and Himalayan monal Lophophorus impejanus have distributions that overlap the known range of the Himalayan Quail and utilise areas of dense grass that are potentially similar to the habitat requirements of the Himalayan Quail (Kaul Reference Kaul1992, Del Hoyo et al. Reference Del Hoyo, Elliott and Sargatal2001, BirdLife International 2014a). Although extensive knowledge of the Himalayan Quail’s habitat preferences is absent, the proxy species used to identify suitable habitat must be well justified, especially if they occupy a much wider geographic distribution than that of the target species.
We attempt to update our knowledge of the potential distribution of the Himalayan Quail by generating a climatic and topographic model for the Himalayan Quail, the Himalayan Monal and Cheer Pheasant using maximum entropy (Maxent) species distribution modelling software (Phillips et al. Reference Phillips, Anderson and Schapire2006, Phillips and Dudík Reference Phillips and Dudík2008). Generating such a model for the Himalayan Quail is useful because these parameters are less likely to have changed than land cover and are therefore likely to represent a realistic niche for this species. We next test the overlap of the Himalayan Monal and Cheer Pheasant climate models with that of the Himalayan Quail to assess whether these two species can be used as potential proxies. Overlap between the climatic and topographic model of the Himalayan Quail and those of the proxy species would indicate that they occupy similar climate spaces and thus reflect the suitability of the proxy species as surrogates. Subsequently, we generate full species distribution models including habitat for the proxy species, which are reliable because there are modern records for these species. Overlap between the full habitat models for the proxy species and the climate model for Himalayan Quail will indicate potentially suitable areas for the species. We refine this area of suitable habitat further by incorporating measures of previous survey efforts to identify areas that have been poorly searched.
Methods
Bird records
We took data from the GALLIFORM: Eurasian Database V.10 (Boakes et al. Reference Boakes, McGowan, Fuller, Ding, Clark, O’Connor and Mace2010) that contained point locality data accurate from 0.62–30 miles (1–48.3 km; for further details see Tables S1-S3 in the online supplementary material), collected from a wide-range of sources including museum specimens, ringing records, biological atlas data and trip reports (Figure 1a and Table S1).
We omitted 19 records that lacked latitude and longitude coordinates and a date. Given the uncertainty around the reliability of the contemporary records, we also undertook a preliminary vetting procedure to ensure that the Himalayan Quail data used in the creation of our species distribution models were suitable (two more records were omitted in this way; for more details of the number and nature of records considered, see Table S4). All other records were used along with the historical specimen records in the climate model. This reduced our sample size for Himalayan Quail from 55 to 34 records.
We used all available post-1980 records in the Cheer Pheasant and Himalayan Monal models, as these species are well known and readily identified. The Cheer Pheasant has already been identified as a potential proxy for the Himalayan Quail (Kaul Reference Kaul1992) but the use of the monal in this way is new and is based upon field insights into its habitat use (R. J. Cuthbert unpubl. data). In contrast to the Cheer Pheasant, the monal is known to be an altitudinal migrant (BirdLife International 2014b), leading us to believe that the summer distribution of the monal might be most similar to that of the Himalayan Quail. Accordingly we first attempted to model the summer distribution of the monal but also modelled the annual distribution of the monal in case the resulting summer climate model overlapped poorly with that of the Himalayan Quail. Summer was defined two ways: 1) on the basis of elevation i.e. records below 3,800 m (the mid-point of the Himalayan Monal’s altitudinal range) were classified as winter records (n = 249) and those above 3,800 m (n = 68) were classified as summer records (Ali and Ripley Reference Ali and Ripley1983); 2) on the basis of specimen labels recorded as summer (n = 12) or winter (n = 54). We used both approaches because although the former was more subjective, not all records were labelled according to the season when they were collected and we wanted to include as much point locality data in our models as possible to ensure they gave both accurate and precise predictions.
Calculating the likelihood of extinction
We used Optimal Linear Estimation (OLE; Cooke Reference Cooke1980) to assess the probability that the Himalayan Quail was globally extinct. OLE is a technique that is commonly used to assess the likelihood of extirpation (Roberts and Solow Reference Roberts and Solow2003, Solow Reference Solow2005). OLE is a non-parametric method that allows extinction dates to be estimated based on the distribution of the most recent sightings. The main assumption is that the sighting effort never falls to zero over an annual time step, particularly around the time of extinction. OLE is known to be sensitive to the number of records included in the calculation, but for species with more than five sightings, it is thought to provide robust estimates for the time of extinction (Collen et al. Reference Collen, Purvis and Mace2010). To ensure our estimates were robust, we included all records, records from the last five years before (and including) the final sighting and specimen records only (see Table S5). Calculations were undertaken in R (R Development Core Team 2012) using the R package ‘sExtinct’ (Clements Reference Clements2013).
Modelling approach for climate envelopes and species distribution models
We created two sets of species distribution models: climatic and topographic models for the Himalayan Quail and proxy species and then full models also including land cover for the proxy species only. We used Maxent (version 3.3.3k; Phillips et al. Reference Phillips, Anderson and Schapire2006) to model the likelihood of occurrence of the species using presence locations of each species in turn as a function of topography and climate for the first set of models and topography, climate and land cover for the second set of models. Maxent has been found to perform well against other distribution models (Elith and Graham Reference Elith and Graham2009) and produces models that have particularly high accuracy in the case of species with small sample sizes and restricted geographic locations. Rather than use Maxent’s default setting of selecting 10,000 randomly generated pseudo-absences, we used locations from which there were records from other Galliformes (n = 820 records) to generate a targeted set of pseudo-absences. Thus, our pseudo-absences were chosen from sites with the same sampling bias as the presence sites for a suite of species that were observed with similar sampling techniques. This “target group” approach (Phillips and Dudik Reference Phillips and Dudík2008) reduces the potential for species distribution model outputs to be affected by sampling biases in study species records in both time and space (see Boakes et al. Reference Boakes, McGowan, Fuller, Ding, Clark, O’Connor and Mace2010 for a description of the sampling biases in the dataset that we use).
For each species, we omitted duplicate records to identify a subset suitable for inclusion in our final species distribution models. Analysis was restricted to occupied and immediately surrounding ecoregions for each modelled species. By restricting the analysis in this way, we ensure that we obtain outputs that are sensible in their geographic distribution (for further details of the Himalayan ecoregions used in our analysis and the distribution of records see Table S6).
Climate variables were downloaded from www.worldclim.org/bioclim (Hijmans et al. Reference Hijmans, Cameron, Parra, Jones and Jarvis2005) and were approximately 1 km2 in resolution. In an effort to reduce the number of variables considered (and thus reduce the risk of model over-fitting), we only considered four of the 17 potential climate variables relating to temperature and precipitation that described the major axes that are likely to affect the distribution of our three bird species. These were annual mean temperature, temperature seasonality, average precipitation and precipitation seasonality. Elevation was downloaded from the 90 m Shuttle Radar Topography Mission (SRTM) at 30 arc seconds (Jarvis et al. Reference Jarvis, Reuter, Nelson and Guevara2008) and was used to assess altitude, slope and aspect. Slope and aspect rasters were created and standardised in ArcMap version 10.2.
For the full proxy species distribution models we also incorporated a measure of Normalised Difference Vegetation Index (NDVI), a continuous measure of habitat type. NDVI data collected by the SPOT-Vegetation platform from 1999–2007 were downloaded from www.vito.be and the middle dekad (10-day period) of each month was extracted and averaged across years using raster calculator in ArcGIS version 10.2 and clipped to the modelled area using VGT-Extract. The rasters were standardised to the same resolution (1 km2) as before. No attempts were made to omit collinear variables as machine learning methods have been shown to still perform well with such variables, especially when the study goal is predictive accuracy (Elith et al. Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011).
Variables selected for inclusion in the final proxy models were those that contributed > 3% to the maximal model to avoid over-fitting the models while maximising their predictive power (this threshold was chosen after trial and error; for further model details see Table S7). Optimal feature functions were chosen based on sample size (Phillips and Dudík Reference Phillips and Dudík2008) and regularisation parameters were chosen based on AICc (Warren and Seifert Reference Warren and Seifert2011; for further details see Tables S8 and S9). All final models (see Table S10) were clipped to the occupied and neighbouring ecoregions of Himalayan Quail to maintain a focus on the Himalayan Quail climate map.
The ability of each model to discriminate between occupied and unoccupied areas was estimated from the area under the curve (AUC) of the receiver operating characteristics (ROC) (Phillips et al. Reference Phillips, Anderson and Schapire2006). We used cross-validations to generate 10 folds of randomly-selected presence data and ran the model 10 times, excluding each fold in turn and using the fold to validate the model (Phillips and Dudík Reference Phillips and Dudík2008). This uses all the data for validation and allows the standard deviation of the mean AUC to be assessed.
A minority of the point locality data were collected at a larger spatial resolution than our predictor variables (Tables S1-S3) which might potentially bias our outputs. However, the use of data with 1 km2 accuracy did not change the high AUC value for our SDMs (the difference in mean AUC values for models based on all data vs. models based on accurate data only were + 0.01 for the full model of the Cheer Pheasant, 0.00 for the full model of the Himalayan Monal and 0.00 for the climate model of the Himalayan Quail).
Comparing climate envelopes
The continuous model outputs for each species were categorised as suitable or unsuitable based on a threshold derived from the equal training sensitivity and specificity score (Tables S9 and S10). This threshold was chosen because it has been shown to minimise the rate of false positives and negatives (Pearson et al. Reference Pearson, Dawson and Liu2004).
We compared climate space between the Himalayan Quail and proxy species by identifying thresholded climate space for Himalayan Quail from the western part of the range and adding a minimum convex polygon to create a comparison area of approximately 20,000 km2. This was done to ensure a focus on the Himalayan Quail’s climate space and to ensure a fair comparison was made. We then clipped our Himalayan Monal and Cheer Pheasant climate outputs to this area and calculated a Spearman rank correlation for the continuous logistic output using ENMTools and the Kappa statistic, a measure of spatial agreement for the categorical thresholded output using Map Comparison Kit 3.2.3 (RIKS BV 2013). Finally, we assessed whether the Himalayan Quail was restricted in its climatic envelope by examining the shape of the modelled response curves.
Identifying overlap between climate envelope and species distribution models
The climate map of the Himalayan Quail and the full ENM maps for the proxy species were combined using raster calculator in ArcGIS 10 version 10.2. The maps were projected using a South Asia Albers equal area projection. Areas of suitable climate for the Himalayan Quail that overlapped with suitable habitat for the respective proxy species indicated potential areas for surveys.
Additionally, areas that overlapped with the proxy species were further refined by a measure of search effort. This was undertaken by creating a kernel density raster based on the number of post-1980 records of all galliform species in the area per square kilometre and divided into five geometric intervals. The localities previously identified by the overlap analysis were then multiplied by this raster to further refine priority sites. Thus, of the suitable areas based on climate and habitat suitability, those with the lowest search effort may be those in highest need of surveying.
Results
When all reliable records were used, our Optimal Linear Estimation calculation estimated the time of extinction to be the year 2023 (CI 1999-2120; Figure 1b and Table S5).
Based on the four climate variables used, Himalayan Quail had the smallest modelled climate distribution, followed by Cheer Pheasant and Himalayan Monal (Figure 2a). Variation in temperature and precipitation best predicted Cheer Pheasant potential occupancy; for Himalayan Monal it was highest based on variation in temperature, elevation, mean annual temperature, mean annual precipitation, temperature variation and slope; for Himalayan Quail it was highest based on variation in temperature, mean annual temperature and elevation (Table S11). The Himalayan Monal climate envelope created from all records was more tightly associated with the Himalayan Quail’s climate envelope than the Himalayan Monal climate envelopes based on summer records (see Table 1), so all records were used for the Himalayan Monal’s final species distribution model.
The likelihood of occupancy from the Cheer Pheasant and Himalayan Monal climate model was positively correlated to that of the Himalayan Quail’s climate model (Table 1) suggesting the likelihood of occupancy in shared locations increased in a similar way for both species. The Kappa statistics revealed slight (< 0.2) spatial agreement between the Himalayan Quail and Cheer Pheasant and moderate (0.4–0.6) spatial agreement between the Himalayan Quail and Himalayan Monal model. Taken together, Himalayan Monal appears to be a better proxy for Himalayan Quail by way of its climate envelope than Cheer Pheasant.
Creating species distribution models
We produced species distribution models (Figures 2b and 2c) for the two proxy species that incorporated additional land cover variables (NDVI). They represented suitable areas of 104,228 km2 for Cheer Pheasant and 162,249 km2 for Himalayan Monal. In comparison to the climate/topography models, this corresponded to an increase in suitable area of 1.2% and 150% for Cheer Pheasant and Himalayan Monal respectively. The mean AUC for Cheer Pheasant was 0.89 (SD = 0.02) and for Himalayan Monal was 0.80 (SD = 0.06). For full details of the relative importance of model covariates see Table S12.
We produced a combined map (Figure 3a) to show locations where the proxy species’ SDMs overlapped with the climate envelope of the Himalayan Quail (an area of 8,607 km2). Of the total Himalayan Quail cells with a modelled value exceeding the logistic threshold of equal training sensitivity and specificity (9,734 km2), 43.5% overlapped with both proxy species, 1.6% with Himalayan Monal only, 43.3% with Cheer Pheasant only and 11.6% had no overlap with either of the proxy species.
Of the localities that overlapped with both proxy species, 924 km2 (23%) were in areas with low levels of survey effort (6–40 records/km2), particularly in north-west India near to, but outside of, Mussoorie, Uttarachand and next to the Nepalese border (Figure 3b). As a starting point, we subsequently identified five large patches of suitable habitat (Table 2) that we feel should be targets for survey efforts based on these outputs and have included in the supplementary materials scalable maps in .kml format, which can be examined with Google Earth. By producing maps based on both proxy species and also on individual proxy species, we have given field scientists the flexibility to choose where to survey based on their own expert knowledge.
Discussion
The extinction status and distribution of the Himalayan Quail are uncertain. Clarification of both is urgently needed to determine the most appropriate conservation action for the species. We found that the species may still be extant, albeit with considerable uncertainty surrounding the precise extinction date. This may be because the OLE technique assumes search effort never falls to zero for each annual time step, particularly around the time of extinction (Collen et al. Reference Collen, Purvis and Mace2010). In our study, this assumption is often violated (Figure S1), possibly extending the upper confidence interval of our extinction estimates. However, even if this is true, it is still important to search for the Himalayan Quail to confirm extinction as the costs of incorrectly declaring extinction can be large. For example, giving up prematurely may doom the species to extinction (the ‘Romeo error’; Collar Reference Collar1998) and re-appearances (the ‘Lazarus effect’; Keith and Burgman Reference Keith and Burgman2004) may waste conservation resources if costly and extensive surveys are undertaken.
Having established there is reason to consider the species as extant, and hence there is the possibility that it will be rediscovered, we need to identify the potential distribution. However, identifying such habitat through species distribution models is problematic due to data paucity and a temporal mismatch between reliable specimen records and the habitat covariates available for analysis. We identified the habitat requirements for Cheer Pheasant and Himalayan Monal, which taken together should encapsulate those of the Himalayan Quail. Our results suggest that areas other than the Indian localities of Mussoorie and Nainital should be searched and that western Nepal appears less likely to contain suitable habitat for the Himalayan Quail, although there is a large area of suitable habitat on the Indian side of the Nepalese border. We identified new areas based on our models that have a high potential occupancy likelihood and where intensive survey efforts have yet to be applied. There is an urgent need to conduct better surveys in some parts of the Himalaya, a fact that is reinforced by the recent discovery of two new bird species in the Eastern Himalaya (WWF 2009).
While proxy species have been used relatively extensively in conservation biology (Caro et al. Reference Caro, Eadie and Sih2005), they have not yet been used to potentially delimit areas of suitable habitat for rare species as far as the authors of this paper are aware. We consider the approach undertaken in this paper to be applicable for identifying macro-habitat preferences for other species, where data are limited and proxy species with similar ecological requirements and habitat use are available. While this method may be of use for identifying broad areas where species are likely to occur, it is unlikely that we can use the habitat preferences of the Cheer Pheasant and Himalayan Monal to provide information on the Himalayan Quail’s micro-habitat preferences. Even at the broadest scale caution must be used given that Himalayan Monal are altitudinal migrants: in the summer they appear on edges of alpine pastures near treeline, but in winter, are driven down by snow to mid-altitude forests (Del Hoyo et al. Reference Del Hoyo, Elliott and Sargatal2001). Conversely, the available evidence suggests that this quail species is not an altitudinal migrant (Das Reference Das1995, Talwar Reference Talwar1995). Thus, there could be a potential difference in habitat use in the winter, even if the Himalayan Monal and Himalayan Quail inhabit the same altitude. Similarly, Cheer Pheasant might occupy shorter, less dense grassland at the microhabitat scale than the Himalayan Quail. As a result, it is difficult to recommend discounting suitable areas that have been better surveyed entirely, as previous searches might have been in the wrong microhabitats given the likely specificity of the Himalayan Quail’s habitat preferences.
In addition to these microhabitat considerations, there are also a number of specific methodological caveats for our analysis. First, it is possible that the Himalayan Quail had a large historical distribution in the lowlands, but had been forced upwards into the mountains by historical habitat conversion (Rieger and Waltzhony Reference Rieger and Waltzhony1993). If so, our species distribution models might be incorrect, giving a false impression of the quail’s preferred climate envelope. However, the data we have do not support this hypothesis with no records of any kind from the lowlands. Second, it is highly likely that there is heterogeneity in detectability for each species depending on habitat. Investigating this further through surveys is important, as species occurrences and recorder biases could vary with covariates in unexpected ways (Yackulic et al. Reference Yackulic, Chandler, Zipkin, Royle, Nichols, Grant and Veran2013).
A potential weakness in our approach was the possible effect of locational error from some of our point locality data on our modelling procedure. These locational errors could have made our models less accurate, both in terms of habitat associations and the resulting species distribution maps. In order to investigate this, we rebuilt our models based on a subset of data that all had the same locational error of 1 km2. We found that the resulting AUC values were similarly high, suggesting that locational error had not compromised the accuracy of our models, which is reinforced by the fact that we did not use the Maxent default setting of 10,000 pseudo-absences in our modelling procedure. Similarly, given that the majority of records used in each of our models are accurate to 1 km2 resolution and that the underlying land cover, topography and climate data were not at a particularly coarse resolution, we believe our approach is robust.
An alternative approach to the one outlined in this paper could be to use a rule-based assessment to identify potentially suitable habitat as has been used previously for both bird and mammal species (Hansen et al. Reference Hansen, McComb, Vega, Raphael and Hunter1995, Knick and Dyer Reference Knick and Dyer1997). For example, we could identify the Himalayan Quail’s altitudinal range from suitable historical records and use descriptions of long grass being a preferred habitat to identify a subset of localities for searches. However, this approach may be difficult to implement and may provide biased maps. This is because grass is included as a broad habitat class in many typical land cover maps such as the GLC2000 map (Bartholomé and Belward Reference Bartholomé and Belward2005), making it difficult to identify the long grass the Himalayan Quail is thought to prefer and which allows for potential confusion with agricultural land (Fritz and See Reference Fritz and See2008). As a result, we feel that our approach that makes use of NDVI data is likely to offer a better solution.
Our results indicate there is a large potential area of suitable habitat that has had little search effort applied to it, making ground-truthing some of the areas we have identified in our results the next logical step. To aid and encourage this, we include scalable .kml maps of our results in the supplementary material that identify potential areas for searching to enable field scientists and bird-watchers to plan effectively. As a starting point, we also include a list of five top priority localities with latitude and longitude information that should be surveyed first.
As our results indicate similarities in macro-habitat rather than micro-habitat, we recommend the following should be conducted prior to field surveys: 1) interviews with local herders, hunters, villagers and the State Forest Department staff, 2) a poster-mediated information campaign. Once this information is available, the exact nature of the field-surveys must account for local context and any searches for the Himalayan Quail in suitable habitat must be undertaken with the collaboration and cooperation of local communities and landowners. Field surveys could potentially include the use of camera traps, trained dogs and audio-surveys/playback in order to try to find this cryptic and hard to flush species. We recognise that continued, unsuccessful searching for the Himalayan Quail may eventually become cost-ineffective. As such, we suggest that the decision-theory framework provided by Chadès et al. (Reference Chadès, McDonald-Madden, McCarthy, Wintle, Linkie and Possingham2008) is used to determine how best to allocate limited conservation resources in the face of continued uncertainty and the authors are willing to provide analytical support to facilitate this. Once searches have been conducted, there is the potential to refine our species distribution models using the resulting presence/absence survey data (a similar approach has been taken with Gurney’s Pitta Pitta gurneyi in Myanmar; Donald et al. Reference Donald, Aratrakorn, Htun, Eames, Hla, Thunhikorn, Sribua-Rod, Tinun, Aung, Zaw and Buchanan2009).
In conclusion, the new data, while imperfect, suggest the Himalayan Quail may remain extant despite a lack of confirmed sightings for over 130 years. By using the habitat distributions of the Cheer Pheasant and Himalayan Monal as proxies, we identify new areas that could be used to guide search efforts on the macro-habitat scale, circumventing the minimal raw data available for the Himalayan Quail. While robust, our maps are built on some questionable assumptions, which require further testing. Overall, this is a novel approach to the problem of identifying survey areas for critically endangered species.
Supplementary Material
The supplementary materials for this article can be found at journals.cambridge.org/bci
Acknowledgements
Many thanks go to Dr Matthew Grainger for helpful advice and discussion and to the anonymous referees whose comments greatly improved this manuscript. This work was supported by National Environmental Research Council grant number 3000021024.