Using a capture-recapture approach for modelling the detectability and distribution of Houbara Bustard in southern Tunisia

Summary

Modelling the distribution of species of conservation concern is an important issue in population ecology. Classically, logistic regression analyses are conducted to estimate species’ distributions from detection/non-detection data in a sample of sites and to test for the significance of several environmental variables in predicting the probability of occurrence. These modelling approaches assume that species detection probability is constant and equals one in all sampled sites, which is critical, notably in the case of rare, shy and cryptic species. The capture-recapture-like approach developed by Mackenzie et al. (2002, 2003) provides a reliable tool that accounts for imperfect detection when estimating species occurrence, as well as for assessing the relevance of site features as predictors of species occurrence probability. The aim of this study was to explore the possibility of using this approach in the context of Houbara Bustard Chlamydotis undulata in southern Tunisia. Our results show once more the low detectability of this emblematic species and stress the need to take this factor into account when estimating Houbara spatial distribution. The distribution of Houbara in southern Tunisia is more likely to be shaped by human-related than by habitat factors. In particular, Houbara occurrence was positively associated with site remoteness and camel numbers. Houbara seemed to avoid areas with high human presence and shared the most remote and agriculture-free zones with free-ranging camels.

Introduction

Statistical modelling and analyses of species’ distributions using biological survey data and biotic and abiotic predictors are widely used in conservation planning (Guisan and Zimmermann 2000, MacNally 2000, Lehmann et al. 2002, Scott et al. 2002). Mackenzie et al. (2002, 2003) proposed use of a capture-recapture-like approach that can provide an appropriate solution for species that are difficult to detect and when many populations are inevitably missed. This approach controls for heterogeneity in species detectability among sampling locations when estimating species occupancy from count replicates. It also permits an estimate of the effects of various environmental factors on both detectability and occurrence probability of a species by incorporating them as covariates in the estimation (Mackenzie et al. 2002). This provides a less biased and/or more precise estimate of species occupancy than the naive estimate provided by the simple proportion of sites where the species was actually observed (Mackenzie et al. 2002, 2003). The approach has proved to be particularly useful in the case of rare and cryptic species (Pellet and Schmidt 2005, Schmidt and Pellet 2005, Buij et al. 2007, Kroll et al. 2008, Srinivas et al. 2008, Watson et al. 2008).

The Houbara Bustard Chlamydotis undulata is one of the most threatened birds inhabiting desert areas across North Africa, the Middle East and western Asia (Collar 1980, Johnsgard 1991, Goriup 1997). Over the past few decades, Houbara populations have suffered a serious decline, and the species is now classified as ‘Vulnerable’ (IUCN 2010). Severe habitat loss and/or degradation in all parts of its geographic range, due to increasing human pressures on natural habitats, as well as overhunting and poaching, are believed to be the main causes of Houbara decline (Collar 1980, Schulz and Seddon 1995, Goriup 1997, Combreau et al. 2001, Birdlife International 2004, Le Cuziat et al. 2005, Tourenq et al. 2005). In Tunisia, the North African Houbara Bustard Chlamydotis undulata undulata was largely distributed in desert areas of the central and southern parts of the country, where it inhabited steppes with sparse shrubs and grasses (Etchécopar and Hüe 1964, Anonymous 1980, Goriup 1997, Chammem et al. 2003, Isenmann et al. 2005). Nowadays, only small relict populations still survive in the most remote Saharan areas in the south (Chammem et al. 2003, Isenmann et al. 2005).

Recently, considerable attention has been given to the conservation of Houbara in Tunisia (Chammem et al. 2003, Azafzaf et al. 2005). Houbara Bustard is now protected by Tunisian law, and some protected areas have been created within its area of distribution in southern Tunisia. However, the status of this species in the wild is poorly known, and there is a lack of detailed information on the ecology of remaining wild populations in southern Tunisia. For instance, the spatial distribution of Houbara in southern Tunisia and the ecological factors affecting Houbara occurrence probability in a given site have never been investigated appropriately. Given their rarity and extreme vigilance, Houbara are difficult to detect, so their occurrence at a given site can easily be underestimated, resulting in an incomplete spatial distribution. This could lead to biased estimates and erroneous conclusions about the effects of environmental factors on Houbara occurrence probability.

