Combining radio-telemetry and random observations to model the habitat of Near Threatened Caucasian grouse Tetrao mlokosiewiczi

Abstract

The distribution of the Near Threatened Caucasian grouse Tetrao mlokosiewiczi, endemic to the Caucasus, was examined to model the species’ nesting habitat, and thus facilitate its conservation and the identification of Key Biodiversity Areas in the Caucasus. The species’ occurrence was defined by field surveys and radio-telemetry. Data were managed and analysed using a geographical information system and various modelling techniques. Grouse locations were divided into training and testing datasets. Habitat variables measured at training locations were used to develop models, and testing locations were used to validate the models. The final best-fit model suggested that Caucasian grouse prefer open habitat, and the most important independent variables accounting for the species' distribution were annual mean temperature, mean temperature of warmest quarter, precipitation seasonality and proximity to deciduous broad-leaf forest. The incorporation of human disturbance and ruggedness into the final model significantly increased its predictive power. This model provides a tool to improve search effectiveness for Caucasian grouse in the Caucasus and for the conservation and management of the species. The model can predict the probable distribution of Caucasian grouse and the corridors between known populations. Threatened and endemic species are often used as species for setting site-based conservation priorities, and this habitat model could help to identify new Key Biodiversity Areas for protection in the Caucasus. The Ministry of Environmental Protection and Natural Resources of Georgia is going to use the results of this study to reshape existing protected areas and identify new ones.

Introduction

Researchers and natural resources managers have long recognized that the best way to save rare and threatened species is to protect the places where they live. Protection or restoration of wildlife communities and their interactions requires prioritization of species (Langhammer et al., 2007). Such species are typically threatened and/or have a restricted range. Focusing on the protection of Key Biodiversity Areas is one effective way to conserve biodiversity globally (Myers et al., 2000). One way to identify such Areas is to model and overlay habitats of threatened and restricted-range species. Prioritization of species or wildlife communities for protection is required because resources available for conservation are limited and should therefore be invested in strategic ways to ensure the greatest contribution to preserving global biodiversity (Pressey et al., 1993).

The Caucasian grouse Tetrao mlokosiewiczi is endemic to the Caucasus and categorized as Near Threatened on the IUCN Red List (BirdLife International, 2008). It has the smallest range of any grouse species (c. 12,000 km²) and this range is highly fragmented and confined to the Greater and Lesser Caucasus mountains. Based on surveys and some habitat mapping the total population is believed to be 68,000–84,300 (BirdLife International, 2008) across Georgia, the Russian Federation, Turkey, Azerbaijan, Armenia and Iran. Georgia and the Russian Federation have the largest proportion of the species’ range and Iran supports the smallest fragment of the range, just to the south of Armenia. Gokhelashvili et al. (2003) summarize previous publications on the species’ biology. It is generally found at or near the tree line in meadows on slopes with Rhododendron thicket and sparse birch Betula forest at elevations of 2,000–3,300 m. Habitat loss and transformation are thought to be the major threats to the species, with many of the subalpine meadows within its range being used for intensive grazing (BirdLife International, 2008). In the Lesser Caucasus, gentler terrain, road construction and changes in land use provide relatively easy access for developers and hunters, contributing to increased habitat disturbance and degradation.

The Caucasian grouse is a key habitat quality indicator species that can be used to monitor the effectiveness of conservation actions. The species is sensitive to deteriorating habitat and to the quality of montane forests and subalpine and alpine ecosystems (particularly Rhododendron cover endemic to the Caucasus). Many other species of conservation concern could benefit from protection of the quality of the grouse’s habitat. The grouse could also serve as a flagship species for use in environmental awareness-raising and for generating support for conservation amongst the public, and the species has economic value because it attracts foreign birdwatchers and this benefits local economies.

