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 Reference Guisan and Zimmermann2000, MacNally Reference MacNally2000, Lehmann et al. Reference Lehmann, Overton and Leathwick2002, Scott et al. Reference Scott, Heglund, Haufler, Morrison, Raphael, Wall and Samson2002). Mackenzie et al. (Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002, Reference Mackenzie, Nichols, Hines, Knutson and Franklin2003) 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. Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002). 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. Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002, Reference Mackenzie, Nichols, Hines, Knutson and Franklin2003). The approach has proved to be particularly useful in the case of rare and cryptic species (Pellet and Schmidt Reference Pellet and Schmidt2005, Schmidt and Pellet Reference Schmidt and Pellet2005, Buij et al. Reference Buij, McShea, Campbell, Lee, Dallmeier, Guimondou, Mackage, Guisseougou, Mboumba, Hines, Nichols and Alonso2007, Kroll et al. Reference Kroll, Risenhoover, McBride, Beach, Kernohan, Light and Bach2008, Srinivas et al. Reference Srinivas, Venugopal and Ram2008, Watson et al. Reference Watson, Weckerly, Hatfield, Farquhar and Williamson2008).
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 Reference Collar, Coles and Collar1980, Johnsgard Reference Johnsgard1991, Goriup Reference Goriup1997). 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 Reference Collar, Coles and Collar1980, Schulz and Seddon Reference Schulz, Seddon, Saint Jalm and Van Heezik1995, Goriup Reference Goriup1997, Combreau et al. Reference Combreau, Launay and Lawrence2001, Birdlife International 2004, Le Cuziat et al. Reference Le Cuziat, Lacroix, Roche, Vidal, Médail, Orhant and Béranger2005, Tourenq et al. Reference Tourenq, Combreau, Lawrence, Pole, Spalton, Xinji and Launay2005). 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 Reference Etchécopar and Hüe1964, Anonymous Reference Coles and Collar1980, Goriup Reference Goriup1997, Chammem et al. Reference Chammem, Khorchani, Boukhris, Combreau, Chniti and Hammadi2003, Isenmann et al. Reference Isenmann, Gaultier, El Hili, Azafzaf, Dlensi and Smart2005). Nowadays, only small relict populations still survive in the most remote Saharan areas in the south (Chammem et al. Reference Chammem, Khorchani, Boukhris, Combreau, Chniti and Hammadi2003, Isenmann et al. Reference Isenmann, Gaultier, El Hili, Azafzaf, Dlensi and Smart2005).
Recently, considerable attention has been given to the conservation of Houbara in Tunisia (Chammem et al. Reference Chammem, Khorchani, Boukhris, Combreau, Chniti and Hammadi2003, Azafzaf et al. Reference Azafzaf, Sande, Evans, Smart and Collar2005). 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 Reference Cramp and Simmons1983, Launay et al. Reference Launay, Roshier, Loughland and Aspinall1997, Osborne et al. Reference Osborne, Launay and Gliddon1997, Chammem et al. Reference Chammem, Khorchani, Boukhris, Combreau, Chniti and Hammadi2003) and because of its high sensitivity to human disturbance (Lavee Reference Lavee1985, Reference Lavee1988, Seddon et al. Reference Seddon, Saint Jalme, Van Heezik, Paillat, Gaucher and Combreau1995, Osborne et al. Reference Osborne, Launay and Gliddon1997, Van Heezik and Seddon Reference Van Heezik and Seddon1999, Le Cuziat et al. Reference Le Cuziat, Lacroix, Roche, Vidal, Médail, Orhant and Béranger2005) 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 Reference Floret and Pontanier1982). The climate is arid and hot, with annual rainfall < 200 mm and average annual temperature exceeding 21°C (Floret and Pontanier Reference Floret and Pontanier1982). 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. Reference Akrimi, Ben Hammouda, Abaab, Neffati and Sghaier2000). 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).
Map of the study area.
Data collection
The study area was divided into a grid consisting of 32 cells of 400 km2, 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. Reference Launay, Roshier, Loughland and Aspinall1997, Yang et al. Reference Yang, Qiao, Combreau, Gao and Zhong2003). 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 Reference Le Houérou and Le Floc’h2001). 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. (Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002, Reference Mackenzie, Nichols, Hines, Knutson and Franklin2003) 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. Reference Hingrat, Saint Jalme, Ysnel, Lacroix, Seabury and Rautureau2004, Reference Hingrat, Saint Jalme, Ysnel, Le Nuz and Lacroix2007), 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. Reference Bailey, Simons and Pollock2004, MacKenzie Reference MacKenzie2006). 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. (Reference Betts, Rodenhouse, Scott Sillett, Doran and Holmes2008) and integrated in the list of site covariates an autocorrelation term (AUTO) that we calculated according to Augustin et al. (Reference Augustin, Mugglestone and Buckland1996) and Klute et al. (Reference Klute, Lovallo, Tzilkowski, Scott, Heglund, Morrison, Haufler, Raphael, Wall and Sampson2002). For a square i, AUTOi = ΣWijYj/ΣWij, where Wij is the weight given to square j, calculated as the inverse of the geographical distance between grid squares i and j, and Yj 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 Reference Burnham and Anderson2002). We also considered model weight as an indication of the relative support for each model (Burnham and Anderson Reference Burnham and Anderson2002, Wintle et al. Reference Wintle, McCarthy, Volinsky and Kavanagh2003). 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 Reference Burnham and Anderson2002, Wintle et al. Reference Wintle, McCarthy, Volinsky and Kavanagh2003). 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.
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.
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.
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. (Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002, Reference Mackenzie, Nichols, Hines, Knutson and Franklin2003) 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 Reference Dominguez-Casanova1989, Martín et al. Reference Martín, Nogales, Hernández, Lorenzo, Medina and Rando1996, Medina Reference Medina1999, Le Cuziat, et al. Reference Le Cuziat, Lacroix, Roche, Vidal, Médail, Orhant and Béranger2005) for Chlamydotis undulata, and also in China for Chlamydotis macqueenii (Combreau et al. Reference Combreau, Qiao, Lawrence, Gao, Yang and Launay2002, Tourenq et al. Reference Tourenq, Combreau, Lawrence, Pole, Spalton, Xinji and Launay2005). 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 Reference Abaab1986, Nasr Reference Nasr1993, Bourbouze Reference Bourbouze2000, Nasr et al. Reference Nasr, Ben Salem, Rachidi and Medouni2000). 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. Reference Floret, Le Floc’h and Pontanier1987, Nasr Reference Nasr1993, Nasr et al. Reference Nasr, Ben Salem, Rachidi and Medouni2000). 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 Reference Sghaier2004). 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 Reference Le Houérou, Di Castri, Goodall and Specht1981, Abaab et al. Reference Abaab, Bedrani, Bourbouze and Chiche2000, Bourbouze Reference Bourbouze2000). 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.
Acknowledgements
We thank M. Abid, H. Tir, S. Sekrafi, F. Abdelli, A. Ouled Belgacem and A. Barmat for helping with data collection. We are also grateful to J. Hines for his help and for instructions regarding the use of program PRESENCE for data analyses and interpretation of results.