In this work, we used survey data on Houbara in El-Ouara region in south-eastern Tunisia to estimate the occupancy of the species and to assess the relevance of several habitat and human variables as predictors of Houbara occurrence probability. Program PRESENCE was used (1) to estimate the detection probability of Houbara, (2) to account for the presumed low detection probability in estimating its occurrence, and (3) to test the significance of the effects of habitat features and human pressures on the occurrence of Houbara in this area. Our main goal was to assess the relative importance of habitat features versus human pressures in shaping Houbara occurrence in the study area. Given that Houbara is a habitat generalist species (Cramp and Simmons 1983, Launay et al. 1997, Osborne et al. 1997, Chammem et al. 2003) and because of its high sensitivity to human disturbance (Lavee 1985, 1988, Seddon et al. 1995, Osborne et al. 1997, Van Heezik and Seddon 1999, Le Cuziat et al. 2005) we expected Houbara occurrence more likely to be affected by the intensity of human presence and land use rather than by habitat variables.

Methods

Study area

Our work was carried out in El-Ouara area, along the Libyan border in south-eastern Tunisia, between 32.00° and 33.20°N and from 10.26° to 11.58°E (Fig. 1). This area is characterised by flat to gently undulating topography with a sand and gravel substrate and containing flat saline depressions (sebkha and chott) (Floret and Pontanier 1982). The climate is arid and hot, with annual rainfall < 200 mm and average annual temperature exceeding 21°C (Floret and Pontanier 1982). The vegetation generally consists of patches of short perennial grasses and scattered dwarf shrubs (chamaephytes), such as Stipagrostis pungens, Anthyllis sericea, Gymnocarpos decander, Hammada schmittiana, Hammada scoparia, Traganum nudatum, Limoniastrum guyonianum Retama raetam and Ziziphus lotus (Akrimi et al. 2000). Urban and agricultural developments border the study area, and are mainly concentrated around the road joining Ben Guerdane in the north to Tataouine in the west and Dhehibat in the south (Figure 1). This area is increasingly used as a permanent grazing zone for camels, as well as for sheep and goats (hereafter called small ruminants).

Figure 1. Map of the study area.

Data collection

The study area was divided into a grid consisting of 32 cells of 400 km², each measuring 20 km x 20 km. During March and April 2004, these 32 sites were surveyed for Houbara. Within each site, five sampling stations 3 km apart were randomly selected along a line transect. Each sampling station was visited once for a single 30 min survey occasion, early in the morning or at the end of the afternoon, when Houbara are most active and probability of detection is highest. For a given site, all sampling occasions were carried out on the same day. The time of each sampling occasion was recorded (HOUR) and used as a sampling covariate as it can be assumed to affect Houbara detection probability. All sampling was conducted by the same four-observer team (M. Chammem, M. Abid, H. Tir and S. Sekrafi) and only under good meteorological conditions, notably when wind velocity was light (< 3 knots). During each sampling occasion, the observers walked in different directions and searched for Houbara at unlimited distances using binoculars. Because it is very difficult to observe this vigilant, shy and cryptic species directly, observers concentrated on searching for indices of Houbara presence, notably tracks, that are easily recognisable (Launay et al. 1997, Yang et al. 2003). From their long field experience with the species and the study area, the observers were familiar with Houbara footprints and where to search for them, in particular by orienting search effort towards the vicinity of plants appreciated by Houbara. For each sampling occasion we retained a value of 1 if Houbara were detected and a value of 0 if not detected. This provided us with five detection/non-detection (1/0) records from each of the 32 grid cells that we used to estimate Houbara occupancy. However, the data obtained are not appropriate for estimating species abundance.

During each sampling occasion, the observers also recorded the number of small ruminants (SRUM) and the number of camels (SCAM) in the surrounding area, as they were assumed to compete with Houbara, and consequently affect its occurrence probability. For each grid cell, five values of SRUM and SCAM were obtained and were used to calculate the mean value of each variable for the site. These mean values were then used as descriptors of human pressure, as well as two further variables, derived from official maps. We used 1:100,000 scale maps of Medenine and Tataouine governorates (Cartes Agricoles 2002) to measure the intensity of urbanisation for each grid cell, calculated as the percentage cover of urbanized area (URBA), and the intensity of agricultural activity (percentage cover of agricultural land; AGRI). We also used these maps to collect data on the percentage cover of each of the main geomorphological units identified in the region: plains (PLAI), sandy areas (SAND) and salty areas (SALT). These data were complemented by those concerning vegetation cover that we extracted from a 1:1,000,000 scale vegetation map (Le Houérou and Le Floc’h 2001). We used this map to measure the proportion of each site covered by the main vegetation types in the region: shrubs (SHRU), desert woody plants (WOOD) and halophytes (HALO). Thus we obtained four descriptors of human pressure (SCAM, SRUM, URBA and AGRI) and six habitat variables (PLAI, SAND, SALT, SHRU, WOOD and HALO) in each of the 32 grid cells.