In 2003–2005 large-scale conservation projects, funded by BP, were implemented by the Georgian Center for the Conservation of Wildlife (GCCW) and Doğa Derneği (the Turkish Nature Society) to improve knowledge of the species’ biology, develop monitoring and management activities and promote public awareness (BP, 2006a,b). In 2002 the Azerbaijan Ornithological Society and the German BirdLife partner NABU (Nature and Biodiversity conservation Union) started similar work on the species, financed by NABU, the UK Royal Society for the Protection of Birds, and the Dutch Ministry of Environment through the BirdLife European Division. Under these projects most of the species' range was surveyed and the presence of the grouse was recorded in many areas, some previously unknown. However, political unrest in the Caucasus and remote and rugged terrain hampered the collection of more recent data on the species' distribution and abundance and on threats. In these projects future work to develop a conservation strategy and create a potential distribution map for all range countries was planned. Gottschalk et al. (2007) developed a grouse habitat model for Turkey.

The species’ Near Threatened status, the lack of an accurate standardized habitat model for the whole species’ range, and the recent discovery of the species in new areas motivated us to improve modelling techniques for the species and examine its habitat selection. Our objectives were to predict the global distribution of Caucasian grouse and estimate its population size based on knowledge of the species' density. Because search for, and census of, Caucasian grouse is difficult over rugged and climatically inclement areas, our model will allow detection and population estimates to be more efficient and will facilitate planning and zoning efforts in national park systems in the Caucasus.

Methods

Sampling

Our study area was the entire Caucasus. We captured eight adult males and six adult females at active leks in Lagodekhi National Park and Erusheti Key Biodiversity Area (which lie in Georgia’s Greater and Lesser Caucasus, respectively) during the breeding season (April–June) in 2004–2005 using walk-in traps (Schroeder & Braun, 1991) and rubberized leg-hold traps originally designed for martens Martes martes and Martes foina. Each captured bird was fitted with a 10.7-g necklace-style radio-transmitter (A3950 transmitters; Advanced Telemetry Systems, Isanti, USA). These individuals were located every month during 2004–2006, except for January–March (because of deep snow cover), using a 3-element folding Yagi antenna and an FM100 scanning telemetry receiver (Advanced Telemetry Systems, Isanti, USA). Standard telemetry techniques were used to determine locations of the birds. Locations were recorded, using a global positioning system and a rangefinder, when birds could be flushed or directly observed using a 22 × telescope.

To avoid the repeated sampling of habitat variables (see below) we used presence locations that were > 1,415 m from neighbouring locations (our analyses were performed on habitat variable grids of 1,000 × 1,000 m cells). We did not use presence locations recorded at leks because a study based on these locations (A. Gavashelishvili, unpubl. data) had provided a model that was appropriate only for mapping leks and produced meaningless results for habitat where grouse occur for most of the year. Radio-telemetry provided 200 presence locations. To include a wider geographical range we also used 150 additional presence locations that we had recorded during wildlife surveys in various parts of Georgia and the Caucasian part of the Russian Federation in 2004–2005 and from Armenia and Turkey (see Acknowledgements). Thus we obtained a total of 350 presence locations (Fig. 1) based on observation of signs (faeces, feathers, snow burrows, egg shells) and on sightings. Grouse faeces were identified from their mtDNA. Data were mapped and managed using the geographical information system (GIS) ArcGIS v. 9.2 (ESRI, Redlands, USA).

Fig. 1 (a) The training locations obtained to develop Caucasian grouse Tetrao mlokosiewiczi models, and (b) the test locations used to validate the models. +, presence; •, absence.

In addition to these 350 presence locations we obtained additional presence locations and absence locations from various sources (see below). We used the 350 presence locations and some of the absence locations to fit models (training), and the other presence and absence locations to evaluate the predictive accuracy of the models (testing).

Habitat variables

