Populations at the periphery of a distribution often are genetically and ecologically distinct but tend to be less represented in conservation planning. Peripheral populations may exhibit different behaviors and persist longer than core populations as landscapes change, especially because of human-induced habitat loss and climate change. We examined how landscape factors influence occupancy patterns of Toad Headed Agamas (Phrynocephalus versicolor) along the northern periphery of their range in Asia, a region experiencing increasing development and warming conditions. We collected detection/nondetection data during surveys of 180 sites in Mongolia in 2016. We then developed a set of 70 candidate models that included the single, additive, and interacting effects of nine covariates on occupancy and detection probability and used model selection techniques to determine the best model in the set. We detected agamas at 89 sites (49.4%) and during 141 (39.2%) of 360 surveys. Only one model had strong empirical support, one which included the additive effects of forb cover, grass cover, and ruggedness on occupancy probability and wind speed on detection probability. All four covariates had negative effects, suggesting that ideal conditions for occupancy were areas with little vegetation and topographically flat or gently rolling and that detection was higher in low wind. Average predicted occupancy across all sites was 55%. Our results indicated agamas were more sensitive to vegetation cover in this area than in other parts of their range. Agamas may benefit from future climate conditions that reduce vegetation but face negative impacts from increases in landscape ruggedness because of development activities.
Populations at the periphery of a species distribution often are less represented in conservation planning than are core populations (Channell and Lomolino, 2000). Peripheral populations typically occupy less-favorable habitats, exhibit lower and more-variable densities, and can be more isolated and fragmented than are core populations and, as a result, may be considered less valuable (Brown et al., 1996; Channell and Lomolino, 2000). Conversely, peripheral populations are often genetically and morphologically distinct, indicating they can be valuable in the long-term conservation of a species, especially in rapidly changing landscapes (Safriel et al., 1994; Lesica and Allendorf, 1995; Hamilton and Eckert, 2007). Peripheral populations may also persist longer than those in the core, and in some cases have become the only remaining populations of once widespread species (e.g., Giant Panda Ailuropoda melanoleuca, California Condor Gymnogyps californianus; Channell and Lomolino, 2000).
Phrynocephalus versicolor (Toad Headed Agama) is a wide-ranging member of the Agamidae family that ranges throughout the desert environments of central and northern Asia (Wang and Fu, 2004). Their distribution covers portions of western China (Xinjiang, Gansu, Ningxia, and Inner Mongolia), extends northward into Mongolia, and includes some areas in Kazakhstan and southeastern Russia (Tuv region; Terbish et al., 2006a). Mongolia largely represents the northern edge of their global range (Ananjeva et al., 1997; Terbish et al., 2006b). Toad Headed Agamas are considered common and widespread and probably represent the most abundant vertebrate in some regions. Surveys in Ikh Nart Nature Reserve in Mongolia (Dornogobi Aimag) indicated that they can reach densities of up to 112/ha (average = 66/ha; Murdoch et al., 2010b); however, density varies considerably and has been reported as 37/ha in northern Mongolia (Rogovin et al., 2001) and 0.6/ha in Xinjing Zhunge'er Basin, China (Quan and Zhang, 2006).
Despite being common, we know relatively little about the behavior and ecology of the species or its ecological roles in the arid ecosystems of Asia. Toad Headed Agamas are arid adapted and occupy a variety of desert, semidesert, and steppe habitats including grasslands, forblands, and shrublands in flat or gently rolling terrain (Rogovin et al., 2001; Terbish et al., 2006a; Murdoch et al., 2010b, 2013). They are most active during the daytime hours and spend nights (and periods beyond their active temperature range) in subterranean burrows. During winter months, agamas hibernate; in Mongolia, the species typically emerges in April/May and is active until September/October. Toad Headed Agamas also represent an important food source for predators, including threatened species such as Lesser Kestrels (Falco naumanni) (Ganbold et al., 2017), and their burrows probably provide shelter for insects and other lizards and contribute to soil aeration and regeneration (Ananjeva et al., 1997; Terbish et al., 2006b).
The northern peripheral region of the range of Toad Headed Agamas is experiencing a relatively high rate of landscape change, mainly because of human activities. Mongolia is rapidly developing industry (including mining and other forms of resource extraction) and building infrastructure (including roads and railways), and large increases in livestock numbers affect habitat quality and availability for many species (Farrington, 2005; Berger et al., 2013; Ito et al., 2013; Batsaikhan et al., 2014). Climatic conditions have also changed, resulting in noticeable shifts in the distribution of some major ecosystems (Dagvadorj et al., 2009). Since 1940, average annual temperature has increased >2°C and precipitation has declined ∼0.1–0.2 mm/year in some steppe regions in the country (Dagvadorj et al., 2009). As Mongolia continues to change, “how species will respond” is key, especially those sensitive to temperature and habitat disturbance (Huey et al., 2012).
Our goal in this study was to examine how a variety of biotic and abiotic elements of a landscape influence the distribution of Toad Headed Agamas in a peripheral region of their global range. Our objective was to build an occupancy model that accounted for detection probability to provide a quantitative tool for estimating the impact of future landscape change on the species. Our broad hypothesis was that occupancy of Toad Headed Agamas is shaped by five factors: amount of vegetation cover, substrate type, vegetation structure, topographic ruggedness, and elevation. We examined a range of specific hypotheses in our modeling process that included different combinations of factors.
Materials and Methods
We conducted the study in central and western Mongolia including the aimags (or provinces) of Khovd, Gobi-Altai, Bayankhongor, Uvurkhangai, and Tuv (Fig. 1). The study area was a semiarid and desert steppe landscape that included three primary ecosystems: Mongolian Altai Mountain Steppe, Dzungarian Gobi Desert, and Gobi-Altai Mountainous Desert-Steppe (Grubov and Junatov, 1952). Study sites were in mostly semiarid and arid regions that had sparse vegetation cover. Among the primary habitat types included, forblands dominated by nongraminoid species such as Allium polyrrhizum, Anabasis brevifolia, Artemisia frigida, Caragana leucophloea, Nitraria sphaerocarpa, and Reaumuria songarica and grasslands dominated by graminoid species such as Agropyron cristatum, Festuca lenensis, Setaria viridis, Stipa glareosa, and Stipa gobica. Major soil types included brown desert steppe and grey-brown desert soils (Tsegmid, 1969). Soils were grayish brown, brown carbonaceous sand, and coarse gravel. The region was characterized by a cold semiarid and desert climate with annual average precipitation that typically varied from 119 to 220 mm, and average monthly temperature ranged from 19°C in July to −24°C in January.
We collected detection/nondetection data during surveys of 180 sites (Fig. 1). Sites were selected to maximize variability in ecosystem type, vegetation and habitat type, soil conditions, and geographic breadth within the known distribution of the species in central and western Mongolia (based on Terbish et al., 2006a). Sites were separated by a minimum of 10 km for independence. We did not select sites randomly because of the constraints of land use access, ownership, and status, but think that variability and geographic extent of sites reflected the range of conditions found at the periphery of the species distribution and that any biases were limited. We defined a site as a 25-m radius circular plot, and a survey involved walking in a zig-zag fashion through a given plot for 5 min and recording whether we detected an agama (Murdoch et al., 2013). We conducted two surveys of each site, with each survey separated by 1 h. This period was meant to ensure independence of our surveys, and we based the time on preliminary field observations of the species. We scored a detection as a ′1′ and nondetection as a ′0′ for each survey, resulting in a capture history for each site. Possible capture histories were 11, 10, 01, and 00. All surveys occurred during June and July and between 0800 h and 2100 h when agamas are typically active (Borkin and Semenov, 1986; Murdoch et al., 2013). At the start of each survey we measured the temperature (°C) and wind speed (m/s) using a handheld Kestrel 2500 weather meter (Nielsen-Kellerman Co., Boothwyn, Pennsylvania, USA). We also recorded the proportion of major habitat types (forbland and grassland), proportion of major soil types (sand and gravel), average vegetation height, and elevation. In addition, we estimated site ruggedness, which we indexed on a scale from 1 (least rugged, flat terrain) to 5 (most rugged, highly variable terrain).
We developed a model that predicted probability of occupancy from the detection/nondetection data. We used an occupancy modeling approach that estimates occupancy probability (ψ) while accounting for detection probability (p) using maximum likelihood methods (MacKenzie et al., 2002). Occupancy modeling involved developing a set of a priori models, fitting each to the data, and then using model selection techniques to determine the best model in the set (MacKenzie et al., 2006).
Our model set consisted of 70 total models that explored a range of covariate effects. Models in the set included those with single covariates, additive covariates, and interactions between covariates on occupancy probability and detection probability as well as a null model (Table 1; Appendix 1). Covariates included two habitat types (forb cover and grass cover), two substrate types (sand and gravel), vegetation structure (height), topography (ruggedness), and elevation (Table 1). Habitat types represented the two major land cover types in the region, and we predicted that forb-dominated areas would positively affect occupancy because they offer greater opportunities for vigilance against predators and concealment, which would probably not be the case for grass cover (Ananjeva et al., 1997; Rogovin et al., 2001; Terbish et al., 2006b). Similarly, we expected taller vegetation would limit vigilance, especially against aerial predators such as Lesser Kestrels and other diurnal raptors, and negatively influence occupancy (Ganbold et al., 2017). Soil type represents an important consideration for burrowing, and both sand and gravel are frequently used, so we assumed both would positively affect occupancy but to varying degrees (Shenbrot et al., 1991). We defined soil type based on top soil particle size: sand <2 mm and gravel 2–16 mm. Lastly we included ruggedness, which negatively influenced occupancy of Toad Headed Agamas elsewhere (Ikh Nart Nature Reserve, Dornogobi Aimag, Mongolia; Murdoch et al., 2013), and elevation to account for its potential effect. Our model set included models with each covariate as a single effect on occupancy, models with additive combinations of habitat type, soil type, and ruggedness (covariates we believed had the strongest effects) on occupancy, and interactions between habitat type, soil type, and ruggedness on occupancy (Appendix 1). Every model also included a parameter for detection probability. We then replicated the set twice: the first replicate had all of the same models but included the effect of temperature on detection probability—we modeled temperature as a polynomial (temperature + temperature2) because agamas are ectothermic and temperature does not affect activity linearly. The second replicate included all of the same models but included the effect of wind speed on detection probability (Appendix 1).
Variables used to describe occupancy probability (ψ) and detection probability (p) of Toad Headed Agamas (Phrynocephalus versicolor) in central and western Mongolia in 2016.
We used Akaike's Information Criterion (AIC) to rank each model and considered models with <2 ΔAIC units from the top-ranked model to have strong empirical support (Burnham and Anderson, 2002). We also calculated evidence ratios to quantify the relative support of models (ratio of the weight of the top ranking model to the weight of a given model; Burnham and Anderson, 2002), We made occupancy and detection probability predictions from a supported model using the logit-link function. We performed all analyses in Program Presence v11 (Hines, 2006).
Assessing Model Fit
The AIC provides a measure of the support of a given model relative to others in the set. To estimate how well the data fit the top-ranking model, however, we used a bootstrapping technique described by MacKenzie and Bailey (2004). Briefly, this involved generating 100 bootstrapped data sets predicted by the most-parameterized model in the set and calculating a goodness-of-fit chi-square statistic for each by comparing observed (bootstrapped) to expected capture histories. Chi-square values were plotted to create a distribution, then we estimated where the chi-square value for the actual observed data fit in this distribution. If the value fit in the tail of the distribution (P < 0.05), we considered the data not to have good fit with the data. Unless otherwise indicated, we report summary statistics as mean ± SD.
We detected agamas at 89 sites and during 141 of 360 (39.2%) surveys, resulting in a naïve occupancy of 49.4% (89 sites detected/180 total sites). Among the four possible capture histories for a site, we recorded 52 ′11′, 22 ′10′, 15 ′01′, and 91 ′00′ scores. Temperature varied from 16.0 to 33.1°C (mean = 25.1 ± 3.6) and wind speed varied from 0 to 16.9 m/s (mean = 4.9 ± 3.6) for surveys. Mean values for site covariates were: proportion of forb = 0.19 ± 0.13, proportion of grass = 0.19 ± 0.15, vegetation height = 13.0 ± 15.2 cm, proportion of sand substrate = 0.37 ± 0.33, proportion of gravel substrate = 0.50 ± 0.33, and ruggedness = 1.6 ± 0.9. Elevation varied from 743 to 2,146 m across sites and averaged 1,489 m.
Model selection results indicated that only one model had strong empirical support, so we considered this model to be the best in the set: ψ(Forb + Grass + Ruggedness) p(Wind) (Table 2). All other models, including those that were simpler with fewer parameters, had little to no support based on ΔAIC values and evidence ratios (Table 2). In the top model, forb, grass, and ruggedness all negatively influenced occupancy probability (Figs. 2, 3; Table 3). As the amount of forb and grass cover and topographic ruggedness increased, occupancy decreased (Figs. 2, 3; Table 3). The odds ratio for each covariate, which represents the predicted change in the odds of occupancy in response to a one-unit increase in a given covariate, was: forb = 0.0052 (95% confidence interval [CI]: 0.0001–0.1827), grass = 0.0006 (95% CI: 0.00001–0.0242), and ruggedness = 0.2369 (95% CI: 0.1186–0.4733). Detection probability in the top model was also negatively affected by wind: as wind speed increased, detection probability decreased (Fig. 4; Table 3). The odds ratio for wind was 0.8772 (95% CI: 0.8017–0.9598). The top model predicted detection probability to be 0.84 with no wind and 0.42 at high wind speed (15 m/s) (Fig. 4). All covariate parameter (β) estimates in the top model had 95% CIs that did not cross zero, suggesting that effects (and their directions) were meaningful (Table 3). Similarly, odds ratio CIs for each covariate did not include 1.0. The detection/nondetection data had good fit with the model. The chi-square statistic for the observed data was 0.39, and the probability of this value in the distribution of bootstrapped values was 0.51. Across all survey sites, mean occupancy predicted from the top model was 0.55 ± 0.02 SE (range = 0–0.93), higher than the naïve occupancy estimate.
Model selection results of occupancy probability (ψ) of Toad Headed Agamas (Phrynocephalus versicolor) for the 10 top-ranking models out of a set consisting of 70 total models based on detection/nondetection data collected in central and western Mongolia in 2016. Candidate models in the set considered the single, additive, and interacting effects of seven covariates on occupancy probability (proportion of forbland, grassland, sand substrate, and gravel substrate; average vegetation height; topographic ruggedness; and elevation) and two covariates on detection probability (p; temperature modeled as a polynomial: temperature + temperature2; and wind) at a site.
Parameter (β) estimates with standard errors (SE) and upper (UCI) and lower (LCI) 95% confidence intervals for the top-ranking model of occupancy of Toad Headed Agamas (Phrynocephalus versicolor) based on detection/nondetection data collected in central and western Mongolia in 2016.
Toad Headed Agamas range across the desert environments of central and northern Asia, including the Taklamakan Desert in northern China and Gobi Desert in northern China and southern Mongolia (Terbish et al., 2006a; Ananjeva, 2010). Despite occupying a large distribution and being considered a relatively common species, little ecological information exists on how landscape conditions shape patterns of occupancy, which makes comparisons of our results with other regions challenging. Our study along the periphery of the species range in Mongolia revealed that the amount of vegetation and degree of topographic ruggedness influence presence. Ideal conditions appear to be areas with little vegetation that are topographically flat or gently rolling. These conditions characterize much of the core areas of the species distribution in China and Mongolia, indicating that our model may more-broadly reflect occupancy, not only in the periphery but also in core areas.
Our results found average occupancy was 55% across our study sites, considerably lower than occupancy found at another site with similar habitat conditions in central-eastern Mongolia (Ikh Nart Nature Reserve), where occupancy was examined using a similar approach and estimated to be ∼85% (Murdoch et al., 2013). We expected occupancy across our sites to be similar to this site, and the difference may be explained by the much larger geographic scale of our study, which probably captured more variability in landscape conditions; the Ikh Nart study occurred in and around a single, relatively small nature reserve. For example, several of our survey sites occurred immediately along the edge of the known, current distribution, where drier deserts transition into other distinct ecosystems including mountain steppe and forests. These sites had low predicted occupancy, presumably because they included suboptimal conditions for agamas. Together, these two studies appear to confirm the known distribution of the species in Mongolia (estimated from expert opinion), which covers ∼40% of southern and western Mongolia (Terbish et al., 2006a). They also suggest that the relative quality of the landscape within this distribution is highly variable.
When considering occupancy at the site level, we expected that habitats (defined by vegetation cover and structure) would influence predation risk by providing some level of cover and concealment as well as opportunity for vigilance from aerial and ground-based predators. Instead, we found that both main habitats, grassland and forbland, exerted negative effects on occupancy, suggesting that agamas probably prefer more-open, unobstructed habitats for greater vigilance despite the trade-offs associated with greater cover and concealment. Toad Headed Agamas use subterranean burrows during the night and in unsuitable conditions (e.g., temperatures outside of their normal activity range), and individuals appear to maintain a network of several burrows (Murdoch et al., 2013). Agamas probably rely on these burrows as safe refuges when predators are detected, which may partly offset the need for vegetation for cover and concealment.
Topographic ruggedness also negatively affected occupancy, reflecting patterns of occupancy observed in the Ikh Nart study (Murdoch et al., 2013). That study examined the influence on agama occupancy of several habitat and landscape conditions similar to those in our study, including the amount of open plain vegetation types (grassland and forbland) and rugged terrain, and found that only one had a meaningful influence: highly rugged, rocky terrain (Murdoch et al., 2013). Agamas rarely occurred in rocky areas and lived mainly in flat, gently rolling terrain. Highly rugged terrain probably limits the ability of agamas to quickly flee from a predator to a burrow. Rugged, rocky terrain may also absorb, retain, and reflect temperature differently, which may make thermoregulation more difficult in the Mongolian landscape and provide substrate less suitable for burrowing. Our results differed from the Ikh Nart study, as agamas were also influenced by the presence of two other habitats. Together, results from both areas reflect the ecological variability of the species and indicate that occupancy patterns are not always consistent across a geographic range. Similar studies in other parts of the species distribution (e.g., core areas) should be conducted to provide a more complete picture of how landscape conditions affect distribution.
Other factors considered in our analysis, including soil type and elevation, did not impact occupancy. Soil type affects the distribution of Phrynocephalus and other arid-adapted lizards elsewhere and provides important resources such as burrowing substrate and particular prey communities (Shenbrot et al., 1991; Murdoch et al., 2013). Soil probably has some influence on Toad Headed Agamas, but that influence was less supported than were the other covariates in the top model. Observations during surveys indicated that agamas generally occupied and burrowed in both types of soil substrate relatively equally. We also found no influence of elevation, indicating that at least among the range of elevations of our survey sites, the species was not highly sensitive to this factor. In Mongolia, Toad Headed Agamas have been recorded from 600 to 2,000 m above sea level (Borkin et al., 1990; Ananjeva et al., 1997). We recorded no agamas at sites below 800 m or above 1,900 m.
Our analysis involved accounting for detection probability when estimating occupancy. Detection was most influenced by wind speed and less by temperature. Wind affects temperature experienced by agamas and reduces visibility, especially at ground level because of the movement of soil and sand particles. Temperature clearly has an effect on the activity of agamas because they are ectothermic (Shenbrot et al., 1991; Murdoch et al., 2013), but in our study area the effect of wind was more influential. Most of the surveys occurred within the temperature range of normal activity. Studies across a broader range of temperatures may reveal a greater effect of temperature on activity.
We considered a variety of biotic and abiotic factors in our analysis of occupancy which we thought, a priori, had the greatest potential to impact distribution; however, other factors may influence occupancy such as competition from other sympatric lizards (e.g., Eremias multiocellata), predator abundance, and food abundance. The effect of competition was probably negligible, as Toad Headed Agamas are typically dominant to other lizards in the study region (Terbish et al., 2006a,b), and we rarely recorded the presence of other species during our surveys, suggesting that other lizards were far less common. Aerial predators were also rarely seen during surveys and presumably were low in density, and ground-based predators such as foxes (Vulpes spp.) and badgers (Meles leucurus) infrequently consume the species (<5% of diet by volume; Murdoch and Buyandelger, 2010; Murdoch et al., 2010a). The abundance and distribution of food could impact agama distribution, but we were unable to estimate food sources in this study, partly because very little information exists on the diet of the species.
Effective conservation planning requires information on how landscape conditions shape the abundance and distribution of species. Mongolia is experiencing rapid economic change that has and will continue to affect the natural environment, especially of desert and steppe regions of the country (Batsaikhan et al., 2014). Changes have led to increases in livestock for meat and other goods, including cashmere for the global market (Reading et al., 2006; Berger et al., 2013). Livestock reduce vegetation cover through grazing and browsing that may benefit agamas, but also trample the vegetation, which affects soil conditions and represents a source of mortality for agamas (Rogovin and Semenov, 2004). Other forms of landscape change include the construction of roads, railways, and other infrastructure (Ito et al., 2013). Road and rail density was fairly low in our study region, but future developments are planned that could have localized effects on agama populations. Lastly, climate change represents a concern for species conservation. Mongolia has experienced a relatively high increase in average annual temperature, and reduction of precipitation in some regions, which has led to noticeable northward shifts of some major ecosystems (Dagvadorj et al., 2009). As Mongolia becomes warmer and drier, species distributions will probably change in response, and populations at the leading edge of distribution will be important sources of variability that may facilitate persistence and range shifts (Gibson et al., 2009). Toad Headed Agamas represent a relatively dominant reptile in arid regions of Mongolia and probably are not highly vulnerable to population declines (Terbish et al., 2006a), but may be susceptible to changes (positive or negative) in distribution as conditions change because they are sensitive to conditions in the landscape (in our case, habitat and ruggedness).
Our occupancy model represents a model of distribution that provides a quantitative tool for estimating how occupancy will change under different scenarios of landscape change which may be helpful to wildlife and ecosystem managers in desert and steppe communities. In some instances, scenario-based assessments of distribution changes have been used for other species in the country that led to changes in resource management (e.g., Argali Sheep, Ovis ammon, and Corsac Fox, Vulpes corsac; Lkhagvasuren et al., 2016; Murdoch et al., 2017). Our model provides a means of assessing localized impacts, such as habitat loss or alteration from resource extraction activities, and broader-scale impacts because of habitat shifts associated with climate change on the species that may be helpful for conservation planning.
The project received financial support the Mongolian Ministry of Education, Culture, Science and Sports supported by the Asian Development Bank (Grant No. L2766:HERP) within the framework of higher education reform. We also thank the support of the Mongolian National University of Education and National University of Mongolia and Batbayar, Solongo, and Javkhlanzol for their fieldwork assistance. JDM received support from the Core Fulbright US Scholar Program for part of this work. The Mongolian Ministry of Nature, Environment and Tourism and General Authority for Border Protection granted permission for the study under permits: 2016.06.10 #09_3657 and 2016.06.02 #0445.
Model set used to evaluate occupancy probability (ψ) of Toad Headed Agamas (Phrynocephalus versicolor) based on detection/nondetection data collected from 180 sites in central and western Mongolia in 2016. Models below include 23 base models with single covariates, additive combinations of covariates, and interactions of covariates. Each of these models included a parameter for detection probability. We then created two replicates of the set: the first included all the models with the addition of the effect of temperature (modeled as a polynomial: temperature + temperature2) on detection probability; the second included all of the models with the addition of the effect of wind on detection probability. Lastly, we also included a null model with no effects on occupancy or detection. Total models in the set was 70.