Data analysis

Given that the original human and habitat variables were correlated, we first conducted two separate principal component analyses (PCA) to summarise each subset of variables (human and habitat) into a smaller number of independent factors. The factors extracted from the PCAs were then used as site covariates when estimating Houbara occupancy.

We used the approach of Mackenzie et al. (2002, 2003) as implemented in program PRESENCE (available for download from http://www.proteus.co.nz) to estimate the proportion of sites occupied by Houbara from the recorded lists of detection/non-detection over the five spatial replicates within each cell. This approach takes into account the possible effects of different parameters on both occupancy (Ψ) and detection (p) probabilities by incorporating them as covariates in the estimations. Given that our survey work was conducted during the same breeding season (March–April), and as Houbara are known to be sedentary with relatively limited home-ranges during this period (Hingrat et al. 2004, 2007), we assumed that Houbara occupancy was invariable during the period of study. This satisfied the assumption of “closed distribution” required when applying this approach.

We employed a two-step process to estimate site occupancy parameters. We first selected the best model for detectability by completing a model assuming detectability was constant (p(.)) with a second model assuming that detectability varies with HOUR (p(HOUR)), while holding occupancy constant (Ψ(.)) for both models. We then used the best detectability model to estimate occupancy by defining a set of a priori models we believed might explain site occupancy of Houbara in our study area (Bailey et al. 2004, MacKenzie 2006). The simplest model assumed that occupancy was constant throughout the study area (Ψ(.)). Then we developed candidate models accounting for the factors extracted from the PCAs of the original human and habitat variables as site covariates. Because our observation units are spatial locations (32 grid squares), a possible problem of spatial autocorrelation may have occurred in our data. In order to take into account any possible spatial covariance, we followed Betts et al. (2008) and integrated in the list of site covariates an autocorrelation term (AUTO) that we calculated according to Augustin et al. (1996) and Klute et al. (2002). For a square i, AUTO_i = ΣW_ijY_j/ΣW_ij, where W_ij is the weight given to square j, calculated as the inverse of the geographical distance between grid squares i and j, and Y_j is the response in square j(0 if the species is absent and 1 if it is present). For a given square, AUTO may take a value near 1 if the Houbara occurs in all neighbouring squares and a value near 0 if it occurs only in the more distant squares. However, given that AUTO was calculated without taking into account the problem of detectability, it should be viewed as a simple index of Houbara records in neighbouring areas. Because of the relatively low number of sampled sites (32) and given that the number of sites with multiple detections was small (see Results), our data were insufficient to build models with many covariates (J. Hines pers. comm.). For this reason, we only employed models with zero or one site covariate, giving us a total of six candidate models.

The competing models were ranked according to their AICc (Akaike’s Information Criterion for small sample size). The most parsimonious model was that with the smallest AICc, but if the absolute difference between two models was < 2, they were not regarded as different (Burnham and Anderson 2002). We also considered model weight as an indication of the relative support for each model (Burnham and Anderson 2002, Wintle et al. 2003). Because no covariate was entered in more than one model (see above), there was no need to use a model-averaging technique for the estimation of the effects of covariates on occupancy (Burnham and Anderson 2002, Wintle et al. 2003). For this reason, we considered the results of the best-ranked model in investigating the significance and direction of the relationship between the corresponding covariate and Houbara occupancy. However, we stress that our models were used to investigate factors that we thought would influence Houbara occupancy, and these do not necessarily represent the complete range of factors that could influence Houbara occupancy in our study area.

Results

The first PCA summarised the four measured human variables into two independent factors accounting for 67.50% of the variance in the original human dataset. The first factor (HUM1) represents an axis of decreasing remoteness, as it was positively correlated with URBA (r = 0.671, P < 0.0001) and AGRI (r = 0.761, P < 0.0001), and decreasing camel abundance, as it was negatively correlated with SCAM (r = -0.772, P < 0.0001). The second factor (HUM2) represents an axis of increasing small ruminant abundance, as it was positively correlated with SRUM (r = 0.945, P < 0.0001). The second PCA summarised the original habitat variables into two independent factors, accounting for 62.26% of the original variance. The first factor (HAB1) is positively correlated with SALT (r = 0.743, P < 0.0001) and HALO (r = 0.922, P < 0.0001). It provides a measure of the importance of salty areas covered by halophytes. The second factor (HAB2) is positively correlated with SAND (r = 0.763, P < 0.0001) and WOOD (r = 0.489, P = 0.0046) but negatively correlated with PLAI (r = -0.736, P < 0.0001). It represents an axis of increasing importance of sandy areas covered by desert woody shrubs. The descriptors of human and habitat variables derived from the PCAs of the original variables were not correlated (HUM1 and HAB1: r = -0.172, P = 0.3470; HUM1 and HAB2: r = 0.250, P = 0.1671; HUM2 and HAB1: r = -0.224, P = 0.2169; HUM2 and HAB2: r = -0.119, P = 0.5143). These descriptors were thus entered as site covariates in the occupancy models to check for their relative importance in predicting Houbara distribution in our study area.

Houbara were recorded on only 13 of the 160 sampling occasions (32 grid squares x five replicates). These 13 detections were recorded in nine different grid squares, providing a naive occupancy estimate of 0.28 (9/32). All these records corresponded to fresh tracks, and no direct observation was made, which shows the extreme shyness of the species and at the same time highlights the need to account for detectability issues in the estimation of Houbara occupancy. We first found that the constant detectability model, i.e. the model assuming that detection probability was constant across sampled sites (p(.)), was better ranked, and hence more supported by the data, than the model assuming that detection probability varied with time of day ((p(HOUR)) (Table 1). Both models showed that Houbara detection probability was low (Table 1). The constant detectability model was then used to estimate and model Houbara occupancy. Among the six occupancy models we tried, no model received exclusive support from the data, as no model had an AIC weight > 0.90 (Table 2). Nonetheless, there was a tendency for the model including HUM1 as a site covariate to be better supported by the data than the other models. Indeed, the results of AIC selection showed that the HUM1 model was the best-ranked compared to those including HUM2, HAB1, HAB3 or AUTO as site covariates (Table 2). This model received the highest support from the data (weight = 55%). However, this does not necessarily mean that this model is the most statistically significant. For instance, there was no substantial difference between the HUM1 model and the model assuming occupancy constant, as the ∆AICc was slightly less than 2 (Table 2). All models provide occupancy estimates varying between 0.41 and 0.44, higher than the naive estimate (0.28), which assumes that detection probability was equal to 1, suggesting that it underestimated Houbara in the study area.

Table 1. Results of the comparison between the constant-detectability model and the model assuming that detectability varies with time of day (HOUR). AICc = Akaike’s Information Criterion corrected for small sample size. P = detection probability.

Table 2. Summary of model selection results regarding the occupancy of Houbara Bustard in the study area. AICc = Akaike’s Information Criterion corrected for small sample size. p = estimated detection probability. Ψ = estimated occurrence probability. HUM1 and HUM2 = first and second factors derived from the PCA of the original human variables. HAB1 and HAB2 = first and second factors derived from the PCA of the original habitat variables. AUTO = index of presence spatial autocorrelation.

Table 3. Beta estimates, standard errors, and 90% confidence intervals of the covariates used in modeling Houbara occupancy.

Discussion

Our work aimed to examine the potential of the approach of Mackenzie et al. (2002, 2003) in estimating Houbara occupancy in one of the most remote areas in Tunisia where this species still survives in the wild, and at the same time to assess the relative importance of habitat factors versus human-induced disturbance in shaping Houbara occurrence probability. The results show once more the low detectability of this species, underlining the need to take detectability issues into account when estimating the spatial distribution of Houbara and when investigating the association between Houbara occurrence and environmental factors. Our results suggest that the proportion of all the sampled sites where Houbara were actually recorded underestimated Houbara occupancy in the study area. This naïve estimate assumes detection probability to be equal to one, which was not supported by our data. Furthermore, our results suggest that the distribution of Houbara in the study area is more likely to be shaped by human-related, rather than habitat factors. In particular, Houbara occurrence was negatively associated with the first factor derived from the PCA of the original human variables (HUM1), which represents an axis of increasing urban and agricultural land cover and decreasing camel numbers. This would mean that Houbara occurrence is negatively associated with land transformation into urban and agricultural areas but positively associated with camel numbers.

Our results suggest that habitat factors have no significant effect on the occurrence of Houbara in the study area as they seemed to use the different habitat types available in this area in similar ways. Even if some habitat preferences exist, they seem to have a negligible effect compared to human-related factors. However, the non-significant relationship between habitat and Houbara occurrence that we found cannot be generalised to other areas within the species’s geographic range. For instance, this lack of association could simply be due to the fact that the relationship was investigated at an inappropriate scale. Habitat parameters do not vary much across the sampled sites because of the relatively small study area and this low level of habitat heterogeneity may explain why we could not detect potential associations between habitat features and Houbara occurrence. Further large-scale surveys across its geographic range are likely to provide better information about the role of habitat features in shaping the large-scale spatial distribution of this species.

The negative relationship we found between Houbara occurrence probability and HUM1 suggests that Houbara avoid areas with relatively high human presence and agricultural land and inhabit more remote areas. This result is consistent with several previous studies on desert bustards, such as in the Canary Islands and Morocco (Dominguez-Casanova 1989, Martín et al. 1996, Medina 1999, Le Cuziat, et al. 2005) for Chlamydotis undulata, and also in China for Chlamydotis macqueenii (Combreau et al. 2002, Tourenq et al. 2005). In our study area, the last decade has been accompanied by an increasing intensity of land transformation for agriculture and urbanisation that is likely to result from changes in human lifestyles from a non-sedentary to a sedentary system (Abaab 1986, Nasr 1993, Bourbouze 2000, Nasr et al. 2000). The change from a nomadic lifestyle in southern Tunisia follows the recent transition from an old subsistence economy to a new market economy, mainly involving the privatisation of land and the creation of farmland for cereal and olive crops, as well as urbanisation and development of road networks (Floret et al. 1987, Nasr 1993, Nasr et al. 2000). The noise generated by humans and their vehicles and the presence of dogs are likely to disturb Houbara and keep them away from human settlements. Illegal hunting may also constitute a serious threat to the remaining Houbara. Although hunting of Houbara is prohibited by Tunisian law, poaching is carried out by local poachers as well as Arab falconers, who regularly visit southern Tunisia. Knowing the extreme vigilance of Houbara, this growing human presence and land transformation constitutes a serious threat that could rapidly lead to the total disappearance of the remaining populations of the species in southern Tunisia.

Interestingly, our results show a positive relationship between Houbara occurrence and camel numbers. Two hypotheses can be proposed to explain such an association. First, both species are constrained to use the most remote and human-free areas where they share the same areas without any tendency to live together. This is the consequence of the camel breeding system adopted in the region. Camels are not allowed to graze around villages and cultivated areas, but are released into more remote zones where they are free to roam and live with almost no human control for long periods of time (Sghaier 2004). Within these areas, camels do not seem to negatively affect Houbara, as the two species do not compete for food. The alternative explanation could be that Houbara tend to follow camels because they are beneficial. Houbara may be attracted by the abundance of prey associated with camels and their dung (insects and other invertebrates) or they may use camel herds to help avoid predators, such as golden jackal Canis aureus. Sheep and goats were not found to influence Houbara occurrence probability. This could be due to the traditional small ruminant breeding system, based on seasonal movements and nomadism, and which may not strongly affect habitat quality, as it does not lead to overgrazing or disturbance (Le Houérou 1981, Abaab et al. 2000, Bourbouze 2000). The compatibility of Houbara with livestock that we report in this study has important conservation implications. It suggests that future Houbara conservation programmes in southern Tunisia should aim to integrate different land uses and focus on sustainable rangeland management rather than exclude livestock grazing. Further fine-scale investigations of the interactions between Houbara and livestock and their possible repercussions on Houbara abundance, habitat use, feeding success, health and stress are needed to make a more accurate evaluation of the impact of livestock on long-term survival of the studied Houbara population.

Overall, our work suggests that the distribution of Houbara in southern Tunisia is more strongly affected by human factors, especially agricultural development than by natural habitat characteristics. Habitat features do not seem to account for Houbara occurrence, likely because the species is now forced to use the most remote areas which seem to provide sub-optimal habitats. Continuing the current policy of land privatisation and agricultural development in this area, through the creation of water points and the development of agricultural projects without taking into account the ecological requirements and behavioural features of Houbara, is likely to result in the rapid and complete extinction of this species in southern Tunisia. Urgent conservation programmes dealing with these new human pressures are needed if we are to conserve the remaining populations of this threatened species.