We considered bioclimatic variables related to climate, terrain, human disturbance and land cover (Table 1). Extracting terrain, anthropogenic and land cover data from topographic maps and atlases of large areas is time-consuming and expensive. Instead, we used free online digitized data and managed them using ArcGIS and ERDAS IMAGINE v. 9.1 (Leica Geosystems Geospatial Imaging, Norcross, USA). The variables used were based on documented species-habitat associations (Averin, 1938; Kutubidze, 1961; Tkachenko, 1966; Sikharulidze, 1974; Potapov, 1985; Atkinson et al., 1995) and our field experience, and with regard to their availability. Before analyses all data were re-projected into a world sinusoidal projection and re-sampled to a 1,000-m cell size, using the bilinear interpolation for continuous data and the nearest neighbour method for categorical data. We obtained climatic data from WorldClim v. 1.4 (WorldClim, 2010). This is a set of global climate layers with a spatial resolution of 1-km cells (Hijmans et al., 2005; Table 1). Terrain data were measured from the Shuttle Radar Topography Mission (SRTM) elevation data (Global Land Cover Facility, 2010) at a resolution of 1-km cells.

Table 1 Bioclimatic variables used for modelling Caucasian grouse Tetrao mlokosiewiczi habitat throughout the Caucasus.

For land cover data we processed the MODIS Land Cover Yearly L3 Global 500 m dataset coded MCD12Q1 (EOS Data Gateway, 2010). From this dataset we used the primary land cover scheme that identifies 17 land cover classes defined by the International Geosphere Biosphere Programme, which includes 11 natural vegetation classes, three developed and mosaiced land classes, and three non-vegetated land classes (Table 2). We derived a distance grid to each land cover class and used these distances for model development rather than land cover class categories themselves. We did this to eliminate any potential sampling bias that could have been caused by researchers’ avoidance of dense land cover during fieldwork. This kind of bias in habitat models results in erroneous species' preference for open habitats. This approach also controls placement errors caused by the cell size resolution of land cover layers.

Table 2 MODIS land cover types based on the International Geosphere Biosphere Programme classification scheme (EOS Data Gateway, 2010).

As a surrogate for human disturbance and development we derived least-cost distances from polygons of urban areas. Urban areas were extracted from the MODIS Land cover dataset (Class 13, Table 2). The computation of least-cost distances was based on the cost distance algorithm implemented in the ArcGIS module Spatial Analyst. This algorithm considers a friction or cost grid that is a raster map where each cell indicates the relative difficulty (or cost) of moving through that cell. A least-cost path minimizes the sum of frictions of all cells along the path, and this sum is the least-cost distance (Adriaensen et al., 2003). In the calculation of cost distances we incorporated information about the terrain to provide more realistic distances regarding the movement of humans. We used slope derived from the SRTM 1-km grid of elevation, thus obtaining terrain-adjusted distance from urban areas. The terrain-adjusted distances between an urban area and a certain point account for not only straight-line distances but also the additional effort humans have to make to move through rugged terrain to reach a point.

Model development and validation

The best-fit habitat model was selected using several modelling techniques: logistic regression (Hosmer & Lemeshow, 1989; Menard, 2002), classification and regression tree (CART; Breiman et al., 1984), Mahalanobis distance (Clark et al., 1993) and maximum entropy (Maxent; Phillips et al., 2006; Phillips & Dudik, 2008). Mahalanobis distance and Maxent do not require use of absence locations for model development. Maxent is an effective tool for estimating a large number of parameters with a small sample size. It eliminates problems associated with data endogeneity and collinearity (Golan et al., 1997). The software used was Maxent v. 3.2.19 (Maxent, 2010). Although Maxent is relatively new in habitat modelling it shows better accuracy and more robustness than some well-established methods (Elith et al., 2006; Phillips et al., 2006; Hernandez et al., 2008).

The Maxent model uses the equation

where P = probability of occurrence of the study object, in our case the presence of Caucasian grouse, and

where λ = a coefficient, the calculation of which is the primary task of the Maxent method, V(x) = habitat variable at an individual cell x of the study area (variables are scaled so that their values lie between 0 and 1 in the training data, using minimum and maximum values), C = linear predictor normalizer chosen so that the exponent is always negative (for numerical stability), Z = density normalizer calculated over the background (Z ensures that q(x) sum to one over the study area), and entropy = -Σq(x)*ln[q(x)].

Using all the habitat variables measured at the 350 training presence locations, we estimated a Mahalanobis Distance model and a Maxent model (hereafter the Maxent single-model approach). We also developed a stepwise Maxent model for which we first developed four Maxent models for grouse distribution: a climate model based on the climate dataset (Table 1), a disturbance model based on human disturbance (i.e. terrain-adjusted distance from urban areas), a land cover model based on distances to land cover classes, and a terrain model based on slope derived from the SRTM dataset. We then multiplied these four models to obtain the final Maxent model (hereafter the Maxent four-model approach).

Our analyses were based on 350 presence/350 absence training locations and 229 presence/229 absence test locations (Fig. 1). Presence test locations in Turkey were adopted from Gottschalk et al. (2007) and Isfendiyaroğlu et al. (2007). Our Armenian and Azeri colleagues provided presence test locations in Armenia and Azerbaijan (see Acknowledgements). We randomly generated the absence locations, as follows: (1) test locations in mountain ranges with no records of the species’ historical presence and (2) training locations within some areas where our repeated surveys over the last 10 years revealed no signs of grouse. On some occasions we surveyed these areas using dogs. The predictive accuracies of the models were defined using the area under the receiver operating characteristic (ROC) curve (AUC; Hanley & McNeil, 1982; Zweig & Campbell, 1993). The ROC curve analysis was performed using SPSS v. 11 (SPSS, Chicago, USA). We applied the final best-fit model to the entire Caucasus and some parts of the Middle East to generate a predictive map of grouse distribution.

The final best-fit model was used to estimate the species’ global population. We did so by weighting the species’ known maximum densities of 3.0–4.8 birds km^-2 (Sikharulidze, 1974; Potapov, 1985) and 2.3–3.9 birds km^-2 (Drovetski & Rohwer, 2000; Baskaya, 2003) with the probabilities of the grouse’s occurrence over the study area:

where N = total population size, d = grouse density (km^-2), a = cell area (1 km² = 1,000 × 1,000 m), P_i = a value of the best-fit model at an i-th cell of the study area, n = number of all cells in the study area, and b = a constant chosen so that P varies between 0 and 1 (being a product of four models whose values range between 0 and 1, the final model is condensed between 0 and 1).

Results

Caucasian grouse were found on north-facing slopes in open habitat (Table 3) at 1,191–3,011 m (mean = 2,494 ± SD 255 m). Student’s t-tests showed no statistically significant differences in habitat variables between training presence locations from radio-tracked birds and the additional training presence locations from birds that were not radio-tracked. Models other than Maxent had considerable omission and commission errors. The Mahalanobis distance model under-predicted the species’ range; i.e. it omitted many of the species' presence locations not used for training. In contrast, logistic regression and classification tree analyses predicted the species’ presence in many areas where Caucasian grouse have never to our knowledge been recorded (e.g. in Iran’s Elburz, or Alborz, Mountains, Turkey’s midland and areas east of the Azov Sea, Russian Federation). The Maxent four-model approach performed best (Table 4). Visual examination of probability maps derived by applying these models to the entire study area showed that only the Maxent four-model approach was strongly consistent with the species' known distribution. Therefore we present the final best-fit model only (Fig. 2). Appendices 1–4 summarize the models that the Maxent four-model approach was derived from.

Fig. 2 Probability of Caucasian grouse occurrence throughout the Caucasus, derived from Maximum Entropy modelling (see text for details).

Table 3 Distribution of Caucasian grouse at 350 training locations by aspect (derived from 1:50,000 Soviet military topographic maps) and land cover classes (Table 2).

Table 4 Measure of predictive accuracy of models for habitat use by Caucasian grouse using the AUC and a test dataset (n of presence and absence = 229).

The Maxent procedure estimates relative contributions of variables to the model by adding the increase in regularized gain to the contribution of the corresponding variable in each iteration of the training algorithm, or by subtracting from it if the change in absolute value of λ is negative. For the climate model annual mean temperature (Bio1) had 67% of the total contribution and each of the other bioclimatic variables had < 6%. The results of the Jackknife test of variable importance performed on presence/absence test locations suggested that the habitat variable with the highest AUC, when used alone, was mean temperature of the warmest quarter (Bio10). The habitat variable that decreased the AUC the most when it was omitted was precipitation seasonality (Bio15).

The analysis of the land cover model suggested that distance from deciduous broad-leaf forest (Distance4) and from grasslands (Distance10) had the largest contributions, 30 and 25% respectively, to the model. The results of the jackknife test of variable importance performed on test locations suggested that the habitat variable with the highest AUC when used alone was Distance4. The habitat variable that decreased the AUC the most when it was omitted was Distance10.

Maxent modelling creates response curves for each variable. These plots reflect the dependence of predicted suitability both on the selected variable and on dependencies induced by correlations between the selected variable and other variables (i.e. marginal response curves). The analysis of these showed that the probability of the grouse occurrence reached its maximum at values of Bio1 of c. 0^o, Bio10 c. 8^o and Bio15 c. 55 mm (Fig. 3). The probability had a bell-shaped response to terrain-adjusted distance from urban and built-up areas as well as ruggedness, reaching its maximum at a cost distance of c. 700 and slope c. 18^o. The probability of grouse presence was highest in grassland near the edge of deciduous broad-leaf forest (see the curves of Distance4 and Distance10 in Fig. 3). Based on the Maxent four-model approach and known species' densities, the grouse population was estimated to be 30,203–63,034, which is less than the 68,000–84,300 estimated by BirdLife International (2008).

Fig. 3 Probability of the occurrence of Caucasian grouse in relation to the most important variables identified by Maximum Entropy modelling (Bio1, annual mean temperature; Bio10, temperature of warmest quarter; Bio15, precipitation seasonality; Cost distance, terrain-adjusted distance from urban and built-up areas; Distance4, distance from deciduous broad-leaf forest; Distance10, distance from grasslands).

Discussion

In the Caucasus there are ongoing efforts to protect and restore wildlife communities and thus the areas where they occur. Because resources available for conservation efforts are limited there is a need to decide which areas require priority attention. One of the best approaches is to identify such areas by outlining and overlaying the distributions of threatened and restricted-range species (Langhammer et al., 2007). Spatial quantitative assessment of habitat requirements or habitat modelling is an effective way to outline species' distributions. Our study species, the Caucasian grouse, is both threatened and has a restricted range, and our study quantified habitat suitability for the species and used the habitat model to estimate population size.

Our objective was to provide an accurate method to quantify habitat requirements for the Caucasian grouse and estimate its total population. The use of Maximum Entropy modelling considerably improved model accuracy and the inclusion of radio-telemetry for data collection enabled us to obtain a large sample of observations in habitats where otherwise it would have been difficult to detect this secretive species. Radio-telemetry also helped us acquire full temporal and spatial coverage of the habitats where the species occurs. In comparison with previous studies (e.g. Gottschalk et al., 2007) our model is applicable to Caucasian grouse on a broader scale because the sample for model development was obtained from almost the entire range of the species and across a broader spectrum of landscape types. The utility of our final model is demonstrated by its identification of small patches of habitat suitable for the Caucasian grouse in Iran, just south of the southernmost part of Armenia. These patches are close to the only areas in Iran where the Caucasian grouse has been recorded (Gharamani, 1971; Scott, 1976). The Iranian population is the southernmost point of the species’ known range.

Ecological studies of the Caucasian grouse have generally shown that the species is associated with rugged areas at 1,300–3,000 m altitude on northern slopes covered in Rhododendron and birch thicket, and has negative correlations with proximity to roads and human activity (Averin, 1938; Kutubidze, 1961; Tkachenko, 1966; Sikharulidze, 1974; Johnsgard, 1983; Potapov, 1985; Atkinson et al., 1995; Baskaya, 2003; Gokhelashvili et al., 2003; Klaus & Storch, 2003; Gottschalk et al., 2007; Isfendiyaroğlu et al., 2007). Even though the species is more often encountered on northern slopes, it uses south-facing grassy slopes, especially south-west-facing slopes, as lek sites during the breeding season (Gokhelashvili et al., 2003). The species’ aspect preference is possibly explained by the fact that livestock grazing has less impact on the species’ cover (tall grass, Rhododendron thicket, shrubland and forest) and food on north-facing slopes where vegetation re-growth and recovery rates are usually greater than on south-facing slopes because of the higher moisture index. Furthermore, Caucasian grouse make burrows in snow cover, where they pass nights and wait out bad weather. When predators approach the burrow entrance the grouse breaks through the snow layer over the burrow end to escape. This behaviour is possible only in fluffy snow cover. On southern slopes snow melts during the day and freezes at night whereas it remains fluffy on north-facing slopes.

In our analysis we did not use altitude as a surrogate for climate because we wanted the models to be applicable to a wide geographical range, over which climate will be determined by additional factors (e.g. latitude, terrain, air currents, coastal streams). Analysis of contributions of climatic variables to the Climatic model suggests that annual mean temperature (Bio1), mean temperature of the warmest quarter (Bio10) and precipitation seasonality (Bio15) are the variables that have the greatest influence on the distribution of the species. To see how these results relate to the findings of previous ecological studies, we used ArcGIS to identify the range between two minimum point values on each axis of these variables (Bio1, Bio10, Bio15; Fig. 3). The range of Bio1 describes the elevation range where the species generally occurs. Overlaying the ranges of Bio10 and Bio15 shows the potential distribution of land cover types suitable for the species within this elevation range. Some areas where suitable land cover types were absent within this suitable climatic range were generally areas modified by human activities.

Even though our encounters with Caucasian grouse were more frequent in open than in dense habitat (Table 3), our land cover model suggested that the species prefers grassland in proximity to the edges of deciduous broad-leaf forest. Gottschalk et al. (2007) also indicated the importance of proximity to forest in Turkey but did not specify which forest type is most important.

The species' occurrence had a bell-shaped response to terrain-adjusted distance from urban and built-up areas instead of the expected linear or near linear relationship. The reason for this bell-shaped relationship is that areas at a cost distance of > 700 lie above 3,000 m, where the climate is harsh and the land cover unsuitable for Caucasian grouse. The reason for the species' bell-shaped response to slope is simply that we did not have locations in areas with slope > 18^o.

Our estimate of the entire Caucasian grouse population is less than previous estimates (BirdLife International, 2008), in the calculation of which we participated. We believe the estimate presented here is more accurate because previous estimates were based on either less accurate regional models or expert judgements. Actual population size could be lower than our estimate because our results are based on a relatively coarse 1,000-m spatial resolution.

Our model provides a tool for identification of potential Caucasian grouse populations, and habitat for the species, throughout its range. The model predicts the probable distribution of grouse and, based on knowledge of the species’ density, facilitates estimation of population size. Searching for the grouse and estimating numbers is difficult over vast, rugged and climatically inclement areas; our model will facilitate improved detection and more accurate population estimates. Another practical use of the model will be for prediction of connectivity between populations and thus planning of corridors for conservation purposes. The research also provides an algorithm that could be used for modelling habitats for similar species. Finally, as threatened and endemic species are often used for setting site-based conservation priorities, this habitat model could help to identify new Key Biodiversity Areas in the Caucasus. The Ministry of Environmental Protection and Natural Resources of Georgia is going to use the results of this study to reshape existing protected areas and identify new ones.