Mapping Habitat for Multiple Species in the Desert 

Southwest

Open-File Report 2014–1134

U.S. Department of the Interior

U.S. Geological Survey

Mapping Habitat for Multiple Species in the Desert 

Southwest

By Richard D. Inman, Kenneth E. Nussear, Todd C. Esque, Amy G. Vandergast, Stacie A. Hathaway, Dustin A. 

Wood, Kelly R. Barr, and Robert N. Fisher 

Open-File Report 2014–1134

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2014 

For more information on the USGS—the Federal source for science about the Earth,

its natural and living resources, natural hazards, and the environment—visit 

http://www.usgs.gov or call 1–888–ASK–USGS

For an overview of USGS information products, including maps, imagery, and publications,

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Suggested citation:

Inman, R.D., Nussear, K.E., Esque, T.C., Vandergast, A.G., Hathaway, S.A., Wood, D.A., Barr, K.R., and Fisher, R.N, 

2014, Mapping habitat for multiple species in the Desert Southwest: U.S. Geological Survey Open-File Report 2014-

1134, 92 p., http://dx.doi.org/10.3133/ofr20141134.. 

ISSN 2331-1258 (online)

Any use of trade, firm, or product names is for descriptive purposes only and does not imply 

endorsement by the U.S. Government. 

Although this information product, for the most part, is in the public domain, it also may contain

copyrighted materials as noted in the text. Permission to reproduce copyrighted items must be 

secured from the copyright owner. 
iii

Contents

Abstract ......................................................................................................................................................................... 1

Introduction.................................................................................................................................................................... 1

Study Area ................................................................................................................................................................. 2

Study Species............................................................................................................................................................ 2

Habitat Suitability Modeling ........................................................................................................................................... 4

Results..........................................................................................................................................................................12

Environmental Covariates Summary.........................................................................................................................12

Elevation ...............................................................................................................................................................17

Winter Precipitation ...............................................................................................................................................17

Surface Texture.....................................................................................................................................................17

Topographic Position.............................................................................................................................................18

Seasonal Thermal Change....................................................................................................................................18

Seasonal Greenness Timing .................................................................................................................................19

Slope.....................................................................................................................................................................19

Soil Adjusted Greenness.......................................................................................................................................19

Soil Water Stress...................................................................................................................................................19

Sandy Surfaces.....................................................................................................................................................20

Mountain Bases.....................................................................................................................................................20

Washes .................................................................................................................................................................20

Alluvial Fan Proximity ............................................................................................................................................20

Rocky Surfaces .....................................................................................................................................................21

Seasonal Greenness Ratio....................................................................................................................................21

Cliffs ......................................................................................................................................................................21

Vegetation Greenness Coefficient of Variation......................................................................................................21

Habitat Model Summaries.........................................................................................................................................22

Anaxyrus punctatus...............................................................................................................................................22

Chaetodipus penicillatus .......................................................................................................................................29

Chionactis occipitalis .............................................................................................................................................33

Crotaphytus bicinctores and C. collaris .................................................................................................................38

Dipsosaurus dorsalis .............................................................................................................................................43

Homalonychus selenopoides and H. theologus.....................................................................................................48

Lichanura trivirgata................................................................................................................................................51

Perognathus longimembris....................................................................................................................................56

Plestidon gilberti....................................................................................................................................................61

Sceloporus magister..............................................................................................................................................65

Thomomys bottae..................................................................................................................................................70

Uma scoparia ........................................................................................................................................................75

Xantusia vigilis.......................................................................................................................................................80

Discussion ....................................................................................................................................................................85

References Cited..........................................................................................................................................................86
iv

Figures

Figure 1. Representative species response curves across hypothetical environmental covariate gradients .............. 14

Figure 2. Map showing standard error (σ/N) among the selected models for Anaxyrus punctatus 

(Red-spotted Toad)...................................................................................................................................................... 24

Figure 3. Map showing maximum habitat suitability among the selected models for Anaxyrus punctatus 

(Red-spotted Toad)...................................................................................................................................................... 25

Figure 4. Generalized response curves and histograms for Anaxyrus punctatus (Red-spotted Toad) for 

environmental covariates............................................................................................................................................. 28

Figure 5. Map showing average habitat suitability among the selected models for Chaetodipus penicillatus

(Desert Pocket Mouse)................................................................................................................................................ 30

Figure 6. Generalized response curves and histograms for Chaetodipus penicillatus (Desert Pocket Mouse) for 

environmental covariates............................................................................................................................................. 32

Figure 7. Map showing standard error (σ/N) among the selected models for Chionactis occipitalis 

(Shovel-nosed Snake). ................................................................................................................................................ 34

Figure 8. Map showing average habitat suitability among the selected models for Chionactis occipitalis 

(Shovel-nosed Snake). ................................................................................................................................................ 35

Figure 9. Generalized response curves and histograms for Chionactis occipitalis (Shovel-nosed Snake) for 

environmental covariates............................................................................................................................................. 37

Figure 10. Map showing standard error (σ/N) among the selected models for the Crotaphytus spp. group 

(collared lizards) .......................................................................................................................................................... 39

Figure 11. Map showing average habitat suitability among the selected models for the Crotaphytus spp. group

(collared lizards) .......................................................................................................................................................... 40

Figure 12. Generalized response curves and histograms for the Crotaphytus spp. group (collared lizards) for 

environmental covariates............................................................................................................................................. 42

Figure 13. Map showing standard error (σ/N) among the selected models for Dipsosaurus dorsalis 

(Desert Iguana)............................................................................................................................................................ 44

Figure 14. Map showing average habitat suitability among the selected models for Dipsosaurus dorsalis (Desert 

Iguana) ........................................................................................................................................................................ 45

Figure 15. Generalized response curves and histograms for Dipsosaurus dorsalis (Desert Iguana) for

environmental covariates............................................................................................................................................. 47

Figure 16. Generalized response curves and histograms for Homalonychus genus (Homalonychus spiders) for 

environmental covariances.......................................................................................................................................... 49

Figure 17. Map showing average habitat suitability among the selected models for Homalonychus spiders.............. 50

Figure 18. Map showing standard error (σ/N) among the selected models for Lichanura trivirgata (Rosy Boa)......... 52

Figure 19. Average habitat suitability among the selected models for Lichanura trivirgata (Rosy Boa)...................... 53

Figure 20. Generalized response curves and histograms for Lichanura trivirgata (Rosy Boa) for environmental 

covariates .................................................................................................................................................................... 55

Figure 21. Map showing standard error (σ/N) among the selected models for Perognathus longimembris 

(Little Pocket Mouse)................................................................................................................................................... 57

Figure 22. Average habitat suitability among the selected models for Perognathus longimembris (Little Pocket 

Mouse)......................................................................................................................................................................... 58

Figure 23. Generalized response curves and histograms for Perognathus longimembris (Little Pocket Mouse) for 

environmental covariates ............................................................................................................................................. 60

Figure 24. Map showing average habitat suitability among the selected models for Plestidon gilberti 

(Gilbert’s Skink) ........................................................................................................................................................... 62

Figure 25. Generalized response curves and histograms for Plestidon gilberti (Gilbert’s Skink) for environmental 

covariates .................................................................................................................................................................... 64
v 

Figure 26. Map showing standard error (σ/N) among the selected models for Sceloporus magister 

(Desert Spiny Lizard)................................................................................................................................................... 66

Figure 27. Average habitat suitability among the selected models for Sceloporus magister (Desert Spiny Lizard) .... 67

Figure 28. Generalized response curves and histograms for Sceloporus magister (Desert Spiny Lizard) for 

environmental covariates............................................................................................................................................. 69

Figure 29. Map showing standard error (σ/N) among the selected models for Thomomys bottae (Botta’s Pocket 

Gopher). ...................................................................................................................................................................... 71

Figure 30. Map showing average habitat suitability among the selected models for Thomomys bottae (Botta’s Pocket 

Gopher) ....................................................................................................................................................................... 72

Figure 31. Generalized response curves for Thomomys bottae (Botta’s Pocket Gopher) 

for environmental covariates........................................................................................................................................ 74

Figure 32. Map showing standard error (σ/N) among the selected models for Uma scoparia

(Mojave Fringe-toed Lizard)......................................................................................................................................... 76

Figure 33. Map showing average habitat suitability among the selected models for Uma scoparia 

(Mojave Fringe-toed Lizard)......................................................................................................................................... 77

Figure 34. Generalized response curves and histograms for Uma scoparia (Mojave Fringe-toed Lizard) for 

environmental covariates............................................................................................................................................. 79

Figure 35. Map showing standard error (σ/N) among the selected models for Xantusia vigilis 

(Desert Night Lizard).................................................................................................................................................... 81

Figure 36. Map showing average habitat suitability among the selected models for Xantusia vigilis 

(Desert Night Lizard).................................................................................................................................................... 82

Figure 37. Generalized response curves and histograms for Xantusia vigilis (Desert Night Lizard) for

environmental covariates............................................................................................................................................ 84

Tables

Table 1. Species for which habitat models were derived, the number of observation records used, the Area

Under the receiver operator characteristic Curve (AUC), the Boyce Index (BI), and the resulting amount of 

suitable habitat predicted from the averaged model for each species ........................................................................... 3

Table 2. Environmental covariates used in developing habitat suitability models (mean, standard deviation, 

minimum and maximum for study area)......................................................................................................................... 5

Table 3. Correlation coefficients for positive association (lower left) and significance (p-values; upper right) "among 

habitat suitability models.............................................................................................................................................. 13

Table 4. Relative contributions, prevalence, and average importance for environmental covariates included in 

habitat models ............................................................................................................................................................. 15
vi

Conversion Factors and Datums

Conversion Factors

Multiply By To obtain

Length

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

kilometer (km) 0.5400 mile, nautical (nmi) 

meter (m) 1.094 yard (yd) 

Area

square kilometer (km2) 247.1 acre

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Temperature in degrees Fahrenheit (°F) may be converted to degrees Celsius (°C) as follows:

°C=(°F-32)/1.8

Datums

Vertical coordinate information is referenced to the North American Vertical Datum of 1988 (NAVD 88). 

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83). 

Elevation as used in this report, refers to distance above the vertical datum.
1 

Mapping Habitat for Multiple Species in the Desert 

Southwest

By Richard D. Inman, Kenneth E. Nussear, Todd C. Esque, Amy G. Vandergast, Stacie A. Hathaway, Dustin A. 

Wood, Kelly R. Barr, and Robert N. Fisher

Abstract

Many utility scale renewable energy projects are currently proposed across the Mojave 

Ecoregion. Agencies that manage biological resources throughout this region need to understand the 

potential impacts of these renewable energy projects and their associated infrastructure (for example, 

transmission corridors, substations, access roads, etc.) on species movement, genetic exchange among 

populations, and species’ abilities to adapt to changing environmental conditions. Understanding these 

factors will help managers select appropriate project sites and possibly mitigate for anticipated effects of 

management activities. We used species distribution models to map habitat for 15 species across the 

Mojave Ecoregion to aid regional land-use management planning. Models were developed using a 

common 1 × 1 kilometer resolution with maximum entropy and generalized additive models. 

Occurrence data were compiled from multiple sources, including VertNet (http://vertnet.org/index.php), 

HerpNET (http://www.herpnet.org), and MaNIS (http://manisnet.org), as well as from internal U.S. 

Geological Survey databases and other biologists. Background data included 20 environmental 

covariates representing terrain, vegetation, and climate covariates. This report summarizes these 

environmental covariates and species distribution models used to predict habitat for the 15 species 

across the Mojave Ecoregion. 

Introduction

Many species of endemic plants and animals inhabit the Desert Southwest of the United States

and compete for resources and space in a rapidly changing landscape. Some of these species are 

federally or locally listed as sensitive, threatened or endangered in the region, primarily due to loss and 

fragmentation of habitat. Continued pressures from urbanization, utility-scale renewable energy 

development, altered landscape conditions and climate space resulting from climate change may cause 

additional declines in biodiversity and species abundance. When tasked with managing biological 

resources throughout this region, agency personnel must understand the potential impacts of these 

pressures on the biological communities under their stewardship to make appropriate decisions. 

Tools, such as habitat suitability models, can help inform managers about some of the potential 

impacts of these future pressures by identifying areas where one or more key species exist, and can 

provide insight to potential outcomes of land-use changes. Habitat suitability models, or species 

distribution models (SDMs), as they are often called, can be quantitative, spatially explicit models that 

relate species occurrence to environmental predictors thought to influence or define the suitability of 

habitat for the species (Franklin 2010). These models can be ecologically complex or very simple, and 

can be modeled at many scales. Habitat suitability models also are used for tasks such as the design of 

conservation and monitoring programs, evaluating the efficacy or potential effects of proposed 
2 

management actions, and in recovery planning for threatened or endangered species (Graham and 

others, 2004; Elith and Leathwick, 2009). Additionally, SDMs are being used for species richness 

assessments (Elith and others, 2006; Franklin, 2010) and the evaluation of potential changes in species 

distributions resulting from climate change (Araújo and Williams, 2000; Guisan and Zimmermann, 

2000; Pearson and Dawson, 2003; Franklin, 2010).

Habitat models are now often built with presence-only data (as absence data for many species 

are often unavailable) because of the increased availability of historical records of species occurrence 

(Elith and Leathwick, 2007; Frey, 2009; Newbold, 2010), recent adaptations of software that use only 

presence data to model habitat (for example, MaxEnt; Phillips and others, 2006), and studies using 

different methods of generating pseudo-absence that are used in more traditional modeling frameworks 

(Engler and others, 2004). Although models using presence-only observation records are generally 

thought to be less accurate in defining species-environment relationships (Guisan and Zimmermann, 

2000; Austin, 2002), new methods can be very accurate at discriminating habitat from non-habitat 

across many spatial scales (Elith and Graham, 2009). Our objectives were to show the geographic 

distributions of habitat for a broad taxonomic range of species and for species-environment 

relationships. We used two methods for deriving presence-only SDMs, and averaged them for each 

species to derive predictions of habitat suitability across our study area. This report provides the 

methods and results for the habitat models developed as part of a scientific response to management 

needs for evaluating areas where species and their potential habitats overlap in context of genetic 

diversity and renewable energy development throughout the Mojave Desert (Vandergast and others, 

2013). 

Study Area

Our study area included the Mojave Desert (as defined by the U.S. Environmental Protection 

Agency (EPA) Mojave Basin and Range Level III Ecoregion [Commission for Environmental 

Cooperation Working Group, 1997]) and a 10 km buffer area around its periphery. The total area, 

150,000 km2

, includes elevations ranging from 146 m below sea level (Death Valley, Calif.) to more 

than 3,633 m above sea level (Mount Charleston, Nev.), and is characterized by broad valleys separated 

by rugged mountain ranges and alluvial fan complexes. Precipitation ranges from 100 to 350 mm/yr and 

shows a strong longitudinal gradient with most precipitation in the Western Mojave Desert occurring in 

the winter, and most precipitation in the eastern and southern parts occurring in the late summer and 

winter (Turner, 1982). Summer precipitation, however, is highly variable and normally does not exceed 

winter precipitation. Extreme variation in air temperature is evident throughout the year with observed 

temperatures well over 54 °C in summer months and below 0 °C in winter months. Vegetation 

communities vary widely from barren salt flats and playas at the lowest elevations, through a variety of 

shrublands, chaparral, woodlands, forests, and a fractional amount of alpine habitat on the highest peaks 

(Brown, 1982; Sawyer and others, 2009). 

Study Species

The 15 species covered in this report are: Anaxyrus punctatus (Red-spotted Toad), Chaetodipus 

penicillatus (Desert Pocket Mouse), Chionactis occipitalis (Shovel-nosed Snake), Crotaphytus 

bicinctores (Great Basin Collared Lizard), Crotaphytus collaris (Eastern Collared Lizard), Dipsosaurus 

dorsalis (Desert Iguana), Homalonychus selenopoides and Homalonychus theologus (Homalonychus

spiders), Lichanura trivirgata (Rosy Boa), Perognathus longimembris (Little Pocket Mouse), 

Plestiodon gilberti (Gilbert’s Skink), Sceloporus magister (Desert Spiny Lizard), Thomomys bottae

(Botta’s Pocket Gopher), Uma scoparia (Mojave Fringe-toed Lizard), and Xantusia vigilis (Desert 
3 

Night Lizard). Because of the recent reorganization of the Crotaphytus genus, we assumed that many of 

the observation records for C. collaris were incorrectly classified (many C. collaris are now C. 

bicinctores). However, we assumed that in terms of habitat, few, if any, differences would be observed 

between the species within our study area. We therefore combined the two species into a single taxon

for modeling. We also combined the two species of Homalonychus spiders into a single taxon to 

increase sample size, thus reducing the 15 species to 13 taxa for which we modeled habitat (table 1). In 

general, these 13 taxa represented small, ground dwelling animals, which may be impacted by 

renewable energy development at a population scale (for example, potential loss of populations or loss 

of connectivity among populations). 

Table 1. Species for which habitat models were derived, the number of observation records used, the Area Under 

the receiver operator characteristic Curve (AUC), the Boyce Index (BI), and the resulting amount of suitable habitat 

predicted from the averaged model for each species. 

Species

(abbreviation)

Occurrence in 

model AUC BI

Suitable

habitat

(square kilometers)

Anaxyrus punctatus (ANPU) 150 0.83 0.89 87,510

Chaetodipus penicillatus (CHPE) 99 0.87 0.22 53,988

Chionactis occipitalis (CHOC) 231 0.79 0.64 85,397

Crotaphytus bicinctores/collaris

(CRBI/CRCO) 88 0.77 0.58 100,293

Dipsosaurus dorsalis (DIDO) 343 0.82 0.59 95,914

Homalonychus selenopoides/theologus

(HOTH/HOSE) 20 0.82 0.21 62,058

Lichanura trivirgata (LITR) 75 0.91 0.87 59,521

Perognathus longimembris (PELO) 302 0.78 0.75 98,515

Plestidon gilberti (PLGI) 83 0.86 0.59 60,983

Sceloporus magister (SCMA) 472 0.72 0.33 111,731

Thomomys bottae (THBO) 186 0.83 0.51 104,675

Uma scoparia (UMSC) 78 0.90 0.54 30,127

Xantusia vigilis (XAVI) 332 0.78 0.05 106,994
4 

Habitat Suitability Modeling

We compiled observation records from online public databases for museum collections 

including VertNet (http://vertnet.org/index.php), HerpNET (http://www.herpnet.org), and MaNIS 

(http://manisnet.org), as well as other sources (including the California Natural Diversity Database, 

Joshua Tree National Park, Edwards Air Force Base, the Conservation Biology Institute, internal U.S. 

Geological Survey [USGS] databases, and other biologists). Because the taxonomy for many of our 

species was recently revised, we used both historical and current scientific names during searches. 

Duplicate records and those lacking sufficient locality information (that is those with no site name or 

GPS coordinates, or estimated precision error of greater than 1.5 km) were removed. For each species, 

we developed conceptual models of habitat needs to identify the environmental covariates thought to 

influence their geographic distribution. We considered a total of 20 environmental covariates (table 2), 

which were derived prior to developing habitat models. Prior to creating each of the layers, a spatial 

mask was used to eliminate impermeable surfaces, such as urban areas and extensive roadways, as well 

as reservoirs, lakes, and other major water bodies, as all of our species were terrestrial obligates. This 

mask was derived from the National Land Change Database (NLCD) 2006 Percent Developed

Imperviousness layer, downloaded from the Multi-Resolution Land Characteristics Consortium (Fry and 

others, 2011) and the National Hydrography Dataset (at a 1:24,000 scale; Simley and Carswell, 2009). 

Models were conducted using a common 1 × 1 km resolution raster grid. Grid cells with more than 20

percent of their surface area covered by at least 20 percent imperviousness were categorized as urban. 
5 

Table 2. Environmental covariates used in developing habitat suitability models (mean, standard deviation, minimum and maximum for study area).

Environmental 

covariate

layer

Description Mean Standard 

deviation Minimum Maximum

alluvial fan 

proximity

Alluvial fans were defined as areas with less than 10 degrees slope and topographic position less than 11.5 at 

a 30 m spatial resolution. We used a fine scale filter to remove areas that were comprised of three or fewer 30 

m cells, and calculated Euclidean distance to any 30 m cell defined as an alluvial fan. The value given to each 

1 km cell was the mean of all Euclidean distances for the 30 m cells contained within it.

4,717.7 5,098.2 0.0 39,221.6

cliffs Cliffs were defined at a 30 m spatial scale as areas with slope greater than 25 degrees. The value given to each 

1 km cell was the percent of 30 m cells contained within it defined as cliffs. Prior to use in habitat models, this 

layer was normalized using an arcsine-square root transformation.

0.07 0.17 0.0 1.0

elevation Elevation in meters above sea level, derived from a 30 m National Elevation Dataset Digital Elevation Model 

(USGS). The value given to each 1 km cell is the average of all 30 m cells within it.

1,012.51 473.44 -84.8 4,061.8

mountain bases Mountain bases were defined as 30 m cells that (1) had less than 10 degrees slope, (2) had a profile curvature 

value between -6.0 × 10-4 and -1.8 × 10-4, and (3) surface roughness between 0.7 × 10-8 and 5.8 × 10-8. Profile 

curvature is the curvature in the direction of steepest slope and is expressed as meters-1 (a curvature of 0.05 

corresponds to a radius of curvature of 20 m). Negative values are associated with concave terrain features, 

such as canyons and mountain bases. Profile curvature and surface roughness thresholds were derived as the 

mean ±1 standard deviation of values at 50 locations defined a priori as mountain bases. The value given to 

each 1 km cell is the percent of 30 m cells contained within it that were defined as mountain bases. Prior to 

use in habitat models, this layer was normalized using an arcsine-square root transformation.

0.05 0.08 0.0 0.9

rocky surfaces Rocky surfaces were identified from summer night time Land Surface Temperature (LST) and elevation at a 

spatial scale of 250 m. MODIS MOD11A1 nighttime LST 8-day global data were downscaled to 250 m and 

regressed against elevation to remove the elevation trend. Residuals were categorized into "rock" and "nonrock" categories using 100 known sample locations in Nevada and California. The MODIS MOD11A1 Land 

Surface Temperature (LST) 8-day Global 1 km data were averaged for 6 periods during the summers of 2001

– 2010 nighttime LST. The dates used to define the summer periods were July 12, July 20, July 28, August 5, 

August 13, and August 21. The averaged LST product was downscaled by regressing using a Generalized 

Additive Model (GAM) against seven environmental predictor covariates, all of which were derived at a 250 

m spatial scale. These covariates included: solar beam irradiance (W/km2), average of the black sky albedo 

for the same time periods as the LST, maximum horizon level, elevation, aspect, slope and an average of the 

Normalized Difference Vegetation Index (NDVI) for the same time periods. The value given for each 1 km 

cell is the percentage covered by 250 m cells identified as rock. Prior to use in habitat models, this layer was 

transformed using the arcsine-square root transformation.

0.11 0.23 0.0 1.0
6 

Environmental 

covariate

layer

Description Mean Standard 

deviation Minimum Maximum

sandy surfaces Sandy surfaces were identified from summer night time Land Surface Temperature (LST) and elevation at a 

spatial scale of 250 m. MODIS MOD11A1 nighttime LST 8-day global data were downscaled to 250 m and 

regressed against elevation to remove the elevation trend. Residuals were categorized into sand and non-sand 

categories using 100 known sample locations in Nevada and California. The MODIS MOD11A1 Land 

Surface Temperature (LST) 8-day Global 1 km data were averaged for six periods during the summers of 

2001–2010 night LST. The dates used to define the summer periods were July 12, July 20, July 28, August 5, 

August 13, and August 21. The averaged LST product was downscaled by regressing using a Generalized 

Additive Model (GAM) against seven environmental predictor covariates, all of which were derived at a 250 

m spatial scale. These covariates included: solar beam irradiance (W/km2

), average of the black sky albedo for 

the same time periods as the LST, maximum horizon level, elevation, aspect, slope and an average of the 

Normalized Difference Vegetation Index (NDVI) for the same time periods. The predictions from this GAM 

were mapped at a 250 m spatial scale. We used elevation in a separate regression model against the 

predictions from the previous GAM model to remove the trend of elevation on surface temperature. The 

residuals from this regression were mapped, and provided an estimate of the deviation in night LST from what 

would be expected by elevation alone. This deviation captures much of the variability in surface temperature 

due to substrate type and texture. We used known sample locations to derive a threshold in the mapped 

residuals to assign each 250 m cell as either “sand” or “not-sand”. The value given for each 1 km cell is the 

percentage covered by 250 m cells labeled as sand. Prior to use in habitat models, this layer was transformed 

using the arcsine-square root transformation.

0.16 0.35 0.0 1.0

seasonal 

greenness timing

Seasonal greenness timing represents the period of the year when vegetation greenness reaches its maximum, 

and was defined as the date of the maximum value of the Normalized Difference Vegetation Index (MODIS 

MOD13Q1 Global 250 m 16-day index of NDVI; vegetation greenness). Dates ranged from early spring to 

late summer, and the average of 2001 to 2010 for the dates March 6, March 22, April 7, April 23, May 9, May 

25, June 10, June 26, and July 12 were used.

2.53 1.50 1.0 9.0

seasonal 

greenness ratio

Seasonal greenness contrasts spring and winter (growing and senescent) vegetation greenness from MODIS 

MOD13Q1 16-day L3 Global 250 m NDVI data. The yearly spring (average of March 22, April 7, April 23, 

and May 9) and winter (average of December 3, December 19, January 1, and January 17) NDVI were 

contrasted using the formula: (SpringNDVI-WinterNDVI)/(1 + (SpringNDVI+WinterNDVI)). The value 

given to each 1 km cell was the average of 2001 to 2010.

0.04 0.05 -1.0 1.0

seasonal thermal 

change

Seasonal thermal inertia was defined as the difference between summer and winter daytime MODIS 

MOD11A1 Land Surface Temperature (LST) 8-day Global 1 km for six periods in the summer and six 

periods in the winter during 2001–2010. The summer (July 12, July 20, July 28, August 5, August 13, and 

August 21) and winter (December 11, December 19, December 27, January 1, January 9, and January 17h) 

periods were contrasted to identify areas with greater seasonal fluctuation in LST. The value for each 1 km 

cell is the average of the yearly differences (°C) between the summer and winter daytime LST for 2001 to 

2010.

1,130.57 84.85 598.3 1,466.2

slope Slope in degrees was derived from a 30 m National Elevation Dataset Digital Elevation Model (USGS). The 

value given to each 1 km cell is the average of all 30 m cells contained within.

8.24 8.01 0.0 50.0
7

Environmental 

covariate

layer

Description Mean Standard 

deviation Minimum Maximum

soil adjusted 

greenness

Soil adjusted greenness is an index of vegetation greenness that has been modified to reduce the influence of 

soil color on vegetation greenness (Qi and others, 1994), and was derived from Landsat5 Thematic Mapper 

imagery at a 30 m spatial scale for fall (September–November) in 2010 for the NIR and red Landsat5 bands. 

The mSAVI was computed for fall (September–November) in 2010. The value for any given 1 km cell is the 

average of all mSAVI values for all of the overlapping 30 m cells. The value given to each 1 km cell is the 

average of all 30 m cells contained within it.

0.07 0.11 -0.8 0.8

soil water stress Soil water stress is an index of the in situ top layer soil moisture, and was derived as the mean of the 

Shortwave and Infrared Water Stress Index (SIWSI; Fensholt and Sandholt, 2003) for 46 dates (January 1, 

January 9, January 17, January 25, February 2, February 10, February 18, February 26, March 6, March 14, 

March 22, March 30, April 7, April 15, April 23, May 1, May 9, May 17, May 25, June 2, June 10, June 18, 

June 26, July 4, July 12, July 20 July 28, August 5, August 13, August 21, August 29, September 6, 

September 14, September 22, September 30, October 8, October 16, October 24, November 1, November 9, 

November 17, November 25, December 3, December 11, December 19, and December 27) in each year from 

2001 to 2010 using the MODIS MOD09A1 surface reflectance 8-day Global data product. Prior to calculating 

the SIWSI, each of the MODIS MOD09A1 surface reflectance 8-day Global 500 m data products was 

aggregated to a spatial scale of 1 km. The value for each 1 km cell in this layer is the mean across all years 

and dates.

0.10 0.06 -0.2 0.2

soil water stress 

coefficient of 

variation

The soil water stress coefficient of variation covariate layer represents the coefficient of variation in the 

Shortwave and Infrared Water Stress Index (SIWSI) for 46 dates in each year from 2001 to 2010. The value 

for each 1 km cell in this layer is the coefficient of variation across all years and dates.

0.77 1.31 0.1 10.0

summer surface 

temperature

Summer surface temperature was derived from the MODIS MOD11A1 daytime Land Surface Temperature 

(LST) 8-day Global 1 km data product for summer (July 12, July 20, July 28, August 5, August 13, and 

August 21) in 2001 to 2010. The value for each 1 km cell is the average of all dates and years.

surface texture Surface texture was modeled by taking the difference between the mean daytime and nighttime surface 

temperatures for 6 periods in the summer during the years of 2001–2010. Surface temperatures were obtained 

from MODIS MOD11A1 Land Surface Temperature 8-day Global 1 km data products for daytime and 

nighttime during the summer dates of July 12, July 20, July 28, August 5, August 13, and August 21. The 

value for each 1 km cell is the average of the yearly differences between the daytime and nighttime 

temperatures (°C * 100) for the dates listed above.

1,082.99 179.42 145.9 1,624.7

topographic 

position

We used methods described in Morre and others (1991; terrain wetness index) to calculate topographic 

position at a 90 m spatial scale using the area of the hillslope per unit contour length that drains into a given 

area, and the local surface topographic slope, which is estimated as (Δ vertical elevation) / (Δ horizontal 

distance) from a 30 m National Elevation Dataset Digital Elevation Model (USGS). The value for each 1 km 

cell is the average of the 90 m cells contained within it.

8.58 2.19 1.5 19.7
8

Environmental 

covariate

layer

Description Mean Standard 

deviation Minimum Maximum

vegetation 

greenness

Vegetation greenness was defined as the average of the Normalized Difference Vegetation Index (NDVI) 

across 23 dates (January 1, January 17, February 2, February 18, March 6, March 22, April 7, April 23, May 

9, May 25, June 10, June 26, July 12, July 28, August 13, August 29, September 14, September 30, October 

16, November 1, November 17, December 3, and December 19) in each year from 2001 to 2010, obtained 

from the MODIS MOD13Q1 16-day L3 Global NDVI data product. The value for each 1 km cell is the mean 

all dates and years.

1,761.06 810.44 -738.2 7,575.6

vegetation 

greenness 

coefficient of 

variation

Vegetation greenness coefficient of variation is the coefficient of variation in the Normalized Difference 

Vegetation Index (NDVI) across 23 dates (January 1, January 17, February 2, February 18, March 6, March 

22, April 7, April 23, May 9, May 25, June 10, June 26, July 12, July 28, August 13, August 29, September 

14, September 30, October 16, November 1, November 17, December 3, and December 19) in each year from 

2001 to 2010, obtained from the MODIS MOD13Q1 16-day L3 Global NDVI data product. The value for 

each 1 km cell is the mean all years.

0.21 0.18 0.0 10.0

washes The density of washes was estimated by identifying "wash like" features at 30 m spatial scale, and calculating 

the density (m/km2) in each 1 km cell. Wash-like features were delineated using methods described in 

Montgomery and Foufoula-Georgiou (1993) in areas defined as alluvial fans.

597.08 863.21 0.0 5,575.0

winter 

precipitation

Yearly average of the cumulative winter season precipitation, in 100th of an inch, for a 30-year period. We 

derived winter precipitation from average monthly winter (November–April) precipitation climatic normals 

calculated from conditions between 1961 and 1990. An inverse distance weighted interpolation method for 

climate covariates was used from Nalder and Wien (1998) and was conducted at a spatial scale of 90 m. The 

value for each 1 km cell is the average of all overlapping 90 m cells.

198.11 92.63 20.4 867.9
9 

Habitat models were created using Generalized Additive Models (GAM) with the package mgcv 

(version 1.7-6; Wood, 2011) in R (version 2.13.1; R Core Development Team, 2011) and maximum 

entropy (MaxEnt - version 3.3.3e, Phillips and others, 2006) methods from presence-only observation 

records at a spatial resolution of 1 km. Considerable differences exist between the GAM and MaxEnt 

models, each having its own advantages. The GAM modeling framework was chosen because of its ease 

of interpretation, and because GAM models can be evaluated by many well-understood metrics (such as 

generalized cross validation, AIC, BIC, and explained deviance; Wood, 2006), while MaxEnt models 

use machine learning methods and cannot easily be compared with information criteria, such as AIC or 

BIC. However, MaxEnt models have been shown to use correlated input data better than other methods, 

and can incorporate very complex non-linear response functions (Elith and others, 2011). A substantial 

difference between MaxEnt and GAM models is that MaxEnt models use large random selections (that 

is, 10,000) of the study area to characterize the range of background values for each of the covariates, 

while GAM models require the use of true or pseudo-absences (Mateo and others, 2010; Lobo and 

Tognelli, 2011; Stokland and others, 2011). The use of large random samples of background values

enables the MaxEnt software to compare the distributions of environmental covariates for the entire

study area to distributions of the same covariates at locations where the species of interest is present 

(Elith and others, 2011), and allows it to estimate the probability of an occurrence data value at a 

specific location given the environmental covariates. In contrast, the GAM modeling framework is often 

used as a Resource Selection Function (RSF), wherein the use of true absence data are required to 

capitalize on places where the organism of interest does not occur. We used the GAM modeling 

framework to conduct a population level assessment assessing used (presences) and unused areas 

(absences), but instead of using true absences, we derived an equal number of pseudo-absences (Hirzel 

and others, 2002; Phillips and others, 2009; Stokland and others, 2011) that were selected at random 

from all locations where a species was not observed. The binomial GAM models were created using 

Restricted Maximum Likelihood (REML) for each species, and we considered all habitat models that 

contained up to five of the environmental covariates from the suite of covariates identified in the 

conceptual models for each species. Interaction terms or covariates that were highly correlated were not 

included. We kept models simple to maximize their interpretability, to improve their generalization 

across the study area, and to minimize over-fitting (Collingham and others, 2000; Rushton and others, 

2004). We also reduced over-fitting by using the average of 20 replications of sampling without 

replacement of the observation records for each model. For many species, the observation records were 

spatially aggregated, and often confined to smaller portions of the study area, which is often the case for 

museum specimens (Mateo and others, 2010; Newbold, 2010). To avoid model bias caused by 

aggregated data, we used a geographically weighted resampling method for each of the 20 replications 

of sampling without replacement, and used a gamma control parameter of 1.3 (which is used for 

additional smoothing) (Wood, 2006) for all GAM models. The resampling method used a grid (10 km 

cell size) placed across the study area to identify areas where spatially aggregated observations 

occurred. In these cells, random selections of no more than three occurrences were used in each 

replication. 
10

To evaluate model performance, a random sample of 20 percent of the observation records and 

an equally sized random sample of the pseudo-absences were withheld each iteration. This dataset was 

used to derive the Area Under the receiver operating characteristic Curve (AUC) (Fielding and Bell, 

1997; Cumming, 2000) and Boyce Index (BI) (Boyce and others, 2002) test statistics. In addition, AIC 

(Anderson and others, 2000; Burnham and Anderson, 2002) and BIC (Posada and Buckley, 2004) test 

statistics for the 20 replications of each model were calculated, and the number of replications where a 

given model had a ΔAIC or ΔBIC of 0 were summed to identify the models that consistently had high 

AIC and BIC performance. We selected a candidate suite of approximately 9 models for each species by 

assessing the relative performance using the AIC, BIC, explained deviance (Wood, 2006), BI, AUC, and 

the standard deviation in BI and AUC across the 20 replicates. We selected models that consistently had 

the lowest AIC or BIC, had high discrimination ability (AUC), high deviation from randomness (BI), 

high explained-deviance, and low variance in BI and AUC statistics. Models were excluded if the 

environmental covariates showed unrealistic response curves or predictions; for example, the threshold 

feature used by MaxEnt models can yield rapidly reversing responses that are typically uncharacteristic 

of how animals interact with environmental gradients (Shcheglovitova and Anderson, 2013; Syfert and 

others, 2013).

The selected GAM models for each species (approximately 9) also were derived as MaxEnt 

models, thus creating as many as 18 possible models for each species. MaxEnt models were run with 20 

bootstrap replications, while randomly withholding 20 percent of the observation records for each 

replication. The performance and general plausibility of each of the MaxEnt models were evaluated 

using several methods, including AUC, BI, test gain (Phillips and Dudik, 2008), marginal and solitary 

response curves (Phillips and Dudik, 2008), as well as the environmental covariate contributions 

(Phillips and others, 2006). Marginal response curves illustrate how environmental covariates affect the 

predicted habitat suitability score in relation to the other predictors, and are used to assess speciesenvironment relationships for the environmental covariates (Phillips and Dudik, 2008). As with the 

evaluation of GAM models, MaxEnt models with unrealistic marginal and solitary response curves were 

not selected for review. The percent contribution was measured as the relative increase in training gain 

across all bootstrapped replications (Phillips and others, 2006), and provided a general measure of how 

important a given environmental covariate was in terms of the model fitting process. 

Finally, we reviewed each of the possible 18 models for each of the species by examining maps 

of the predicted habitat and response curves, and evaluating them in context of our knowledge of the 

species. We ranked the selected models for their overall accuracy and goodness-of-fit, and chose a 

subset of models from the 18 and averaged their outputs to create a final habitat suitability map for each 

species. For example, we only chose 14 of the 18 models for D. dorsalis to include in the final averaged 

habitat suitability map. One minor exception to this procedure was for A. punctatus, wherein we 

considered each of the selected models to be too restrictive in their range, and thus used the maximum 

of all models to create the final map. For all species, we rescaled each model to the range [0~1] prior to 

averaging to ensure that all models had equal weight in the final map using a linear data transformation. 
11

Selected areas were spotchecked to ensure that the map predictions were not unrealistic, 

especially in areas where few data were available, but where we (or our collaborators) had local 

knowledge. We also assessed the performance of the averaged model with the AUC statistic and BI 

plots derived from the entire set of observation records for each species. For AUC, we created a test set 

of absence points that was a random sample of background points with a size equal to the number of 

observation records, thus keeping the prevalence ratio at 1:1. The BI plots were evaluated to provide 

additional information on model robustness and deviation from randomness. We also assessed a 

measure of model agreement by creating a map of the standard error (σ/N) in habitat suitability scores 

among the selected models for each species. Although not a measure of model error, these agreement 

maps provide a general measure of model uncertainty across the landscape and show areas where the 

selected models agree (or disagree) with respect to habitat suitability. When compared with the final 

averaged habitat suitability maps, the agreement maps highlight relationships that indicate whether our 

selected models agree more on high or low suitability habitat. We assess the similarity in habitat 

suitability among the species by taking the Pearson’s correlation coefficient for positive association of 

habitat suitability across the entire study area for each pair of taxa. Correlation coefficients and their 

significance are reported, and correlations greater than 0.5 are highlighted. 

The relationship between habitat suitability and environmental covariates should ideally be 

explained by known ecological relationships between the organism of interest and its environment 

(Guisan and Zimmerman, 2000; Austin, 2002; Guisan and Thuiller, 2005), and thus we relate each of 

the environmental covariates back to literature where possible, and propose relevant explanations for 

most of the environmental covariates in our averaged habitat suitability maps. We also derive relative 

contributions for each of the covariates to assess the explanatory power of each predictor covariate for 

each species, and compare the strength of each of the covariates among species. The relative 

contributions for each covariate were derived by taking the average percent contribution (see Elith and 

others, 2011) across all MaxEnt models for a given species. Separately, we derive an index of average 

importance, which is the average of all relative contributions for each covariate layer. This index of 

average importance provides a measure with which to rank each of the environmental covariates. We 

also derive averaged response curves by holding each covariate at their mean of the study area, and 

calculating their linear contributions (Wood, 2006) from a saturated GAM model with all covariates 

represented in the selected models. These response curves illustrate the relationship between a given 

covariate and a species’ habitat suitability, and are shown with histograms of the covariate across the 

study area for comparison. Although our primary objective was to derive predictive models of habitat 

suitability, appropriate knowledge of habitat preference and ecological context can inform the model 

building process to reduce overfitting on spurious environmental covariates and can improve the 

interpretability of SDMs (Austin, 2002; Guisan and Thuiller, 2005). 

Finally, we reclassified the averaged habitat suitability map into a binary representation of 

suitable habitat using the fifth quantile of habitat suitability scores from the observation records. This 

definition of suitable habitat assumes that not all observation records are in areas with suitable habitat, 

but rather, a small proportion (that is, 5 percent) may actually be located in unsuitable habitat. This 

method results in a smaller proportion of the study area being classified as suitable habitat than the more 

commonly used method of selecting the minimum habitat score from the observation records (Liu and 

others, 2005). For each species, we report the total square kilometers of suitable habitat under this 

definition, and the habitat suitability score used for this threshold. 
12

Results

Our habitat models used between 20 (H. selenopoides & H. theologus spiders) and 472 (S. 

magister; table 1) observation records and had AUC scores that ranged from 0.72 (S. magister) to 0.91 

(L. trivirgata; table 1). The mean AUC score among all species was 0.82, suggesting high model 

performance overall (Elith and others, 2006). The average of the total area of suitable habitat for the 13 

taxonomic groups was 81,362 km2 (σ = 24,307 km2), and ranged from 30,127 km2 for U. scoparia to 

111,731 km2 for S. magister. Habitat suitability models for C. penicillatus, C. occipitalis, D. dorsalis, 

and U. scoparia were positively correlated with one another (r > 0.5; table 3). Similarly, L. trivirgata, S. 

magister, T. bottae, P. longimembris, and X. vigilis also were correlated with one another (r > 0.5; table 

3). P. gilberti and D. dorsalis were negatively correlated (r = -0.64, N = 1,000, p < 0.001), indicating a 

geographic separation for these two throughout the study area. The remaining species (A. punctatus, and 

the single taxonomic group for the two species of Crotaphytus) were not strongly correlated with any 

others. 

Species response curves for each of the environmental covariates varied substantially among the 

species, but between the two modelling frameworks (GAM and MaxEnt), remained mostly similar. The 

MaxEnt models often had more complex response curves with multiple thresholds and several minima 

and maxima, while the GAM response curves were (by design) more generalized, and were limited to 

four broad patterns: (1) increasing or decreasing, (2) single minima or maxima, (3) single threshold, and 

(4) double threshold (fig. 1). Species with correlated habitat maps tended to have similar species 

response curves, and environmental covariates that were present in multiple models for the same species 

were generally similar in shape.

Agreement maps (σ/N) for each species showed a wide range of patterns between model 

uncertainty and habitat suitability. The geographic areas showing the least agreement among models for 

each of the species tended to be those with high elevations or those in playas. This is likely due to 

several of the vegetation layers, which often had extreme values at high elevations, and on playas during 

periods when they were wet. None of the species included here are reported to occur on playas, and 

therefore we did not attempt to rectify these differences. Only two models were selected for C. 

penicillatus and P. gilberti, thus making the standard error among the selected models incalculable. We 

also did not calculate a standard error map for the Homalonychus spiders because we only used a single 

MaxEnt model to represent habitat. 

Environmental Covariates Summary

Of the 20 environmental covariate layers considered in this study, three [soil water stress 

coefficient of variation (CV), vegetation greenness, and summer surface temperature] were not included 

in any of the selected models for the species covered here. The remaining 17 covariate layers (alluvial 

fan proximity, cliffs, elevation, mountain bases, rocky surfaces, sandy surfaces, seasonal greenness 

timing, seasonal greenness ratio, seasonal thermal change, slope, soil adjusted greenness, soil water 

stress, surface texture, topographic position, vegetation greenness CV, washes, winter precipitation) 

were included in at least one model, although relative contributions (table 4) and response curves varied 

considerably. Each of the environmental covariates are discussed in order of relative contribution. 

Methods used to derive each of the environmental covariates are given briefly, and in more detail in 

table 2. 
13

Table 3. Correlation coefficients for positive association (lower left) and significance (p-values; upper right) among 

habitat suitability models. 

[Correlation coefficients greater than 0.5 are shown in two shades of gray, representing two associated groups. Two groups 

of species had similar habitat suitability distributions: Dipsosaurus dorsalis, Chionactis occipitalis, Chaetodipus penicillatus, 

and Uma scoparia (light gray); Sceloporus magister, Lichanura trivirgata, Thomomys bottae, Perognathus longimembris, 

and Xantusia vigilis (dark gray). Species abbreviations are shown in table 1]

ANPU CHOC CHPE CRBI/CRCO DIDO HOTH/HOSE LITR PELO PLGI SCMA THBO UMSC XAVI

ANPU 1.000 0.001 0.000 0.853 0.000 0.000 1.000 0.001 0.988 1.000 1.000 1.000

CHOC -0.47 0.000 1.000 0.000 0.107 0.035 0.000 1.000 0.574 0.000 0.000 0.000

CHPE 0.10 0.41 1.000 0.000 0.000 0.000 0.000 1.000 0.994 0.000 0.000 1.000

CRBI/CRCO 0.21 -0.51 -0.27 1.000 0.000 0.149 1.000 0.000 0.000 1.000 1.000 0.073

DIDO -0.03 0.70 0.68 -0.38 0.000 1.000 0.896 1.000 1.000 0.999 0.000 1.000

HOTH/HOSE 0.27 0.04 0.21 0.29 0.32 0.000 1.000 1.000 0.000 1.000 0.000 0.000

LITR 0.19 0.06 0.20 0.03 -0.13 0.25 0.000 0.000 0.000 0.000 1.000 0.000

PELO -0.41 0.38 0.11 -0.24 -0.04 -0.11 0.19 0.000 0.000 0.000 0.147 0.000

PLGI 0.07 -0.41 -0.15 0.21 -0.64 -0.32 0.36 0.38 0.000 0.000 1.000 0.006

SCMA -0.07 -0.01 -0.08 0.21 -0.27 0.33 0.51 0.29 0.17 0.000 1.000 0.000

THBO -0.37 0.27 0.12 -0.26 -0.10 -0.31 0.18 0.82 0.45 0.14 0.368 0.000

UMSC -0.16 0.55 0.62 -0.22 0.67 0.15 -0.19 0.03 -0.36 -0.27 0.01 1.000

XAVI -0.32 0.25 -0.29 0.05 -0.26 0.12 0.30 0.47 0.13 0.74 0.30 -0.30
14

Figure 1. Representative species response curves across hypothetical environmental covariate gradients. Species 

response curves generally fell into one of the following broad patterns: (A) increasing or decreasing, (B) single 

minima or maxima, (C) single thresholded, and (D) double thresholded. Positive (A, C and D; solid) and negative 

(A, C, and D; dashed) responses were less common than unimodal with single minima (solid) or maxima (dashed). 

Species response curves represent the generalized relationship between an environmental covariate and species 

occurrence.
15

Table 4. Relative contributions, prevalence, and average importance for environmental covariates included in 

habitat models. [Relative contributions were derived by taking the average percent contribution across all MaxEnt models for a given 

species. The average of these relative contributions was used as the average importance for each species. Species 

abbreviations are shown in table 1. Prevalence represents the number of species models each covariate was included in which 

each covariate appeared] 

Environmental 

covariate

layer

ANPU CHOC CHPE CRBI/CRCO DIDO HOTH/HOSE LITR Prevalence Average 

importance

elevation 19 27 45 10 61 0 33 11 26.3

winter 

precipitation

6 0 8 0 14 0 22 8 18.0

surface texture 40 42 22 4 18 33 0 10 14.6

topographic 

position

0 22 0 29 2 29 4 10 12.8

seasonal 

thermal 

change

0 0 23 0 2 0 22 8 6.3

seasonal 

greenness 

timing

0 6 0 12 2 23 3 9 4.6

slope 0 0 0 35 0 0 1 5 4.5

soil adjusted 

greenness

32 0 0 0 0 0 0 1 2.5

soil water 

stress

3 2 0 5 0 0 0 8 2.4

sandy surfaces 0 0 0 0 1 0 0 2 2.1

mountain 

bases

0 0 2 1 0 0 5 7 2.1

washes 0 0 0 0 0 15 0 8 1.7

alluvial fan 

proximity

0 0 0 4 0 0 8 2 0.9

rocky surfaces 0 0 0 0 0 0 0 1 0.4

seasonal 

greenness 

ratio

0 0 0 0 0 0 0 4 0.4

cliffs 0 0 0 0 0 0 1 2 0.2

vegetation 

greenness CV

0 1 0 0 0 0 0 1 0.1
16

Layer PELO PLGI SCMA THBO UMSC XAVI Prevalence Average 

importance

elevation 0 1 51 5 35 55 11 26.3

winter 

precipitation

49 73 0 61 1 0 8 18.0

surface texture 1 9 0 0 18 2 10 14.6

topographic 

position

36 0 15 14 2 14 10 12.8

seasonal thermal 

change

8 0 11 4 7 5 8 6.3

seasonal 

greenness timing

2 0 4 7 0 1 9 4.6

slope 0 0 0 8 3 11 5 4.5

soil adjusted

greenness

0 0 0 0 0 0 1 2.5

soil water stress 1 5 2 0 6 7 8 2.4

sandy surfaces 0 0 0 0 27 0 2 2.1

mountain bases 1 10 6 0 0 3 7 2.1

washes 1 0 4 0 0 1 8 1.7

alluvial fan 

proximity

0 0 0 0 0 0 2 0.9

rocky surfaces 0 0 6 0 0 0 1 0.4

seasonal 

greenness ratio

2 0 0 2 1 1 4 0.4

cliffs 0 2 0 0 0 0 2 0.2

vegetation 

greenness CV

0 0 0 0 0 0 1 0.1
17

Elevation

Elevation was included in habitat models for all but two of our study taxa, and had the highest 

average importance (26.3 percent; table 4) of the environmental covariates considered here. However, 

the relationship between habitat suitability and elevation is not a causal one in itself, but the result of 

interacting physical and biological factors that covary with elevation. Elevation is related to climate 

(temperature, precipitation, humidity, etc.), soil, substrate type, and the vegetation community of an 

area, and has been shown to be a reliable predictor of species distributions across many taxa (Merriam, 

1890; Franklin, 1995). Elevation is less helpful for identifying ecological relationships between 

organisms and their environment, but because of its explanatory power, we elected to include it. The 

only taxa where elevation was not included in any of the selected models were the Homalonychus spider

group and P. longimembris. Very little information is known about the two spider species, and few 

observations were available for identifying strong covariate relationships. These observations suggested 

that Homalonychus spiders occur between elevations of 300 and 1,400 m. In contrast, P. longimembris

appear to occur over a broad range of elevations from -70 to 2,700 m. This wide range likely obscured 

any explanatory power that elevation may have imparted. 

Winter Precipitation

Winter precipitation represents the cumulative precipitation from October 1 through March 31, 

averaged across 30 years. This covariate was included in habitat models for 8 of the 13 taxa, and 

showed the second highest average importance value (18.0 percent; table 4). Many of the taxa included 

in this study are herbivorous, feeding primarily on herbaceous grasses and forbs. The germination and 

growth of many grasses and forbs in the Mojave Desert are influenced heavily by winter precipitation 

(Beatley, 1976; Turner and Randall, 1989). Therefore, we hypothesized that many of the herbivorous 

species modeled here would be positively associated with winter precipitation. We also expected that 

many of our insectivorous species (such as X. vigilis) might be positively associated with high winter 

precipitation because of the increased abundance of insect forage resulting from spring production of 

annual plants that may drive mass insect production events (Swetnam and Lynch, 1993). The final 

averaged models for all 4 herbivorous taxa and 4 of the 9 facultative and obligate carnivores included 

winter precipitation. P. gilberti, which primarily eats small arthropods, had the highest average relative 

contribution for winter precipitation among the study species. This positive relationship suggests a 

preference for areas with more winter rainfall, although this was unlikely a preference for areas with 

high herbaceous production for foraging, but rather, a preference for areas across a wide elevational 

range (sea level to more than 2,000 m) near areas with high surface moisture, such as near streams or 

springs. 

Surface Texture

Surface texture broadly reflects the dominant surface mechanical composition, and is an 

indication of the effective grain size for soil and rock fragments (Scott Nowicki, University of Nevada, 

Las Vegas, oral commun., 2011). This index is broadly determined by the thermophysical 

characteristics of its surface substrate material, which reflects the material type and composition, 

temperature holding capacity, particle size, and induration of an area. These characteristics are directly 

related to the range of vegetation that can occur in a given area, and as such, have a strong influence on 

the species composition there. In addition, many species appear to prefer specific particle sizes that vary 

among species (for example, lizards; Pianka, 1966). We modeled the size of surface particles to use for 

a proxy for surface texture. Particle size was modeled as a function of the thermal properties of the 
18

substrate surface. In general, areas with smaller particle sizes including little to no bedrock cover show 

greater diurnal thermal variation than rocky areas with extensive bedrock and boulders (Pratt, 1979; 

Kahle, 1987; Wang and others, 2004). Therefore, areas with greater differences between daytime and 

nighttime surface temperatures are more likely to be sandy or have small particle sizes than areas with 

smaller daytime and nighttime summer differences. This property is magnified during the hot summer 

months when the available solar energy is at a maximum. We used this property to create a covariate 

representing surface texture throughout our study area 

The surface texture covariate was included in habitat models for 10 of the 13 taxa, and showed 

the third highest average importance of the environmental covariates (14.6 percent; table 4). In general, 

species that tend to occur on rocky and bouldery surfaces show preferences for areas with a low surface 

texture index, which is indicative of areas that carry greater diurnal thermal inertia. Species that occur in 

areas with sandy or loamy soils tend to show preferences for areas with high surface texture index 

values. 

Topographic Position

Topographic position represents an index of the landscape terrain shape and connectivity of an 

area, and was calculated as the amount of surface water that could drain into a given1 km cell given an 

impermeable surface. High values are indicative of dry lakebeds, valley bottoms, and surface flow pinch 

points, such as the apexes of alluvial fans. Low values represent ridges and mountain tops, or other areas 

where theoretical waterflow would be minimal.

The topographic position covariate showed the fourth highest average importance (12.8 percent; 

table 4), and was included in habitat models for 10 of the 13 taxa. Response curves for topographic 

position varied considerably across taxa, with some taxa showing a preference for specific values (that 

is, C. bicinctores and C. collaris and S. magister), others showing a preference for a broad range of 

values (P. gilberti and T. bottae), and others showed a preference for values above a certain threshold 

(D. dorsalis and X. vigilis), and yet others exhibited a preference for high values (C. occipitalis and C. 

penicillatus). Generally, none of our taxa showed preference for low values, suggesting that most of our 

study species do not occupy mountainous terrain. 

Seasonal Thermal Change

Seasonal thermal change represents the similarity in summer and winter surface temperatures for 

a given area, and is a measure of the intra-year variability in surface temperature. High values, which 

represent areas with greater differences between the summer and winter daytime surface temperature, 

tended to occur on valley bottoms in the northern latitudes of the study area. In contrast, places with low 

values were high elevation, south-facing slopes, where summer and winter daytime temperatures were 

more similar. 

The seasonal thermal change covariate was included in habitat models for 8 of the 13 taxa, and 

showed the fifth highest average importance of the environmental covariates (6.3 percent; table 4). Most 

taxa showed preferences for areas with low seasonal temperature differences, suggesting an avoidance 

of the high latitudes. However, three of the eight species (S. magister, U. scoparia, and X. vigilis) 

showed unique response curves for seasonal thermal change, with preferences for moderate values, and 

avoidance of the two extremes.
19

Seasonal Greenness Timing

Seasonal greenness timing represents the period of the year when vegetation greenness reaches 

its maximum. We used the date where the maximum value of the Normalized Difference Vegetation 

Index (NDVI; vegetation greenness) was observed to represent the date of peak vegetation greenness. 

Vegetation green-up as measured by NDVI has been used to distinguish multiple vegetation 

communities, including shrubland, grassland, and mixed shrub communities (Peters and others, 1997; 

Weiss and others, 2004), and the timing of greenness has been shown to correlate with elevation 

gradients (Crimmins and others, 2009) and vegetation communities (van Leeuwen and others, 2010). 

Riparian areas and pinyon-juniper communities tended to have later dates of peak vegetation greenness. 

Prior to use, this covariate was transformed to a nominal index of 1–9, representing the dates identified 

above. 

The seasonal greenness timing covariate was included in habitat models for 8 of the 13 taxa, and 

showed the sixth highest average importance of the environmental covariates (4.6 percent; table 4). 

Most species had increased habitat suitability in areas where vegetation greenness peaked later in the 

season, although C. occipitalis and the Homalonychus spiders showed preferences for areas where 

vegetation greenness peaked earlier in the season. 

Slope

Slope was included in habitat models for 5 of the 13 taxa, and was the seventh highest in terms 

of average importance (4.5 percent; table 4). Of the five habitat models where slope was included, three

showed relative contributions less than 10 percent, suggesting that slope provides relatively little 

explanatory power for most of our study species. The notable exceptions were the Crotaphytus lizards 

(C. bicinctores and C. collaris) and X. vigilis, which had relative contributions greater than 10 percent. 

Although Crotaphytus lizards prefer higher slopes and areas with greater rock abundance, desert night 

lizards prefer flat areas with little surface relief. 

Soil Adjusted Greenness

Soil adjusted greenness represents vegetation greenness that has been adjusted to minimize the 

influence of soil color and spectral properties (for example, albedo) on vegetation greenness 

measurements (Qi and others, 1994). Soil adjusted greenness had the eighth highest average importance 

of the covariates (2.5 percent; table 4), and was included in habitat models for only 1 of the 13 taxa. 

Although not used in models for any of the other taxa, soil adjusted greenness provided the second 

highest relative contribution for A. punctatus, and the modeled response suggested a prevalence in areas 

with moderate to high vegetation greenness. 

Soil Water Stress

Soil water stress represents an index of the in situ top layer soil moisture (Fensholt and Sandholt, 

2003). Soil water stress was included in habitat models for 8 of the 13 taxa, but only showed an average 

importance of 2.4 percent (table 4), with the highest relative contribution for X. vigilis (7 percent, table 

4), which preferred areas with high soil water moisture. In contrast, the Crotaphytus lizards preferred 

areas with low soil water moisture.
20

Sandy Surfaces

Sandy surfaces represents the percentage of each 1 km cell that was categorized as sand. We 

categorized sand using three separate modeling steps: (1) MODIS MOD11A1 nighttime Land Surface 

Temperature (LST) 8-day data were downscaled from the original spatial resolution of 1 km to 250 m 

using a geographically weighted generalized additive model with seven environmental predictor 

covariates (solar beam irradiance, average of the black sky albedo for the same time periods as the LST, 

maximum horizon level, elevation, aspect, slope and an average of the NDVI); (2) the downscaled 

MODIS LST data were regressed against elevation to remove the elevation trend; and (3) the residuals 

from the elevation regression were categorized into sand and non-sand categories using 100 known 

sample locations in Nevada and California. Sandy surfaces was only represented in habitat models for 2 

of the 13 taxa, and had an average importance of 2.1 percent (table 4) among all species, and a relative 

contribution of less than 1 percent in D. dorsalis, but 27 percent in U. scoparia (table 4). U. scoparia

was the only aeolian obligate in our study species and showed a preference for areas with higher 

coverage of sandy surfaces.

Mountain Bases

Mountain bases represents the transition areas between steep mountain slopes and flatter valleys, 

and include unique combinations of precipitation, temperature, and soil types. Our environmental 

covariate mountain bases represented the percentage of each 1 km cell that was categorized as a 

mountain base. We considered mountain bases to be areas that (1) have less than 10 degrees slope; (2) 

have a slightly negative profile curvature; and (3) had surface roughness between 0.7 × 10-8 and 5.8 × 

10-8. The mountain bases covariate was included in habitat models for 7 of the 13 taxa, and had an 

average importance of 2.1 percent, although the highest relative contribution was for P. gilberti (10

percent; table 4), which showed a preference for areas categorized as the bases of mountains. L. 

trivirgata also showed a preference for areas categorized as the bases of mountains, and had a relative 

contribution from mountain bases of 5 percent (table 4). 

Washes

The washes layer represents the relative density of washes per square kilometer in alluvial fan 

areas. The washes covariate had the twelfth highest average importance (1.7 percent; table 4), and was 

included in habitat models for 8 of the 13 taxa. In the Homalonychus spiders habitat model, washes 

provided a 15 percent relative contribution, although in the remaining 5 taxa, washes provided less than 

4 percent relative contribution (table 4). 

Alluvial Fan Proximity

Alluvial fans are fan shaped deposits of soil and sediments that have been carried by surface 

water from upslope terrain. These deposits tend to form gradients in particle size from coarse to fine 

when moving down slope. This environmental covariate represents the distance in meters to the nearest 

cell categorized as an alluvial fan. The alluvial fan proximity covariate was included in habitat models 

for only 2 of the 13 taxa, and had an average importance of 0.9 percent. However, in the habitat models 

for L. trivirgata, alluvial fan proximity had a relative contribution of 8 percent, which was the fourth 

highest contribution for this species. Across all study species, alluvial fan proximity provided lower 

explanatory power than we had hypothesized, likely due to the coarse spatial scale at which we modeled 

habitat. 
21

Rocky Surfaces

Rocky terrain is a fundamental component of habitat for many of our study species, and includes 

varying landform features that offer a wide range of surface conditions including caves, outcrops, 

ridges, peaks, and boulders. The rocky surfaces covariate represents the percentage of each 1 km cell 

that was categorized as rock. We categorized areas as rock using three separate modeling steps: (1) 

MODIS MOD11A1 nighttime LST, 8-day global data were downscaled from the original spatial 

resolution of 1 km to 250 m using a geographically weighted generalized additive model with seven

environmental predictor covariates (solar beam irradiance, average of the black sky albedo for the same 

time periods as the LST, maximum horizon level, elevation, aspect, slope, and an average of the NDVI); 

(2) we regressed the downscaled MODIS LST data against elevation to remove the elevation trend; and 

(3) we categorized residuals from the elevation regression into rock and non-rock categories using 100 

known sample locations in Nevada and California. The rocky surfaces covariate was included in habitat 

models for only 1 of the 13 taxa, had an average importance of less than 0.4 percent, and was among the 

lowest 5 environmental covariates. 

Seasonal Greenness Ratio

The phenology of desert plants has been related to the composition and resource use of the 

organisms in an ecosystem (Esler and Rundel, 1999), and is highly dependent on winter precipitation. In 

a similar effort, Wallace and Thomas (2008) show that the green-up potential, or relative difference 

between spring greenness in a dry year and a wet year, is related to the ground cover of annuals at 

multiple sites in the Mojave Desert. We contrasted spring and winter vegetation greenness in order to 

capture the geographic variation in greenness between two separate seasons (periods when annuals were 

growing and winter when only perennials were present) in the Mojave Desert. The seasonal greenness 

ratio covariate was included in habitat models for 4 of the 13 taxa; however, in each of the 4 taxa, 

relative contributions were less than 2 percent, resulting in an average importance of 0.4 percent 

(table 4). 

Cliffs

Cliff-like landform features were hypothesized to represent habitat for several of our study 

species, such as L. trivirgata and S. magister. The cliffs covariate represented the percentage of each 1 

km cell that was categorized as containing cliffs and was included in habitat models for 2 of the 13 taxa. 

Due to its low average importance (0.2 percent; table 4), the cliffs covariate was not considered an 

important covariate for predicting habitat of these species.

Vegetation Greenness Coefficient of Variation

The variability of vegetation greenness throughout each year was hypothesized to represent the 

availability of forage for many of our study species, and was calculated as the coefficient of variation in 

the NDVI across 23 dates in each year from 2001 to 2010. The vegetation greenness coefficient of 

variation covariate was included in habitat models for only 1 of the 13 taxa, and had the lowest average 

importance (0.1 percent; table 4) of the 17 environmental covariates. 
22

Habitat Model Summaries

Anaxyrus punctatus

Conceptual Model.—A. punctatus occur in a wide range of semi-aquatic habitats in and around 

springs or water catchments (Stebbins, 1951; Ferguson and Lowe, 1969; Creusere and Whitford, 1976; 

McClanahan and others, 1994), in canyons with rocky substrates (Sullivan, 2005), and in rocky upland 

habitats or floodplains with elevations up to 2,000 m (Stebbins, 1951; Ferguson and Lowe, 1969; 

McClanahan and others, 1994; Sullivan, 2005). A. punctatus often stray away from water sources to 

burrow in soil or rock crevices (McClanahan and others, 1994), and are known to be active during 

precipitation events (Tevis, 1966; Ferguson and Lowe, 1969; McClanahan and others, 1994). These 

observations of habitat preference were translated to covariate layers representing rocky surfaces, 

surface texture, topographic position, mountain bases, and washes to account for terrain characteristics. 

We expected to see high suitability habitat in areas with high values of rocky surfaces, high values of 

surface texture, low values of topographic position, high values of mountain bases, and high values of 

washes. 

We evaluated several covariate layers depicting surface-water characteristics for this semiaquatic species including soil water stress, vegetation greenness, soil adjusted vegetation greenness, 

seasonal greenness ratio, vegetation greenness timing, and winter precipitation. We expected our models 

of habitat suitability to have high suitability in areas with low soil water stress (corresponding to high 

water availability), high vegetation greenness, high soil adjusted vegetation greenness, high seasonal 

greenness ratio (corresponding to high spring and high winter vegetation), low to medium seasonal 

greenness timing (corresponding to early and mid-season maximum greenness), and high winter 

precipitation. Although the water-holding capacity of soil has been shown to be a good predictor of A. 

punctatus habitat suitability (Dayton and Fitzgerald, 2006), we were unable to use any direct measures 

of the water-holding capacity of soil. Instead, we used the covariate layer representing soil water stress, 

which has been correlated with in situ top layer soil moisture in large portions of African semiarid 

deserts (Fensholt and Sandholt, 2003). 

A. punctatus have been documented in areas with early onset of spring (warmer) air 

temperatures (Tevis, 1966; Ferguson and Lowe, 1969). We included the covariate layers summer 

temperature and seasonal thermal change to express these habitat characteristics. We hypothesized that 

areas with low summer temperatures and low seasonal thermal variation would be more suitable for A. 

punctatus because A. punctatus are poikilotherms, and must select microhabitats with tolerable climates 

rather than relying on metabolic thermoregulation. Bradford and others (2003) found habitat of A. 

punctatus to be associated with high vegetation canopy cover over water and adjacent land. We 

represented riparian canopy cover with our soil adjusted greenness, which was derived for autumn 2010, 

during the time of year when riparian vegetation is more pronounced from surrounding areas in order to 

improve the discrimination of riparian areas. 
23

Habitat Model Results.—Although a single GAM model had the highest AUC among the 

selected models, we included two additional GAM models because of their high AIC and BIC scores, as 

well as their high BI scores and ideal BI curves. These three selected models also were run as MaxEnt

models, resulting in a total of six models that we reviewed. All six models were included in the final 

averaged model, which represented five covariate layers—elevation, soil water stress, surface texture, 

soil adjusted greenness, and winter precipitation. The most prevalent covariates were surface texture, 

soil adjusted greenness, and elevation, and each were included in all six selected models. The least 

prevalent covariate layers were winter precipitation and soil water stress, which were included in only 

two of the six models. The surface texture layer provided the highest relative contribution in each of the 

MaxEnt models (40 percent), followed by soil adjusted greenness (32 percent), elevation (19 percent), 

winter precipitation (6 percent), and soil water stress (3 percent). 

The final averaged model had an AUC score of 0.83 and a BI score of 0.89 (table 1). There was 

high agreement (that is, a low standard error) across the six selected models, with a few notable 

exceptions in large dry lakebeds and low elevation areas, such as the lower portions of Death Valley 

National Park and along the Colorado River (fig. 2). The averaged A. punctatus habitat suitability model 

was reclassified into high and low habitat suitability, and resulted in 87,510 km2 of suitable habitat 

based on the fifth quantile of the observation records (break value=0.35; fig. 3). 
24

Figure 2. Map showing standard error (σ/N) among the selected models for Anaxyrus punctatus (Red-spotted 

Toad). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
25

Figure 3. Map showing maximum habitat suitability among the selected models for Anaxyrus punctatus (Redspotted Toad). Values range from from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
26

Habitat Model Discussion.—Available data were insufficient to derive an adequate 

environmental covariate representing all possible springs and water catchments in the Mojave Desert, 

which are thought to be the single best indicators of A. punctatus habitat (Stebbins, 1951). As a proxy 

for surface water presence, we investigated several covariate layers depicting surface-water 

characteristics including winter precipitation, soil water stress, vegetation greenness, soil adjusted 

vegetation greenness, seasonal greenness ratio, and vegetation greenness timing. Of these layers, only 

soil adjusted greenness was included in all selected models and had a relative contribution of 32 percent

(table 4). This layer showed that A. punctatus habitat suitability was lower when soil adjusted greenness 

was extremely low, indicative of dry lakebeds and areas with little or sparse vegetation, and higher 

when values were greater than 0 (fig. 4C). This is congruent with observations of A. punctatus, which 

are not known to occupy habitats with dry lakebeds or areas with sparse vegetation. 

The most influential environmental covariate was surface texture, which also was included in all 

selected models and had a relative contribution of 40 percent. The response curves suggested that A. 

punctatus tend to occur in areas with low values of surface texture, indicative of rocky surfaces, and 

tend to avoid areas with high values, which represent sandy and silty surfaces (fig. 4B). Areas with 

values greater than 20 included surfaces with mixed particle sizes, including bedrock, boulders, and 

rocks, but also likely including smaller patches of coarse sandy soils as well. The surface texture layer 

was included in the selected models although our rocky surfaces layer was not. The rocky surfaces layer 

was a representation of areas containing extensive surficial rock (that is, bedrock), and was an indicator 

of the amount of surficial bedrock present in an area. The selection of surface texture over rocky 

surfaces was likely due to the prevalence of A. punctatus in areas adjacent to extensive bedrock and in 

areas with mixed surface particle sizes. Extensive bedrock can provide extra runoff resulting in areas of 

temporary surface-water accumulation that provide A. punctatus with foraging areas while also 

accommodating cover sites in nearby rocky soils for periods of inactivity. In the Mojave Desert, areas 

with low surface texture values (that is, less than 20) include features such as canyons and rock outcrops 

where solar radiation, wind patterns, and surface material have a strong influence on the surface radiant 

temperature (Carlson and others, 1981) and create highly variable environments such as crevices and 

chasms that are more thermally stable throughout each day during the summer months. This has been 

corroborated by field biologists who have found that A. punctatus tend to prefer canyons with rocky 

substrates (Sullivan, 2005), or rocky upland habitats (Stebbins, 1951; Ferguson and Lowe, 1969; 

McClanahan and others, 1994). 
27

The elevation layer also was included in each of the selected models, had a relative contribution 

of 20 percent, and habitat suitability decreased as elevation increased (fig. 4A). However, the 

relationship between habitat suitability and elevation is likely not a causal one in itself, but instead is the 

result of interacting physical and biological factors that co-vary with elevation. Including elevation as an 

environmental covariate often masks the relationship between habitat suitability and other covariates 

because elevation is related to the climate (temperature, precipitation, humidity, etc.), soil, substrate 

type, and vegetation community of an area. However, because of its relationship to these important 

covariates, elevation can often statistically explain a large portion of a species’ geographic distribution 

(Franklin, 1995). We included elevation because our intent was to create accurate predictions of habitat 

suitability across the Mojave Desert Ecoregion for the purpose of evaluating influences of habitat loss 

due to construction of new alternative energy facilities. 

The covariate layers, winter precipitation and soil water stress, provided the least explanatory 

power among the environmental covariates we selected (6 and 3 percent, respectively), and were 

included in only two of the six habitat models. We expected winter precipitation to have a large 

influence on habitat suitability for A. punctatus due to the transient nature of water exchange between 

amphibians and their environment (Tracy, 1976). However, our winter precipitation layer was correlated 

with elevation (r = 0.71, N = 1,000, p < 0.001), and therefore the low prevalence and low relative 

contribution (6 percent) of winter precipitation in the selected models is understandable; when 

combined with elevation, winter precipitation provided very little additional explanatory power in A. 

punctatus habitat. We also expected the soil water stress layer would have a high relative contribution to 

A. punctatus habitat than it did (3 percent) due to A. punctatus’ need for surface water. However, 

although the soil water stress layer has been shown to be related to the in situ top layer soil moisture in 

semiarid regions of Senegal (Fensholt and Sandholt, 2003), it has not been evaluated in the Mojave 

Desert, and may not be an important predictor of A. punctatus habitat.
28

Figure 4. Generalized response curves and histograms for Anaxyrus punctatus (Red-spotted Toad) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for A. punctatus, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0.
29

Chaetodipus penicillatus

Conceptual Model.—C. penicillatus tends to occur in or around areas with mildly sloping 

valleys, low hills with washes, and open sandy deserts, often with coarse sandy soils (Bailey, 1931; 

Hall, 1981; Hoffmeister, 1986). Others have found that several rodent species (primarily the 

Chaetodipus, Dipodomys, and Peromyscus genera) could be distinguished in terms of habitat use with a 

suite of environmental covariates including slope, rockiness, soil texture, vegetation cover, and 

precipitation (Suarez-Gracida and Alverez-Casteneda, 2009). C. penicillatus are often found in sparsely 

vegetated areas with perennial plant species, such as Larrea tridentata (creosote bush), Cylindropuntia 

spp. (cholla), Cercidium spp. (palo verde), and Ambrosia spp. (burroweed) (Hoffmeister, 1986), or on 

the edges of alluvial fans with Yucca spp. (yucca), Prosopis spp. (mesquite), Bouteloua spp. (grama 

grass), and Argemone spp. (prickly poppy) (Hoffmeister, 1986). We used vegetation greenness and soil 

adjusted greenness to capture the geographic variation in perennial vegetation sparseness, and 

vegetation greenness coefficient of variation and seasonal greenness ratio to capture the yearly and 

seasonal persistence of green vegetation. An additional vegetation covariate, vegetation greenness 

timing, was used to represent the timing of peak of vegetation greenness. The layer alluvial fan 

proximity was used to represent the proximity to alluvial fan edges, and layers slope, surface texture, 

topographic position, and mountain bases were used to represent terrain and landform features. 

Habitat Model Results.—We selected 9 GAM models to run as MaxEnt models, resulting in a 

total of 18 models that we reviewed. Only two of the MaxEnt models and none of the GAM models 

were included in the final averaged model, which represented five covariate layers. Of these five layers, 

surface texture, seasonal thermal change, and elevation were included in both of the selected models, 

and winter precipitation and mountain bases were included in only one model. The relative contribution 

for each layer was 45 percent (elevation), 23 percent (seasonal thermal change), 22 percent (surface 

texture), 8 percent (winter precipitation), and 2 percent (mountain bases). The two models included in 

the final averaged model were selected because they discriminated habitat in areas where no presence 

points were available, but where we had local knowledge of presence. There was general agreement 

across the two selected models, although we did not calculate a standard error map due to low sample 

size. The final habitat suitability model was the average of the two selected models, and had an AUC 

score of 0.87 and BI of 0.22 (table 1). The averaged model was reclassified into high and low habitat 

suitability, which resulted in 53,988 km2

 of suitable habitat based on the fifth quantile of the observation 

records (break value=0.25). Areas of suitable habitat were limited to southeastern Mojave Desert, and 

throughout much of the Baker Sink landform in Death Valley National Park (fig. 5).
30

Figure 5. Map showing average habitat suitability among the selected models for Chaetodipus penicillatus (Desert 

Pocket Mouse). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
31

Habitat Model Discussion.—We hypothesized that habitat suitability of C. penicillatus would 

be related to broad environmental covariates representing vegetation and landforms due to the range of 

habitat types that C. penicillatus are known to occupy. However, the covariate layers in the two selected 

models for C. penicillatus (surface texture, seasonal thermal change, mountain bases, elevation, and 

winter precipitation) did not include any of the layers describing vegetation conditions, suggesting that 

our vegetation layers were not as important for predicting C. penicillatus habitat suitability as we 

expected. Elevation had the greatest explanatory power (45 percent contribution) in both of the selected 

models, and the model response suggested an avoidance of high elevation areas (fig. 6A). Elevation was 

highly correlated with winter surface temperature (r = -0.95, N = 1,000, p < 0.001), (r = -0.95, N = 

1,000, p < 0.001), and therefore we did not include both in any of our models, but instead used one or 

the other. Models using elevation instead of winter surface temperature generally performed better than

the converse. However, examination of models with winter surface temperature indicate that C. 

penicillatus has high habitat suitability in areas with winter temperatures ranging between 12 and 20 °C. 

This range includes temperatures known to induce torpor in this species (10–15 °C; Bartholomew and 

Cade, 1957). 

Our layer representing seasonal change in surface temperature provided the second highest 

relative contribution to C. penicillatus habitat (23 percent), and the modeled response suggested that 

areas with less seasonal change were more suitable than areas with greater seasonal change (fig. 6B). 

Throughout our study area, areas with low seasonal thermal change tended to be areas with very mild 

winters and moderate summers, in contrast to areas with cold winters, which provided the greatest 

seasonal thermal change. The later conditions occurred primarily in the northern latitudes in areas with 

high solar exposure, such as high elevation valley bottoms. 

We hypothesized that surface texture would be an important factor in C. penicillatus habitat due 

to field observations of fossorial activities in coarse sandy soils throughout its range (Bailey, 1931; 

Findley and others, 1975; Hoffmeister 1986). Our models supported this, with surface texture being 

included in each of the two selected models and with 22 percent relative contribution overall. In each 

case, the surface texture layer showed decreased habitat suitability in areas with surface bedrock or 

other rocky surfaces, and an increase in habitat suitability in areas with sandier surfaces (fig. 6C). This 

preference for areas without large rocks or large areas of exposed bedrock also was supported by the 

mountain bases layer, and although it only provided 2 percent relative contribution, with decreased 

habitat suitability in areas characterized by mountain bases (fig. 6D). Combined, the surface texture and 

mountain bases layers provide support for previous assertions that C. penicillatus prefer flat to mildly 

sloping valleys with open, coarse, and sandy soils (Bailey, 1931; Hall, 1981; Hoffmeister, 1986). 

In terms of diet, C. penicillatus depends heavily on spring annual plants for food during the 

breeding season, which is closely tied to the timing and amount of winter precipitation (Beatley, 1976; 

Turner and Randall, 1989). We found a positive correlation between the winter precipitation layer and 

C. penicillatus habitat suitability, with a peak between 76 and 101 mm of winter precipitation (fig. 6E). 

This peak suggests that although more winter precipitation can be better, areas with winter precipitation 

greater than 127 mm, such as mountain ranges, provide lower habitat suitability than valleys. 
32

Figure 6. Generalized response curves and histograms for Chaetodipus penicillatus (Desert Pocket Mouse) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for C. penicillatus, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0. Only the five highest contributing covariates are 

shown.
33

Chionactis occipitalis

Conceptual Model.—C. occipitalis typically occur in dry, sandy areas that are relatively flat and 

low in elevation, such as valley bottoms and alluvial fans (Klauber, 1951; Funk, 1967; Stebbins, 2003). 

Phil Medica (U.S. Geological Survey, oral commun., 2011) stated that sand dunes may be suitable if 

sufficient vegetation is present, as they are not often found in sand dunes or areas with extensive sand 

and void of vegetation altogether (Klauber, 1951). We included the environmental covariate layers

surface texture and sandy surfaces to capture the preference for sandy areas and sand dunes, and the 

layers alluvial fan proximity, topographic position, washes, and mountain bases to represent valley

bottoms and other relatively flat, low elevation areas. 

Like all ectotherms, C. occipitalis behavior is affected by the temperature of its immediate 

environment, and requires specific temperatures for activity (Klauber, 1951). Throughout most of its 

range, C. occipitalis is generally nocturnal and is active when air temperatures are between 70o and 90oF 

(Rosen and others, 1996). We included summer surface temperature to represent the geographic 

variation in summer environmental temperatures. To capture the avoidance of areas with little or no 

vegetation, we used vegetation greenness coefficient of variation and seasonal greenness ratio to capture 

the yearly persistence of green vegetation. We also included vegetation greenness timing, soil water 

stress, and winter precipitation to represent the availability of water.

Habitat Model Results.—We selected 6 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 12 models that we reviewed. All 12 models were 

included in the final averaged model, which represented six covariate layers—surface texture, 

topographic position, elevation, soil water stress, seasonal greenness timing, and vegetation greenness 

coefficient of variation. The three most important environmental covariates were surface texture, 

elevation, and topographic position, with a prevalence in 12, 10, and 10 total models and relative 

contributions of 42, 27, and 22 percent, respectively. The layer seasonal greenness timing was present in 

6 of the 12 selected models, and provided a relative contribution of 6 percent, although vegetation 

greenness was present in only 2 of the 12 models, and provided a relative contribution of only 1 percent.

Both GAM and MaxEnt models were included in the final candidate suite because there were 

few differences among the paired models. Across all models, there was moderate consistency in habitat 

suitability scores throughout the entire study area, with slight disagreement among models in some of 

the low elevation areas, such as Death Valley National Park and outside of Phelan, Calif. (fig. 7). The 

final habitat model was the average of the 12 selected models, which had an AUC score of 0.79 and a BI 

score of 0.64 (table 1). The averaged C. occipitalis habitat suitability model was reclassified into high 

and low habitat, and resulted in 85,397 km2

 of suitable habitat based on the fifth quantile of the 

observation records (break value=0.26; fig. 8). 
34

Figure 7. Map showing standard error (σ/N) among the selected models for Chionactis occipitalis (Shovel-nosed 

Snake). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
35

Figure 8. Map showing average habitat suitability among the selected models for Chionactis occipitalis (Shovelnosed Snake). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
36

Habitat Model Discussion.—We hypothesized that C. occipitalis would be associated with 

covariate layers representing sandy surfaces and flat areas due to their preference for alluvial fans and 

valley bottoms (Klauber, 1951; Funk, 1967; Stebbins, 2003). As expected, surface texture and 

topographic position were among the best predictors of C. occipitalis habitat, with C. occipitalis habitat 

suitability being positively related to fine-textured soils, such as sand, and negatively related to more 

rocky surfaces, such as bedrock and areas with large boulders (fig. 9A). Topographic position showed 

high habitat suitability occurring in areas with high values, such as alluvial fans and valley bottoms (fig. 

9B). Combined, these two covariates suggest that C. occipitalis were associated with sandy valley 

bottoms, alluvial fans, and playa margins (to the exclusion of the playas themselves that are void of 

vegetation). The elevation layer showed a decreasing relationship with C. occipitalis habitat suitability 

at elevations greater than 1,000 m, and a peak in suitability between 400 and 1,000 m (fig. 9C). Above 

1,000 m and below 400 m elevation, habitat suitability for C. occipitalis was minimal. Areas of high 

habitat suitability also were associated with high soil water stress (fig. 9D), indicating a preference for 

areas with more surface water.

Seasonal greenness timing was negatively associated with habitat suitability, suggesting that 

areas where vegetation reaches peak greenness earlier have higher suitability than areas where peak 

vegetation greenness occurs later in the year (fig. 9F). These areas tended to have Mojave mixed scrub 

with a low density of perennial vegetation, and a high density of grasses and other annual vegetation 

with earlier growth seasons. The other vegetation layer represented in the final averaged model was 

vegetation greenness coefficient of variation. High variability in vegetation greenness was negatively 

related to C. occipitalis habitat suitability (fig. 9E), suggesting a preference for areas with more stable 

vegetation. 
37

Figure 9. Generalized response curves and histograms for Chionactis occipitalis (Shovel-nosed Snake) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for C. occipitalis, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
38

Crotaphytus bicinctores and C. collaris

Conceptual Model.— C. bicinctores and C. collaris (collared lizards) are associated with a 

wide range of habitats, but are primarily found in areas that contain boulders, rocky surfaces, and other 

rock-like formations (McGuire, 1996; Stebbins, 2003). We included environmental covariates of surface 

texture, rocky surfaces, mountain bases, topographic position, washes, surface roughness, and slope to 

capture areas with rocky habitat. In many parts of the Mojave, C. bicinctores are found in areas with 

sparse vegetation at elevations less than 3,000 m (McGuire, 1996; Stebbins, 2003), and do not often 

occur in areas with riparian vegetation, alluvial fans, or flat valley bottoms (McGuire, 1996; Stebbins, 

2003). To capture the association with sparsely vegetated areas, we included vegetation greenness, and 

vegetation greenness coefficient of variation and seasonal greenness ratio were used to capture the 

yearly persistence of green vegetation. We also included vegetation greenness timing and soil water 

stress to represent areas without riparian vegetation. 

Habitat Model Results.—We selected five models from among the suite of possible GAM 

models to run as MaxEnt models, resulting in a total of 10 models that were reviewed by our 

collaborators. All 10 models were included in the final averaged model, and in combination, included 

eight covariate layers—surface texture, slope, topographic position, elevation, vegetation greenness 

timing, rocky surfaces, soil water stress, and alluvial fan proximity. The environmental covariate slope

was included in all selected models, and had a relative contribution of 35 percent. Similarly, topographic 

position was included in all selected models and had a contribution of 29 percent, and seasonal 

greenness timing was included in 6 of the 10 models, and had a relative contribution of 13 percent. 

Elevation was included in 4 of the 10 models, and had a relative contribution of 10 percent. The 

remaining four covariates were each included in only 2 of the 10 models, and had contributions of 5 

percent (soil water stress), 4 percent (surface texture), 4 percent (alluvial fan proximity), and 1 percent

(rocky surfaces). 

The final averaged habitat suitability model had an AUC score of 0.77 and a BI score of 0.58 

(table 1). There was substantial disagreement among the selected models across the entire study area 

(fig. 10), with most of the disagreement occurring primarily in areas of high elevation. The averaged 

collared lizard habitat suitability model was reclassified into high and low habitat, and resulted in 

100,293 km2

 of suitable habitat based on the fifth quantile of the observation records (break value=0. 

44; fig. 11).
39

Figure 10. Map showing standard error (σ/N) among the selected models for the Crotaphytus spp. group (collared 

lizards). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
40

Figure 11. Map showing average habitat suitability among the selected models for the Crotaphytus spp. group

(collared lizards). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
41

Habitat Model Discussion.—We hypothesized that C. bicinctores and C. collaris would be 

highly associated with environmental covariates depicting boulders, rocky surfaces, and bedrock, but 

instead we found that the layers slope and topographic position provided the greatest explanatory power 

among the environmental covariates we considered. Collared lizard habitat suitability was positively 

related to slope such that areas with slope greater than 5 degrees showed high habitat suitability (fig.

12A), suggesting an avoidance of flatter areas. The modeled response for the topographic position layer 

suggested an avoidance of low (mountain tops and ridges) and high (high drainage areas, such as dry 

lakebeds and valley bottoms) values (fig. 12B), and a preference for moderate values (transition zones 

between mountain bases and alluvial fans). Although these two covariates describe broad topographic 

and terrain features, they are likely too coarse to identify the rocky features (that is, boulders, small 

cliffs, and outcrops) that collared lizards use as habitat at an individual level. Areas with rocks and 

boulders tend to occur in association with moderate slope or in areas that have high topographic position 

values (that is, alluvial fans, but not dry lakebeds), and as such, collared lizard habitat suitability was 

better near valley edges and at the apexes of alluvial fans. 

Of the layers describing vegetation, seasonal greenness timing had the highest contribution (13 

percent) and showed an increase in habitat suitability in areas where peak vegetation greenness occurred 

later in the season (fig. 12C). These areas tended to be low to mid-elevation mixed desert scrub and 

southwestern warm temperate forests, although the seasonal greenness timing layer did not directly 

model vegetation community types. Instead, it was a proxy for the interaction of vegetation phenology 

and vegetation greenness (van Leeuwen and others, 2010), and was only broadly related to vegetation 

community type. Generally, collared lizards are associated with sparsely vegetated areas, although C. 

collaris are often found in desert scrub, sagebrush, and semi-desert grasslands (Jones and Lovich, 2009).

The response curves for the elevation layer suggested that habitat suitability for collared lizards

was higher at low elevations (fig. 12D), which was supported by a lack of observations in our dataset at 

high elevation areas. However, collared lizards also are not present in the valley bottoms of Death 

Valley and other areas with elevations less than 250 m, suggesting that the linear trend with elevation 

was over-generalized. In several of the MaxEnt models, Crotaphytus habitat suitability was correlated 

with areas between 250 and 1,500 m, suggesting an avoidance of the lowest elevations in the study area. 

The layers soil water stress, surface texture, alluvial fan proximity, and rocky surfaces all had relative 

contributions of less than 5 percent, and as such were not as important in predicting habitat suitability as 

were slope, topographic position, seasonal greenness timing, and elevation. We expected the rocky 

surfaces and surface texture layers to provide greater relative contributions to habitat suitability than we 

observed, although neither of these were able to capture the boulder fields, isolated rock outcrops, or 

other small rocky features at the fine spatial scale needed to represent Crotaphytus habitat. Future 

efforts to model surface texture at finer spatial scales may improve habitat model performance for 

collared lizards, which had some of the poorest performing models (AUC = 0.77, averaged model) 

among the 13 groups we analyzed. 
42

Figure 12. Generalized response curves and histograms for the Crotaphytus spp. group (collared lizards) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for the Crotaphytus spp., where each of the covariates is held 

at their study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing 

covariates are shown.
43

Dipsosaurus dorsalis

Conceptual Model.—D. dorsalis occur in lowland habitats (for example, usually at elevations 

less than 1,200 m) that are open, flat, and dominated by sand or clay (Norris, 1953; Krekorian, 1977; 

Stebbins, 2003). Frequently incised by washes and arroyos, the intervening lower bajadas have only 

minimal amounts of scattered rocks (Norris, 1953; Stebbins, 2003), and allow D. dorsalis to dig 

burrows and underground retreats, which they use as shelter from extreme temperatures and predation 

(Norris, 1953). We included the environmental covariates of surface texture, sandy surfaces, washes, 

seasonal thermal change, and alluvial fan proximity to capture these habitat characteristics. Larrea 

tridentata (Creosote bush) frequently dominates desert iguana habitats (Norris, 1953; Stebbins, 2003), 

where they use the creosote bush as food when in flower and as shelter during periods of excessive heat 

in the spring and summer months. We used our environmental covariates vegetation greenness, seasonal 

greenness ratio, and seasonal greenness timing to capture habitat preferences for areas with open desert 

scrub. D. dorsalis are not usually found in canyons or steep rocky terrain (Norris, 1953). We therefore 

used elevation, slope, and topographic position to identify areas that do not have high elevations, 

canyons, or steep terrain. 

Habitat Model Results.—We selected nine models from among the suite of possible GAM 

models to derive as MaxEnt models, resulting in a total of 18 models that we reviewed. Fourteen models 

were included in the final averaged model, which represented eight covariate layers—elevation, surface 

texture, winter precipitation, topographic position, seasonal greenness timing, seasonal thermal change, 

sandy surfaces, and washes. The environmental covariate elevation was included in all 14 models and 

showed a relative contribution of 70 percent. Surface texture (18 percent) and winter precipitation (14

percent) both showed high relative contributions, and were included in 12 and 8 of the 14 selected 

models, respectively. The remaining five covariate layers, topographic position, seasonal greenness 

timing, seasonal thermal change, sandy surfaces, and washes, showed relative contributions of 3, 2, 2, 1, 

and <1 percent, respectively. The covariate layers, topographic position, washes, and sandy surfaces, 

were each included in only 2 of the 14 models, and seasonal thermal change was included in 6 of the 14 

models, and seasonal greenness timing was included in 8 of the 14 models. 

There was high agreement across all selected models (fig. 13). The final averaged habitat 

suitability model for D. dorsalis had an AUC score of 0.82 and a BI score of 0.59 (table 1), and was 

reclassified into high and low suitability habitat, resulting in 95,914 km2

 of suitable habitat based on the 

fifth quantile of the observation records (break value=0. 21; fig. 14).
44

Figure 13. Map showing standard error (σ/N) among the selected models for Dipsosaurus dorsalis (Desert 

Iguana). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
45

Figure 14. Map showing average habitat suitability among the selected models for Dipsosaurus dorsalis (Desert 

Iguana). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability 

are shown in red; areas with low suitability are shown in blue.
46

Habitat Model Discussion.—We hypothesized that D. dorsalis habitat suitability would be 

related to habitat types including sandy washes, flat terrain, and alluvial fans. Our environmental 

covariates surface texture and topographic position supported this, and modeled responses suggested a 

preference for sandy surfaces (fig. 15B), alluvial fans, and other surface features with high potential for 

drainage (fig. 15D). Slope was not included among the selected models, although elevation was, and its 

response curve suggested a preference for areas at low elevation (fig. 15A). These tend to have low 

slope in the Mojave Desert, and are generally sandier and contain fewer rocks and boulders at their 

surface. The environmental covariate topographic position, although capturing surface features, such as 

alluvial fans, was only included in 2 of the selected models and had a small relative contribution of only 

3 percent. Winter precipitation showed a unimodal relationship to D. dorsalis habitat suitability (fig.

15C), suggesting a preference for areas with extremely low and extremely high winter precipitation. 

However, this likely is due to the broad region-wide relationship between elevation and precipitation, 

wherein high elevations receive more precipitation than low elevations, and not a preference for low or 

high precipitation in itself. Instead, this relationship is likely the result of a preference for a broad range 

of habitats sandy soils with creosote bush (Jones and Lovich, 2009), which occur throughout the Mojave 

Desert. 

Only one of the covariate layers representing vegetation characteristics, seasonal greenness 

timing, was included in the final averaged model and was present in 8 of the 14 selected models. 

Seasonal greenness timing showed a slight negative association with D. dorsalis habitat suitability in the 

MaxEnt models but a slight positive association in the GAM models, suggesting a complex relationship 

with D. dorsalis habitat suitability. Although inconsistent in its relationship to D. dorsalis habitat 

suitability across the modeling platforms, seasonal greenness timing showed similar relative 

contributions among the two platforms, suggesting little overall influence on D. dorsalis habitat 

suitability. We did not explore the relationship between seasonal greenness timing and D. dorsalis

habitat suitability further, but do suggest future work at finer spatial scales to identify and quantify the 

phenology of Mojave Desert vegetation communities. 

The environmental covariate, seasonal thermal change, showed a negative relationship to D. 

dorsalis habitat suitability, suggesting an avoidance of areas with greater seasonal fluctuation in surface 

temperature (fig. 15E), which tended to occur in valley bottoms of the more northern latitudes in the 

study area. D. dorsalis are often found in valley bottoms throughout their range, but are not often found 

at higher latitudes where seasonal fluctuations in surface temperature are high. The environmental 

covariate, sandy surfaces, showed a low contribution to D. dorsalis habitat suitability, and showed a 

negative relationship with D. dorsalis habitat suitability (fig. 15F). This is likely due to the types of sand 

that were modeled, wherein the sandy surfaces covariate layer was primarily aimed at identifying areas 

with aolian sand devoid of vegetation, such as sand dunes. D. dorsalis are not often found on sand dunes 

and, as such, this negative relationship was not unexpected. The relationship between washes and D. 

dorsalis habitat suitability also was negative, suggesting an avoidance of areas with a high density of 

washes. Although this relationship was unexpected, the low relative contribution of washes (<0.4 

percent) and its presence in only two of the selected models suggests that this relationship is marginal, 

and provided little effect on predicting habitat suitability.
47

Figure 15. Generalized response curves and histograms for Dipsosaurus dorsalis (Desert Iguana) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for D. dorsalis, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
48

Homalonychus selenopoides and H. theologus

Summary.—The literature on Homalonychus spiders is sparse, at best. These wandering, nonweb-building spiders inhabit arid environments of the North American Southwest and are found in areas 

with fine sand or other sandy soils with extensive boulders and rocks, which are used as refugia (Roth, 

1984; Domínguez and Jiménez, 2005; Crews, 2009). They are often found under rocks or other debris, 

on hillsides, or in washes (Crews, 2009). In addition to limited information on habitat preference, 

presence data for these species were limited to locality records provided in Crews and Hedin (2006). We 

therefore created a single MaxEnt model with 11 covariate layers that we hypothesized to be important 

to Homalonychus spiders. These layers represented surface texture, seasonal thermal change, elevation, 

slope, topographic position, and washes to represent rocky areas with mixed sand, and to capture any 

geographic variation in terrain that our locality records exhibited. We also included measures of 

vegetation to capture the hypothetically broad range of vegetation types that Homalonychus spiders may 

occupy, because very little is known about any associated vegetation communities. These measures 

included vegetation greenness timing, vegetation greenness, seasonal greenness ratio, and vegetation 

greenness coefficient of variation. 

Covariates contributing less than 10 percent to model explanation were removed, and a reduced 

model was fit with surface texture, topographic position, vegetation greenness timing, and washes, 

resulting in relative contributions of 33 percent (surface texture), 29 percent (topographic position), 23

percent (seasonal greenness timing), and 15 percent (washes). Generally, Homalonychus spider habitat 

suitability was positively related to rocky and moderately sandy terrain (fig. 16A), and Homalonychus

spiders’ localities were not found on very silty areas, such as dry lakebeds and valley bottoms. This also 

was supported by the response curves of the topographic position layer, which also suggested their 

absence on dry lakebeds, valley bottoms (fig. 16B), and a positive association with alluvial fans, valley 

sides, and mild hillslope. Washes were negatively related to Homalonychus spider habitat suitability 

(fig. 16D), such that areas with low wash density generally resulted in higher habitat suitability, 

although the washes covariate layer provided the least contribution of the four covariate layers. The only 

layer representing vegetation characteristic that was related to Homalonychus spider habitat suitability 

was the seasonal greenness timing layer, with response curves suggesting a preference for areas where 

vegetation reached maximum greenness earlier in the growing season (fig. 16C). Areas with earlier 

onset of maximum greenness tended to be areas at low elevations having sparse perennial vegetation, 

and had high densities of annual grasses and forbs. These areas likely have greater proportions of 

moderately sandy textures–another key component to Homalonychus spider habitat. 

The reduced model had an AUC score of 0.82 and BI of 0.21 (table 1), and was reclassified into 

suitable habitat using the fifth quantile of the observation records (break value = 0.21), resulting in 

62,058 km2 (fig. 17). Our use of a single model precluded estimation of standard error in habitat 

suitability for Homalonychus spiders. 
49

Figure 16. Generalized response curves and histograms for Homalonychus genus (Homalonychus spiders) for 

environmental covariances. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for the Homalonychus genus, where each of the covariates is 

held at their study area mean. The dotted line shows linear predictor value of 0. Only the four highest contributing 

covariates are shown.
50

Figure 17. Map showing average habitat suitability among the selected models for Homalonychus spiders. Values 

range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability are shown in red; 

areas with low suitability are shown in blue.
51

Lichanura trivirgata

Conceptual Model.—L. trivirgata occur in areas with granite boulders and chaparral vegetation 

associations (Klauber, 1931). Other habitat associations include non-granite rocky substrates, and areas 

near permanent water sources, such as springs or oases, especially in arid environments (Stebbins, 

2003). We characterized surface material and soil texture by including surface texture, rocky surfaces, 

seasonal thermal change, mountain bases, topographic position, and cliffs. The surface texture and 

rocky surfaces layers were used to identify areas with suitable surface rock material, although seasonal 

thermal change was used to capture areas that were more stable throughout the year with respect to 

surface temperature. We had hypothesized that areas with less thermal stability would be less suitable 

for L. trivirgata because L. trivirgata are poikilotherms, and must select microhabitats with tolerable 

climates rather than relying on metabolic thermoregulation. We used mountain bases, cliffs, and 

topographic position to represent the preference for terrain features, such as low mountain edges, rocky 

canyons, and rocky outcrops. However, L. trivirgata also have been found in rock-free flats with slopes 

less than 5 percent in areas where rodent burrows are common in the Sonoran Desert (Arizona Game 

and Fish Department, 2003). We used alluvial fan proximity, slope, winter precipitation, washes, soil 

water stress, and seasonal greenness ratio to approximate these habitat characteristics. 

Habitat Model Results.—We selected nine models from the suite of possible GAM models to 

run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. 

Twelve of the 18 models were included in the final averaged model, which represented nine covariate

layers—elevation, seasonal thermal variance, winter precipitation, topographic position, vegetation 

greenness timing, washes, cliffs, mountain bases, and alluvial fan proximity. Elevation was in all but 1 

of 12 of the selected models and had the highest relative contribution of 33 percent; seasonal thermal 

change was in 7 of the 12 selected models and had a relative contribution of 22 percent. Winter 

precipitation was in 6 of the 12 selected models and also had a relative contribution of 22 percent. The 

remaining covariate layers were each present in fewer than three models, and each had relative 

contributions less than 10 percent. The covariate layer, alluvial fan proximity, was used in two of the 

selected models and had a relative contribution of 8 percent, and mountain bases was included in 3 of 

the 12 selected models and had a relative contribution of 5 percent. Topographic position also was used 

in 3 of the 12 models and had a relative contribution of 4 percent. The remaining four covariate layers, 

seasonal greenness timing, cliffs, slope, and washes, were included in one model each, with relative 

contributions of 3, 1 , <1, and <1 percent, respectively. 

We found moderate to high agreement among the 12 models in the southern parts of the study 

area, with less agreement in the northern and eastern parts of the study area (fig. 18). In particular, midelevations provided the least agreement among models, while low elevation valley bottoms and high 

elevation mountain peaks had high agreement. The final averaged habitat suitability model for L. 

trivirgata had an AUC score of 0.91 and a BI score of 0.87 (table 1), which was one of the highest 

among our study species. The averaged model was reclassified into high and low habitat suitability, 

which resulted in 59,521 km2

 of high habitat suitability based on the fifth quantile of the observation 

records (break value = 0.34, fig. 19). In general, areas with high habitat suitability were confined to 

mid-elevation, rocky terrain near the base of steep mountains. The Baker Sink landform provides a 

notable break in habitat between the southwestern and northeastern parts of the Mojave Desert Basin.
52

Figure 18. Map showing standard error (σ/N) among the selected models for Lichanura trivirgata (Rosy Boa). 

Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas with low 

similarity are shown in red. Although not a measure of true model error, these agreement maps provide a general 

measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
53

Figure 19. Average habitat suitability among the selected models for Lichanura trivirgata (Rosy Boa). Values range 

from 0~1, and correspond to an index of habitat suitability, where areas with high suitability are shown in red; areas 

with low suitability are shown in blue.
54

Habitat Model Discussion.—We hypothesized that habitat suitability for L. trivirgata would be 

associated with our modeled environmental covariates representing surface substrates (surface texture 

and rocky surfaces), but instead found that elevation provided the greatest relative contribution and was 

the single most prevalent covariate in the selected models. L. trivirgata habitat suitability was negatively 

related to elevation, suggesting a lower suitability of high elevation areas (fig. 20A). As with many of 

the other species, we assume that the importance of elevation is not due to a direct explanation of L. 

trivirgata habitat, but rather a coincidence of the correlation between elevation and other habitat 

covariates that we were unable to measure or model. Surface texture and rocky surfaces covariate layers 

were derived from MODIS Land Surface Temperature data, at a spatial scale of 1 km, which may be too 

coarse to identify rocky features (that is, boulders, small cliffs, and outcrops) that L. trivirgata use as 

habitat. L. trivirgata also tend to occur in areas with relatively stable surface temperatures between 

summer and winter months (fig. 20B). L. trivirgata can inhabit ridge-like features and mountain slopes, 

both of which have more similar summer and winter temperatures than valley bottoms. Latitude, 

however, is the primary determinant of summer and winter surface temperature differences, and in 

general, L. trivirgata are found in the more southern regions of the study area where winter temperatures 

are less extreme. Winter precipitation showed a positive trend with habitat suitability and exhibited an 

apparent threshold near 64 mm of winter precipitation, above which habitat suitability remained high 

(fig. 20C). Alluvial fans do not appear to constitute suitable habitat (fig. 20D), although distances just 

more than 1 km from them were often suitable. The negative relationship of alluvial fans also was 

corroborated by the response curves of the topographic position layer (fig. 20F), which also suggested 

an avoidance of alluvial fans, along with an avoidance of valley bottoms and dry lakebeds. 

The modeled response of the mountain bases layer suggested a preference for areas classified as 

mountain bases (fig. 20E), which often contain granite boulder fields and other rocky terrain. The lack 

of soil water stress or seasonal greenness ratio in any of the selected models was unexpected due to the 

known association of L. trivirgata with permanent water sources in the Mojave Desert (Stebbins, 2003). 

Although not a true representation of permanent water sources, we had hypothesized that these two 

layers would capture some of the geographic variation in L. trivirgata observations near riparian areas. 

The seasonal greenness timing layer was present in one of the selected models, suggesting a preference 

for areas where vegetation achieves peak greenness earlier in the season in this instance. The model 

response to the cliffs layer that suggested a preference for areas without cliffs was unexpected, and the 

model response for slope suggested a preference for areas with high slope. These two layers were 

correlated with one another (r=0.83, N=1,000, p<0.001), but were not included together in any single 

model, and showed very different relationships with L. trivirgata habitat suitability. 
55

Figure 20. Generalized response curves and histograms for Lichanura trivirgata (Rosy Boa) for environmental 

covariates. Histograms are derived from a random sample of the study area and illustrate the distribution of 

covariates across the landscape. Response curves are derived from a single model that is inclusive of all 

covariates represented in the selected models for L. trivirgata, where each of the covariates is held at their study 

area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
56

Perognathus longimembris

Conceptual Model.—P. longimembris occur in areas with relatively flat, sandy, or gravelly 

substrates with widely spaced shrubs (Jameson and Peeters, 1988; Zeveloff, 1988). In some parts of 

their range, P. longimembris are found in gravel washes and in areas with L. tridentata and Yucca spp.

(Chew and Butterworth, 1964; Miller and Stebbins, 1964), and are rarely found on rocky or steep 

slopes. We included washes, surface texture, sandy surfaces, and slope to represent these gravelly and 

wash-like habitat characteristics. Due to their small size, P. longimembris are more susceptible to cold 

weather extremes, and can experience rapid changes in body temperatures if exposed to pronounced 

environmental variation (Bartholomew and Cade, 1957). This may limit their occupied habitat to areas 

without cold extremes during winter months and areas that provide longer durations of acceptable 

surface temperatures. We included winter surface temperature and seasonal thermal change to represent 

areas where winter temperatures are less severe, and areas with more similar summer and winter surface 

temperatures. We hypothesized that these areas would be more suitable for P. longimembris. Additional 

vegetation covariates, vegetation greenness timing, seasonal greenness ratio, and vegetation greenness 

coefficient of variation were used to capture the yearly persistence of green vegetation, as well as areas 

with low shrub density. 

Habitat Model Results.—We selected 9 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. 

Unlike many of the other species, there were very few differences among the 18 models, and therefore 

all 18 models were included in the final averaged model, which represented 10 covariate layers. The 

covariate layers represented in the final averaged model and their relative contributions were winter 

precipitation (49 percent), topographic position (36 percent), seasonal thermal change (8 percent), 

seasonal greenness ratio (2 percent), seasonal greenness timing (2 percent), soil water stress (1 percent), 

surface texture (1 percent), mountain bases (<1 percent), and washes (<1 percent). Elevation and 

topographic position were present in all 18 models; seasonal thermal change was present in 12 models; 

seasonal greenness timing was present in 6 models; seasonal greenness ratio and soil water stress were 

present in 4 models; and surface texture, mountain bases, and washes were present in 1 model each. 

There was high agreement across all selected models (fig. 21), with very few areas showing 

standard errors greater than 0.05. The final averaged model had an AUC score of 0.78 and a BI score of 

0.75 (table 1), and was reclassified into high and low habitat suitability based on the fifth quantile of the 

observation records (break value = 0. 27), which resulted in 98,515 km2 of high habitat suitability. 

These areas were generally in flat, low elevation terrain, although as with L. trivirgata, the Baker Sink 

landform provided a notable break in habitat across the study area (fig. 22). 
57

Figure 21. Map showing standard error (σ/N) among the selected models for Perognathus longimembris (Little 

Pocket Mouse). Areas where the selected models provided similar habitat suitability scores are shown in blue, and

areas with low similarity are shown in red. Although not a measure of true model error, these agreement maps 

provide a general measure of model uncertainty across the landscape and show areas where the selected models 

agree (or disagree) with respect to habitat suitability. 
58

Figure 22. Average habitat suitability among the selected models for Perognathus longimembris (Little Pocket 

Mouse). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability 

are shown in red; areas with low suitability are shown in blue.
59

Habitat Model Discussion.—We hypothesized that P. longimembris would be associated with 

the washes and alluvial fan proximity environmental covariates, but instead found a high association 

with winter precipitation and topographic position (figs. 23A and 23B). Each of these covariates showed 

a peaked response, suggesting a preference for intermediate ranges of winter precipitation and 

topographic position, which generally corresponded to areas with alluvial fans, mildly sloped valleys 

and areas near the bases of mountains, but not valley bottoms or high drainage areas, such as dry 

lakebeds and playas. The apparent avoidance of areas with very low or very high winter precipitation 

suggests a preference for an optimal amount of winter precipitation, likely coinciding with vegetation 

communities that are tied to intermediate levels of winter precipitation (Beatley, 1976; Turner and 

Randall, 1989). Areas with moderate winter precipitation were generally in the Western Mojave Desert 

along valley washes and the bases of mountains, which receive greater winter precipitation than similar 

terrain in the eastern parts of the Mojave Desert, and less winter precipitation than in the high-elevation 

mountains. 

Alluvial fan proximity was not included in any of the selected models, even though the layer on 

which it was based (topographic position) showed the second highest relative contribution, and was 

prevalent in all selected models. The modeled response for the covariate layer seasonal thermal change

suggested decreased habitat suitability in areas with high thermal variation from summer to winter (fig.

23C), such as the valley bottoms of northern latitudes in the study area, which have very cold winters 

and warm summers. However, P. longimembris were not found on large mountain tops or ridge-like 

features at high elevations, both of which also show low thermal variation between summer and winter, 

but are instead relatively cold throughout the year. The remaining 6 covariate layers showed relative 

contributions of less than 2 percent, and each was included in fewer than 6 of the 18 selected models. 

Although these covariates provided relatively little contribution to the final averaged model, each 

showed relationships to P. longimembris habitat suitability that were not unreasonable given our 

expectations (figs. 23D–23F). 
60

Figure 23. Generalized response curves and histograms for Perognathus longimembris (Little Pocket Mouse) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for P. longimembris, where each of the covariates is held at 

their study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates 

are shown.
61

Plestidon gilberti

Conceptual Model.—P. gilberti occur in rocky boulder filled areas with decomposed granite 

(Jones and others, 1985; Stebbins, 2003). They are often found near permanent or semi-permanent water 

sources, and near the base of mountains in the Mojave Desert (Rodgers, 1944; Stebbins, 2003). To 

capture these habitat preferences, we included surface texture, rocky surfaces, cliffs, topographic 

position, and mountain bases. P. gilberti are ubiquitous across their range, occupying a wide range of 

habitats within and outside of the Mojave Desert. We included the covariate layers soil water stress, 

vegetation greenness timing, and seasonal greenness ratio to represent geographic variation in surface 

water and vegetation characteristics, and as a proxy for surface water and riparian vegetation. The 

environmental covariates washes, alluvial fan proximity, and slope were used to capture the broad range 

of desert scrub habitat, while seasonal thermal change was used to represent the thermal stability of an 

area. We hypothesized that areas along the edges of alluvial fans in low to moderate slope with more 

thermal stability would be more suitable for P. gilberti. 

Habitat Model Results.—We selected 8 models from the suite of possible GAM models to run 

as MaxEnt models, resulting in a total of 16 models that were reviewed by our collaborators. Only two 

of the MaxEnt models were included in the final averaged model, which represented six covariate 

layers—soil water stress, surface texture, mountain bases, winter precipitation, elevation, and cliffs. The 

following covariates were included in both of the selected models (relative contributions shown in 

parentheses): winter precipitation (73 percent), mountain bases (10 percent), surface texture (9 percent), 

and soil water stress (5 percent).

The final averaged model had an AUC score of 0.86 and a BI score of 0.59 (table 1). There was 

high agreement between the two selected models, but we did not calculate a standard error map due to 

the low sample size. The averaged P. gilberti habitat suitability model was reclassified into high and 

low habitat suitability, which resulted in 60,983 km2

 of high habitat suitability based on the fifth

quantile of the observation records (break value=0.14, fig. 24).
62

Figure 24. Map showing average habitat suitability among the selected models for Plestidon gilberti (Gilbert’s 

Skink). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability 

are shown in red; areas with low suitability are shown in blue.
63

Habitat Model Discussion.—We hypothesized that the covariate layers rocky surfaces,

mountain bases, and surface texture would be some of the most important predictors of P. gilberti

habitat suitability due to P. gilberti’s prevalence in areas with rocky boulders and decomposed granite. 

However, of these three layers, only mountain bases and surface texture were included in both of the 

selected models, and neither showed a high relative contribution to the final averaged model. The 

modeled response for mountain bases did, however, show a positive relationship to P. gilberti habitat 

suitability (fig. 25B), suggesting an increase in habitat suitability with increased coverage of mountain 

bases. This is concordant with current evidence suggesting that P. gilberti tend to occur in areas with 

decomposed granite near the base of mountains (Rodgers, 1944; Stebbins, 2003). Similarly, P. gilberti

habitat suitability was higher in areas where surface texture showed more rock-like substrates, 

suggesting decreased suitability in silty, sandy and gravelly areas (fig. 25C), which supported 

observations of P. gilberti in rocky habitats (Jones and others, 1985; Stebbins, 2003). 

The surface texture and rocky surfaces covariate layers were derived from MODIS Land Surface 

Temperature data, and therefore have a spatial scales of 1 km, which is likely too coarse to identify 

rocky features (that is, boulders, and decomposed granite) that P. gilberti use as habitat. Additional 

work is likely needed to develop covariate layers that accurately reflect the physical properties of the 

surface substrate at a spatial scale relevant to P. gilberti. 

Instead of describing the physical terrain of P. gilberti, the most influential covariate in the 

habitat suitability model was winter precipitation, which showed positive relationship with P. gilberti 

habitat suitability (fig. 25A). This relationship suggested an avoidance of areas with low cumulative 

winter precipitation, and a preference for areas with more than 50 mm of winter precipitation. The soil 

water stress covariate was negatively related to habitat suitability, which was likely due to the positive 

relationship between soil moisture and conifer forests and standing water (Cheng and others, 2006, 

2008), both of which are not occupied by P. gilberti.
64

Figure 25. Generalized response curves and histograms for Plestidon gilberti (Gilbert’s Skink) for environmental 

covariates. Histograms are derived from a random sample of the study area and illustrate the distribution of 

covariates across the landscape. Response curves are derived from a single model that is inclusive of all 

covariates represented in the selected models for P. gilberti , where each of the covariates is held at their study 

area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
65

Sceloporus magister

Conceptual Model.—S. magister are found in rocky terrain, often with large boulders and large 

sections of exposed rock (Turner and Wauer, 1963; Parker and Pianka, 1973; Stebbins, 2003). They 

often occur in canyons, washes, and foothills, up to areas where pinyon-juniper vegetation occurs 

(Turner and Wauer, 1963; Parker and Pianka, 1973). We included the environmental covariates surface 

texture, rocky surfaces, slope, alluvial fan proximity, washes, mountain bases, cliffs, and topographic 

position to capture desert spiny lizards’ preference for rocky habitats, canyons, and foothills. 

S. magister also can be found in upland habitats as well, usually near trees or other large 

vegetation (Parker and Pianka, 1973) in areas characterized by creosote bush, mesquite, burro bush, and 

salt bush (among others). Although we did not have covariate layers describing vegetation communities 

directly, we used vegetation greenness, vegetation greenness timing, seasonal greenness ratio, and soil 

water stress as proxies for distinguishing vegetation types among Mojave Desert habitats. We also 

included winter precipitation because of the positive relationship between winter precipitation and 

spring plant production (Beatley, 1976; Turner and Randall, 1989). 

Habitat Model Results.—We selected 9 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. 

Seven MaxEnt models were included in the final averaged model, which represents eight covariate 

layers—topographic position, seasonal thermal change, elevation, vegetation greenness timing, 

mountain bases, washes, rocky surfaces, and soil water stress. Relative contributions for the eight 

covariate layers were 51 percent (elevation), 15 percent (topographic position), 11 percent (seasonal 

thermal change), 6 percent (mountain bases), 6 percent (rocky surfaces), 4 percent (vegetation greenness 

timing), 4 percent (washes), and 2 percent (soil water stress), and were included in 7, 4, 5, 2, 3, 4, 3, and 

1 of the 18 possible models, respectively. 

Generally, there was high agreement among the selected models (fig. 26), with few areas 

showing standard errors greater than 0.05. The final averaged model had an AUC score of 0.72 and a BI 

score of 0.33 (table 1), and was reclassified into high and low habitat suitability, and resulted in 111,731 

km2

 of high habitat suitability based on the fifth quantile of the observation records (break value=0. 31, 

fig. 27).
66

Figure 26. Map showing standard error (σ/N) among the selected models for Sceloporus magister (Desert Spiny 

Lizard). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability. 
67

Figure 27. Average habitat suitability among the selected models for Sceloporus magister (Desert Spiny Lizard). 

Values range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability are 

shown in red; areas with low suitability are shown in blue.
68

Habitat Model Discussion.—We hypothesized that S. magister habitat suitability would be 

related to washes, rocky terrain, mountain bases, or other areas with large boulders. Although each of 

these covariate layers were present in up to three of our selected models, none were most prevalent 

among the covariates represented by the final averaged model. This suggests that while washes, rocky 

terrain, and mountain bases impart some explanatory power on S. magister habitat suitability, these 

three layers do not fully explain the geographic variation in habitat suitability. Instead, elevation was the 

single most prevalent covariate, and also had the highest relative contribution of all environmental 

covariates, suggesting that elevation alone accounted for more of the variation in geographic distribution 

of S. magister habitat suitability than any of the other environmental covariates. Habitat suitability 

decreased sharply at elevations greater than 1,500 m, with a peak elevation between 700 and 1,400 m 

(fig. 28A), which was unlike many descriptions of S. magister habitat that suggest high prevalence 

down to sea level (Stebbins, 2003; Jones and Lovich, 2009). The model response for the topographic 

position layer confirmed current hypotheses of S. magister habitat by suggesting and avoidance of areas, 

such as valley bottoms or alluvial fan apexes, as well as an avoidance of high ridge tops and high 

elevation areas (fig. 28B).

Seasonal thermal change also was a strong predictor of S. magister habitat suitability, as 

evidenced by its inclusion in all but one of the selected models and relative contribution of 11 percent. 

As with elevation and topographic position, seasonal thermal change showed a peak response in the 

middle of the range found in the Mojave Desert, suggesting decreased suitability of areas with 

extremely low or extremely high seasonal thermal differences (fig. 28C). These areas were correlated 

with areas on mountain tops in the southern latitudes and valley bottoms in the more northern latitudes. 

Although derived as a measure of the seasonal change in surface radiant temperature, seasonal thermal 

change is related to surface material type, slope, aspect and latitude, and reflects seasonal change in the 

thermal environment that an organism can experience. As with many of the environmental covariates 

used in the habitat models presented here, the spatial scale at which seasonal thermal change was 

derived precludes analyses of microhabitat or fine scale resource selection, such as the preference for a 

particular type of boulder or specific thermal operative temperature. 

Although the mountain bases covariate was present in only two of the selected models, it 

showed a relative contribution of 6 percent, and suggested an avoidance of areas classified as mountain 

bases (fig. 28D). This was unexpected, considering the prevalence of S. magister in areas with boulders 

and other large rocks, which are often located at the bases of mountains. However, many of the areas we 

identified as mountain bases were at elevations greater than 2,000 m, well above the range that S. 

magister are typically found. Similarly, the rocky surfaces layer showed an avoidance of areas with of 

rocky surfaces (fig. 28E), which was contrary to observations of S. magister’s preference for boulders, 

and other large structures, such as trees, shrubs, and debris piles (Jones and Lovich, 2009). The spatial 

scale at which the mountain base and rocky surface layers likely were derived is too coarse to represent 

the habitat features (boulders and other rocky surfaces) that S. magister use for habitat, and as such, did 

not adequately capture S. magister habitat suitability. Habitat models derived at finer spatial scales may 

see covariates, such as rocky surfaces, mountain bases, and surface texture, have greater importance 

than general indicators, such as elevation.
69

Figure 28. Generalized response curves and histograms for Sceloporus magister (Desert Spiny Lizard) for 

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for S. magister, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
70

Thomomys bottae

Conceptual Model.—Thomomys bottae occur in a wide range of habitats, but are primarily 

found in areas that support burrowing, such as sandy or gravelly soils (Zeveloff, 1988). We included 

surface texture and sandy surfaces to capture these habitat preferences. T. bottae are also often found in 

areas with alluvial soils that can support grasses and forbs for forage (Linsdale, 1938; Bond, 1946; Fitch 

and Bentley, 1949; Smallwood and Morrison, 1999). We included alluvial fan proximity, vegetation 

greenness timing, seasonal greenness ratio, soil water stress, and winter precipitation to represent areas 

with high annual plant production. T. bottae are not often found in extremely rocky terrain, but can 

occur in meadows at high elevations (Zeveloff, 1988). We were unable to derive a layer describing high 

elevation meadows, but instead included elevation and slope to capture high elevation areas with low 

slope. We also included washes, topographic position, mountain bases, and seasonal thermal change

capture the broad characteristics of T. bottae habitat. 

Habitat Model Results.—We selected 9 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. 

Sixteen of them were included in the final averaged model, which represented eight covariate layers—

washes, slope, seasonal thermal change, winter precipitation, vegetation greenness timing, seasonal 

greenness ratio, topographic position, and elevation. Relative contributions for the covariate layers were 

61 percent (winter precipitation), 14 percent (topographic position), 8 percent (slope), 7 percent

(seasonal greenness timing), 5 percent (elevation), 4 percent (seasonal thermal change), 2 percent

(seasonal greenness ratio), and <1 percent (washes), with inclusion in 8, 5, 3, 7, 4, 4, 3, and 1 of the 16

models, respectively.

There was moderate agreement across all selected models (fig. 29), although many areas at high 

elevation had standard errors of greater than 0.1. This was caused primarily by the GAM models, which 

tended to predict more suitable habitat at higher elevations than the MaxEnt models. The final averaged 

model had an AUC score of 0.83 and a BI score of 0.51 (table 1), and was reclassified into high and low 

habitat suitability, which resulted in 104,675 km2 of high habitat suitability based on the fifth quantile of 

the observation records (break value=0.14, fig. 30).
71

Figure 29. Map showing standard error (σ/N) among the selected models for Thomomys bottae (Botta’s Pocket 

Gopher). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability.
72

Figure 30. Map showing average habitat suitability among the selected models for Thomomys bottae (Botta’s 

Pocket Gopher). Vaules range from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
73

Habitat Model Discussion.—Contrary to our hypotheses of T. bottae habitat, surface texture 

and sandy surfaces were not among the selected models, and instead, winter precipitation was the most 

prevalent and had the highest relative contribution. T. bottae use spring grasses and forbs for forage, 

which in the Mojave Desert, are influenced heavily by winter precipitation (Beatley, 1976; Turner and 

Randall, 1989). The response curve for winter precipitation suggested a preference for areas with greater 

than 50 mm of winter precipitation (fig. 31A), which were often areas at high elevations within the 

Mojave Desert. The response curves for the topographic position layer suggested an avoidance of valley 

bottoms and playas, as well as an avoidance of mountain tops and ridges (fig. 31B). This supports 

current hypotheses of T. bottae habitat preference for areas with sandy or gravelly soils, as well as 

meadows at high elevations, and decreased suitability of extensively rocky terrain or fine silty soils that 

are found in playas. The modeled response to the slope layer suggested a similar decrease of habitat 

suitability in steep mountainous terrain (fig. 31C), and a preference for flatter terrain, although the latter 

was slightly contradictory to the modeled response of the topographic position layer, which suggested 

an avoidance of certain flat areas (that is, dry lakebeds, and valley bottoms). This apparent difference in 

habitat preference between topographic position and slope is likely due to the fact that not all areas with 

low slope are valley bottoms or dry lakebeds, and that areas adjacent to dry lakebeds often have low 

slope and alluvial soils suitable for fossorial activities. Habitat suitability for T. bottae showed a positive 

trend with seasonal greenness timing (fig. 31D), suggesting that T. bottae occur in areas where peak 

vegetation density occurs in mid-summer rather than in early spring. These areas tended to be high 

elevation vegetation communities as well as riparian communities, which often maintained greenness 

later in the season than surrounding vegetation communities. 
74

Figure 31. Generalized response curves for Thomomys bottae (Botta’s Pocket Gopher) for environmental 

covariates. Histograms are derived from a random sample of the study area and illustrate the distribution of 

covariates across the landscape. Response curves are derived from a single model that is inclusive of all 

covariates represented in the selected models for T. bottae, where each of the covariates is held at their study area 

mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are shown.
75

Uma scoparia

Conceptual Model.—U. scoparia are arenicolous lizards that occur in areas with very sandy 

soils, such as windblown sand dunes (Stebbins, 1944; Pianka, 1967; Fromer and others, 1983) and areas 

with high day-time temperatures, but milder night-time temperatures (Pianka, 1967; Phil Medica, U.S. 

Geological Survey, oral commun., 2011). We included surface texture, sandy surfaces, and seasonal 

thermal change to capture the geographic variability in these habitat preferences. We had hypothesized 

that areas with more thermal stability would be more suitable for U. scoparia due to the additional 

energetic costs (and opportunistic costs) of thermoregulation when temperatures are above or below 

optimal. Across their range, U. scoparia tend to occur in areas with elevations less than 1,000 m, and 

are not found in areas with steep slopes (Stebbins, 1944; Mayhew, 1964; Pianka, 1967; Fromer and 

others, 1983; Jones and Lovich, 2009). We therefore included elevation and slope to represent these 

habitat preferences. We also included washes, mountain bases, topographic position, winter 

precipitation, vegetation greenness timing, seasonal greenness ratio, soil water stress, and alluvial fan 

proximity to capture other potential variation in habitat for U. scoparia.

Habitat Model Results.—We selected 9 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. All 

9 of the MaxEnt models but none of the GAM models were included in the final averaged model, which 

represented 10 covariate layers. The covariate layers represented in the final averaged model and their 

relative contributions were elevation (34 percent), sandy surfaces (27 percent), surface texture (18 

percent), seasonal thermal change (7 percent), soil water stress (6 percent), slope (3 percent), 

topographic position (2 percent), winter precipitation (1 percent), (seasonal greenness ratio (1 percent), 

and washes (<1 percent). Only 1 of the environmental covariates was present in all 9 models (elevation);

surface texture was included in 7 models; sandy surfaces, seasonal thermal change, and soil water stress

were each included in 6 models; slope was included in 4 models; topographic position was included in 3

models; seasonal greenness ratio was included in 2 models; and winter precipitation and washes were 

included in 1 model each. 

The final averaged model had an AUC score of 0.90 and a BI score of 0.54 (table 1). There was 

a fair amount of disagreement among the selected models, especially in areas of low habitat suitability,

such as high elevation or near the base of mountains (fig. 32). The averaged U. scoparia habitat 

suitability model was reclassified into high and low habitat suitability, and resulted in 30,127 km2 of 

high habitat suitability based on the fifth quantile of the observation records (break value=0.17; fig. 33).
76

Figure 32. Map showing standard error (σ/N) among the selected models for Uma scoparia (Mojave Fringe-toed 

Lizard). Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas 

with low similarity are shown in red. Although not a measure of true model error, these agreement maps provide a 

general measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability.
77

Figure 33. Map showing average habitat suitability among the selected models for Uma scoparia (Mojave Fringetoed Lizard). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high 

suitability are shown in red; areas with low suitability are shown in blue.
78

Habitat Model Discussion.—U. scoparia have a fringe of small scales on their third and fourth

hind foot toes, enabling them to run quickly on fine, windblown sand (Stebbins, 1944). They also have 

extensive adaptations enabling them to burrow down into sand for thermal refugia, such as a counter 

sunk lower jaw and valved nostrils (Stebbins, 1944). Given these highly specialized traits for living in 

sandy habitats, the high prevalence and relative contributions of the sandy surfaces and surface texture

layers was expected, and each suggested a preference for areas with sandy and finer grained surfaces 

(figs. 34B and 34C). In terms of elevation, which provided the single largest relative contribution and 

was present in all of the nine selected models, U. scoparia habitat suitability showed a sharp decrease

with increasing elevation (fig. 34A). Specifically, habitat suitability was depressed at elevations greater 

than 1,000 m, which was consistent with the previously suggested cutoff for habitat (Jones and Lovich,

2009). Unlike the other covariates, seasonal thermal change showed a narrow range for suitable habitat 

relative to the range of values available in the study area (fig. 34D), although with a small relative 

contribution (7 percent), seasonal thermal change did not provide much influence on U. scoparia habitat 

suitability. These values corresponded to areas in the southwestern parts of the Mojave Desert, generally 

in valleys distinct from mountainous terrain. The modeled response to soil water stress suggested that U. 

scoparia habitat suitability increases with increased soil water stress (fig. 34E), but as with seasonal 

thermal change, provided little influence on habitat suitability throughout the Mojave Desert. The 

remaining five environmental covariates provided less than 5 percent relative contribution each, and are 

not discussed.
79

Figure 34. Generalized response curves and histograms for Uma scoparia (Mojave Fringe-toed Lizard) for

environmental covariates. Histograms are derived from a random sample of the study area and illustrate the 

distribution of covariates across the landscape. Response curves are derived from a single model that is inclusive 

of all covariates represented in the selected models for U. scoparia, where each of the covariates is held at their 

study area mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are 

shown.
80

Xantusia vigilis

Conceptual Model.—X. vigilis occur in flat areas with an abundance of surface litter, often in 

stands of Joshua trees and in areas where substantial surface litter is present, often comprised of dead 

cactus or other woody debris (Stebbins, 1948; Zweifel and Lowe, 1966), especially in western portions 

of the Mojave Desert. The presence of shade or vertical structures capable of providing shade during the 

hotter summer seasons are usually present in areas occupied by X. vigilis (Zweifel and Lowe, 1966). 

Although we were unable to create a layer characterizing the abundance of surface litter, we included 

several layers representing vegetation characteristics (seasonal greenness ratio, vegetation greenness 

timing, and soil water stress) to approximate the presence of shade vegetation. Generally, X. vigilis tend 

to occur in areas with coarse sand, gravelly soils, or silt-like substrates with organic material where 

topographic slope is minimal (Stebbins, 1948; Miller, 1951; Zweifel and Lowe, 1966) and is often 

associated with valleys, low hills, and in areas where annual precipitation is relatively high for desert 

environments (Miller, 1951; Stebbins, 2003). We included surface texture, sandy surfaces, slope, 

mountain bases, washes, topographic position, alluvial fan proximity, and winter precipitation to capture

these habitat preferences. 

Habitat Model Results.—We selected 9 models from among the suite of possible GAM models 

to run as MaxEnt models, resulting in a total of 18 models that were reviewed by our collaborators. All 

nine of the MaxEnt models were included in the final averaged model, which represented 10 covariate 

layers. The covariate layers represented in the final averaged model and their relative contributions were 

elevation (55 percent), topographic position (14 percent), slope (11 percent), soil water stress (7

percent), seasonal thermal change (5 percent), mountain bases (3 percent), surface texture (2 percent), 

washes (1 percent), seasonal greenness timing (1 percent), and seasonal greenness ratio (<1 percent),

and each covariate layer was included in 9, 4, 4, 4, 2, 3, 1, 2, 2, and 1 model(s), respectively. The large 

number of environmental covariates present in three or fewer of the selected models is indicative of 

poor performing models.

However, even with low performing models, we found high agreement across all selected 

models (fig. 35), and in particular, areas with high habitat suitability showed the most amount of 

agreement. Standard errors ranged from 0.001 to 0.11. Very few areas were present with standard errors 

greater than 0.1. The final averaged model had an AUC score of 0.78 and a BI score of 0.05, which was 

the lowest among our study species. The averaged model was reclassified into high and low habitat

suitability, resulting in 106,994 km2 of high habitat suitability based on the fifth quantile of the 

observation records (break value=0.29, fig. 36). As with many of the other species included in this 

report, a notable break in habitat suitability occurs near the Baker Sink landform. 
81

Figure 35. Map showing standard error (σ/N) among the selected models for Xantusia vigilis (Desert Night Lizard). 

Areas where the selected models provided similar habitat suitability scores are shown in blue, and areas with low 

similarity are shown in red. Although not a measure of true model error, these agreement maps provide a general 

measure of model uncertainty across the landscape and show areas where the selected models agree (or 

disagree) with respect to habitat suitability.
82

Figure 36. Map showing average habitat suitability among the selected models for Xantusia vigilis (Desert Night 

Lizard). Values range from 0~1, and correspond to an index of habitat suitability, where areas with high suitability 

are shown in red; areas with low suitability are shown in blue.
83

Habitat Model Discussion.—We were unable to derive an environmental covariate representing 

woody surface debris or Joshua tree communities, both of which are thought to be important predictors 

of X. vigilis habitat. Instead, we used measures of terrain topography (primarily slope and topographic 

position) and vegetation characteristics (seasonal greenness ratio, vegetation greenness timing, 

vegetation greenness, and soil water stress) to represent these habitats types. Although none of these 

covariates provided the greatest relative contribution or were the most prevalent among our selected 

models, topographic position was present in slightly fewer than one-half of the selected models, and 

provided 14 percent of the relative contribution in the averaged model. The model response to 

topographic position suggested decreased suitability on mountain tops and ridges, and a preference for 

valleys and alluvial fans, with a slight decrease in habitat suitability at high values, representing playas 

(fig. 37B). This decrease in playas is congruent with observations of X. vigilis habitat suggesting a 

requirement of surface litter, which is often not present in playas. The slope layer was used in nearly 

one-half of the selected models, and showed a relative contribution of only 11 percent. Counter to 

hypotheses of X. vigilis habitat, and the response curve suggested a preference for low slope areas, and a 

preference for high slopes instead, although the increase in habitat suitability was less pronounced at 

higher slopes than lower slopes (fig. 37C). 

Instead of the hypothesized covariates describing terrain or vegetation, elevation was the most 

prevalent covariate and provided the highest relative contribution to the averaged model. Additionally, 

one of the selected models contained only elevation as a single covariate, suggesting that elevation was 

by far the single best predictor of X. vigilis habitat in the Mojave Desert. The modeled response to 

elevation suggested a preference for areas at high elevations, and an avoidance of low elevations,

although areas at elevations greater than 3,000 m also showed a decrease in habitat suitability (fig.

37A). The layer soil water stress, although providing only 7 percent relative contribution to the averaged 

model, was used in 4 of the 9 selected models, and the response curve suggested a slight preference for 

areas with less surface water (fig. 37D). This was counter to previous observations of X. vigilis, which 

suggest that habitat suitability is likely higher in areas with greater surface water, although a general 

avoidance of dense riparian vegetation is also known. Winter precipitation was not among the 

environmental covariates included in the averaged model, which also contradicted many previous 

hypotheses of X. vigilis habitat stating that areas with high precipitation will provide better habitat than 

areas with low precipitation (Miller, 1951; Stebbins, 2003). More work is needed to identify how X. 

vigilis respond to winter precipitation and, particularly, why winter precipitation was omitted from the 

suite of selected models.

Our model of X. vigilis habitat had a test AUC of 0.78, which was one of the lowest of the 

species modeled here. In conjunction with elevation as the most prevalent and important predictor of X. 

vigilis habitat, this low AUC is indicative of the complex habitat requirements of X. vigilis, and general 

difficulty in modeling them. Additional covariate layers representing Joshua tree communities and 

surface litter are likely needed to improve models of X. vigilis habitat. However, as is often the case for 

layers representing habitat preferences, spatial scale is an important issue due to the discrepancy 

between microhabitat features (for example, organic top soil and surface litter) and landscape features 

(for example, slope and topographic features) used in these models. 
84

Figure 37. Generalized response curves and histograms for Xantusia vigilis (Desert Night Lizard) for environmental 

covariates. Histograms are derived from a random sample of the study area and illustrate the distribution of 

covariates across the landscape. Response curves are derived from a single model that is inclusive of all 

covariates represented in the selected models for X. vigilis, where each of the covariates is held at their study area 

mean. The dotted line shows linear predictor value of 0. Only the six highest contributing covariates are shown.
85

Discussion

We modeled habitat suitability in the Mojave Desert Basin for 13 species using Generalized 

Additive Models and Maximum Entropy methods to define their potential distributions within this 

region. The primary sources of observation records for this modeling effort were publically available 

museum records or opportunistic collection records and independent databases provided by our 

collaborators. Due to the lack of known absence localities in these sources, we chose to use the two 

afore-mentioned methods to model habitat suitability. However, although these two approaches are 

prevalent in the modeling literature, they are not without error, and have limitations when heterogeneity 

in the source of observation records is large (Phillips and Dudik, 2008; Fanklin, 2010; Newbold, 2010). 

They provide the best quantification of habitat when a single sampling strategy is used, and when 

adequate samples of the geographic and environmental space of interest are available (Elith and 

Leathwick, 2009). We were unable to use a single sampling strategy, and instead used many sources of 

occurrence data, some of which were opportunistic collection records and others of which were 

observations from specific research efforts with independent sampling strategies. This resulted in 

considerable heterogeneity in the sampling design and methodology in the assimilation of our input 

datasets. There also were a wide range of temporal and spatial sampling regimes across the data sources, 

all of which added variability to the species-environment relationships we were modeling. The 

observation records we obtained were often geographically aggregated, resulting in large areas where 

we had no evidence of presence (or absence). In many of these areas, our collaborators had personal 

knowledge of occurrence generally, but no formal observation records were available. For these species, 

we often selected models that predicted habitat broadly in these areas, thus ensuring that our models 

were not over-fit to our limited and geographically biased observation records.

Another drawback to the heterogeneity in our observation records is the time span over which 

data were collected. When this time span is large, a species may have the opportunity to shift its 

distribution (Elith and others, 2010). Similarly, the landscape may change through time (for example, 

disturbance, fire, climate change), such that a species may no longer be present in the location where it 

once was. All these factors contribute noise to the species-environment signal, and make quantifying 

habitat difficult. Where sufficient locality records were available, we limited the temporal period of 

occurrence data to 1950–2011. However, for several species, such as L. trivirgata and P. gilberti, which 

lacked sufficient observations, we chose to include records dating back to 1900 in order to increase our 

sample size. Similarly, for many of the occurrence records we did not know the collection date, but 

chose to include those records to increase sample sizes.
86

In general, most habitat models developed here had moderate to good discrimination ability with 

AUC scores ranging from 0.72 to 0.90. Although some of these values are not as high as models for

other species (see Elith and Graham, 2009, and Jarnevich and Laubhan, 2011, but see Elith and 

Leathwick, 2009, for similar AUC scores), we think that these values are appropriate given the spatial 

scale (1 km) and the inherent sampling bias from the large assortment of museum locality data used 

(Chefaoui and Lobo, 2008; Phillips and others, 2009; Newbold, 2010). Even though we took measures 

to avoid over-fitted models, such as replicating the model training (bootstrapping), limiting the number 

of covariate layers to five per model, and not including interactions among input covariates (although 

MaxEnt inherently includes interactions of all covariates; Phillips and Dudik, 2008), we were cautious 

to avoid inferring too much ecological significance from any single model. We also looked at the 

standard error in habitat suitability scores across the selected models for each of the focal species to 

identify geographic areas where the suite of selected models did not agree. Although our selected habitat 

models tended to disagree more in playas and areas with high elevation, they were, in general, very 

consistent. Other methods have been proposed to avoid over-fitted MaxEnt models (see Anderson and 

Gonzalez, Jr., 2011); however, because we used model averaging methods to construct consensus 

models for our species, we were unable to implement these methods and we believe the methods used 

here are adequately robust and generalized to the Mojave Basin and Range ecoregion in which they 

were developed.

References Cited

Anderson, D., Burnham, K., and Thompson, W., 2000, Null hypothesis testing—Problems, prevalence, 

and an alternative: Journal of Wildlife Management, v. 64, p. 912–923.

Anderson, R.P., and Gonzalez, Jr, I., 2011, Species-specific tuning increases robustness to sampling 

bias in models of species distributions—An implementation with Maxent: Ecological Modelling, v. 

222, p. 2,796–2,811.

Araújo, M.B., and Williams, P.H., 2000, Selecting areas for species persistence using occurrence data: 

Biological Conservation, v. 96, p. 331-45.

Arizona Game and Fish Department, 2003, Lichanura trivirgata trivirgata in Abstract compiled and 

edited by the Heritage Data Management System: Phoenix, Arizona Game and Fish Department, p. 5.

Austin, M., 2002, Spatial prediction of species distribution—An interface between ecological theory and 

statistical modeling: Ecological Modelling, v. 157, p. 101–118.

Bailey, V., 1931, Mammals of New Mexico: North American Fauna, v. 53, p. 1–412.

Bartholomew, G.A., and Cade, T.J., 1957, Temperature regulation, hibernation, and aestivation in the 

little pocket mouse, Perognathus longimembris: Journal of Mammalogy, v. 38, p. 60–72.

Beatley, J.C., 1976, Rainfall and fluctuating plant populations in relation to distributions and numbers of 

desert rodents in southern Nevada: Oecologia, v. 24, p. 21–42.

Bond, R.M., 1946, The breeding habits of Thomomys bottae in Orange County, California: Journal of 

Mammalogy, v. 27, p. 172–174.

Boyce, M., Vernier, P., Nielsen, S., and Schmiegelow, F., 2002, Evaluating resource selection 

functions: Ecological Modelling, v. 157, p. 281–300.
87

Bradford, D.F., Neale, A.C., Nash, M.S., Sada, D.W., and Jaeger, J.R, 2003, Habitat patch occupancy 

by toads (Bufo punctatus) in a naturally fragmented desert landscape: Ecology, v. 84, p. 1,012–1,023.

Brown, D., 1982, Desert Plants; Biotic communities of the American Southwest—United States and 

Mexico: Superior, Arizona, Boyce Thompson Southwestern Arboretum, p. 342.

Burnham, K.P., and Anderson, D.R, 2002, Model selection and multimodel inference—A Practical 

Information-Theoretic Approach (2nd ed.): Springer-Verlag, p. 488, ISBN 0-387-95364-7.

Carlson, T.N., Dodd, J.K., Benjamin, S.G., and Cooper, J.N., 1981, Satellite estimation of the surfaceenergy balance, moisture availability and thermal inertia: Journal of Applied Meteorology, v. 20, p. 

67–87.

Chefaoui, R.M., and Lobo, J.M., 2008, Assessing the effects of pseudo-absences on predictive 

distribution model performance: Ecological Modelling, v. 210, p. 478–486.

Cheng, Y.B., Ustin, S.L., Riano, D., and Vanderbilt, V.C., 2008, Water content estimation from 

hyperspectral images and MODIS indexes in Southeastern Arizona: Remote Sensing of Environment, 

v. 112, p. 363–374.

Cheng, Y.B., Zarco-Tejada, P.J., Riaño, D., Rueda, C.A., and Ustin, S.L., 2006, Estimating vegetation 

water content with hyperspectral data for different canopy scenarios—Relationships between AVIRIS 

and MODIS indexes: Remote Sensing of Environment, v. 105, p. 354–366.

Chew, R.M., and Butterworth, B.B, 1964, Ecology of rodents in Indian Cove (Mojave Desert), Joshua 

Tree National Monument California: Journal of Mammalogy, v. 45, p. 203–225.

Collingham, Y., Wadsworth, R., Huntley, B., and Hulme, P., 2000, Predicting the spatial distribution of 

non-indigenous riparian weeds—Issues of spatial scale and extent: Journal of Applied Ecology, v. 37, 

p. 13–27.

Commission for Environmental Cooperation Working Group, 1997, Ecological regions of North 

America—Toward a common perspective: Montreal, Quebec, Commission for Environmental 

Cooperation, p. 71.

Creusere, F.M., and Whitford, W.G., 1976, Ecological relationships in a desert anuran community: 

Herpetologica, v. 32, p. 7–18.

Crews, S.C., 2009, Assessment of rampant genitalic variation in the spider genus Homalonychus

(Araneae, Homalonychidae): Invertebrate Biology, v. 128, p.107–125.

Crews, S., and Hedin, M., 2006, Studies of morphological and molecular phylogenetic divergence in 

spiders (Araneae: Homalonychus) from the American southwest, including divergence along the Baja 

California Peninsula: Molecular Phylogenetics and Evolution, v. 38, p. 470-487.

Crimmins, T.M., Crimmins, M.A., and Bertelsen, C.D., 2009, Flowering range changes across an 

elevation gradient in response to warming summer temperatures: Global Change Biology, v. 15, p.

1,141–1,152.

Cumming, G.S., 2000, Using between-model comparisons to fine-tune linear models of species ranges:

Journal of Biogeography, v. 27, p. 441–455.

Dayton, G.H., and Fitzgerald, L.A., 2006, Habitat suitability models for desert amphibians: Biological 

Conservation, v. 132, p. 40–49.

Domínguez, K., and Jiménez, M.L., 2005, Mating and self-burying behavior of Homalonychus 

theologus Chamberlin (Araneae, Homalonychidae) in Baja California Sur: Journal of Arachnology, v. 

33, p. 167–174.

Elith, J., and Graham, C.H., 2009, Do they? How do they? WHY do they differ? On finding reasons for 

differing performances of species distribution models: Ecography, v. 32, p. 66–77.
88

Elith, J., Graham, C.H., Anderson, R.P., Dudík, M., Ferrier, S., Guisan, A., Hijmans, R.J., Huettmann,

F., Leathwick, J.R., Lehmann, A., Li, J., Lohmann, L.G., Loiselle, B.A., Manion, G., Moritz, C., 

Nakamura, M., Nakazawa, Y., Overton, J.McC.M., Townsend Peterson, A., Phillips, S.J., Richardson,

K., Scachetti-Pereira, R., Schapire, R.E., Soberón, J., Williams, S.M., Wisz, S., and Zimmermann, 

N.E., 2006, Novel methods improve prediction of species’ distributions from occurrence data:

Ecography, v. 29, p. 129–151, doi: 10.1111/j.2006.0906-7590.04596.x.

Elith, J., Kearney, M., and Phillips, S., 2010, The art of modelling range-shifting species: Methods in 

Ecology and Evolution, v. 1, p. 330–342.

Elith, J., and Leathwick, J.R., 2007, Predicting species distributions from museum and herbarium 

records using multiresponse models fitted with multivariate adaptive regression splines: Diversity and 

Distributions, v. 13, p. 265–275. 

Elith, J., and Leathwick, J.R., 2009, Species Distribution Models—Ecological Explanation and 

Prediction Across Space and Time: Annual Review of Ecology Evolution and Systematics, p. 677-

697.

Elith, J., Phillips, S.J., Hastie, T., Dudik, M., Chee, Y.E., and Yates, C.J., 2011, A statistical explanation 

of MaxEnt for ecologists: Diversity Distributions, v. 17, p. 43–57.

Engler, R., Guisan, A., and Rechsteiner, L., 2004, An improved approach for predicting the distribution 

of rare and endangered species from occurrence and pseudo-absence data: Journal of Applied 

Ecology, v. 41, p. 263–274, doi: 10.1111/j.0021-8901.2004.00881.x.

Esler, K., and Rundel, P., 1999, Comparative patterns of phenology and growth form diversity in two 

winter rainfall deserts—The Succulent Karoo and Mojave Desert ecosystems: Plant Ecology, v. 142, 

p. 97–104.

Fensholt, R., and Sandholt, I., 2003, Derivation of a shortwave infrared water stress index from MODIS 

near- and shortwave infrared data in a semiarid environment: Remote Sensing of Environment, v. 87, 

p. 111–121.

Ferguson, J.H., and Lowe, C. H., 1969, Evolutionary relationships in Bufo punctatus group: American 

Midland Naturalist, v. 81, p. 435–466.

Fielding, A., and Bell, J., 1997, A review of methods for the assessment of prediction errors in 

conservation presence/absence models: Environmental Conservation, v. 24, p. 38–49.

Fitch, H.S., and Bentley, J.J., 1949, Use of California annual plant forage by rodents: Ecology, v. 30, p. 

306–321.

Franklin, J., 1995, Predictive vegetation mapping—Geographic modelling of biospatial patterns in 

relation to environmental gradients: Progress in Physical Geography, v. 19, p. 474–499.

Franklin, J., 2010, Mapping species distributions: spatial inference and prediction. Cambridge 

University Press, Cambridge, UK.

Frey, J.K., 2009, Distinguishing range expansions from previously undocumented populations using 

background data from museum records: Diversity and Distributions, v. 15, p. 183–187.

Findley, J.S., Harris, A.H., Wilson, D.E., and Jones, C., 1975, Mammals of New Mexico, University of 

New Mexico Press, Albuquerque.

Fromer, P.S., Dodero, M., and Patterson, C., 1983, A Population Study of the Mojave fringe-toed lizard 

(Uma scoparia) on the Twentynine Palms MCAGCC: Prepared for the Natural Resources Office 

Marine Corps Air Ground Combat Center Twentynine Palms, CA 92278.

Fry J., Xian, G., Jin, S., Dewitz, J., Homer, C., Yang, L., Barnes, C., Herold, N., and Wickham, J., 2011,

Completion of the 2006 National Land Cover Database for the conterminous United States:

Photogrammetric Engineering and Remote Sensing, v. 77, p. 858-864.
89

Funk, R.S., 1967, A new colubrid snake of the genus Chionactis from Arizona: The Southwestern 

Naturalist, v. 12, p. 180–188.

Graham, C.H., Ferruer, S., Huettman, F., Moritz, C., and Peterson, A.T., 2004, New developments in 

museum-based informatics and applications in biodiversity analysis: Trends in ecology and evolution, 

v. 19, p. 497–503.

Guisan, A., and Thuiller, W., 2005, Predicting species distribution—Offering more than simple habitat 

models: Ecology Letters, v. 8, p. 993–1,009.

Guisan, A., and Zimmermann, N., 2000, Predictive habitat distribution models in ecology: Ecological 

Modelling, v. 135, p. 147–186.

Hall, E.R., 1981, The mammals of North America—Volume 1 (2nd ed.): New York, John Wiley & 

Sons, p. 1,300.

Hirzel, A.H., Hausser, J., Chessel, D., and Perrin, N., 2002, Ecological-niche factor analysis—How to 

compute habitat-suitability maps without absence data?: Ecology, v. 83, p. 2,027–2,036.

Hoffmeister, D.F., 1986, Mammals of Arizona: Tucson, Ariz., University of Arizona Press, 602 p.

Jameson, Jr., E.W., and Peeters, H.J., 1988, California Mammals: Barkeley, Calif., University of 

California Press, 403 p.

Jarnevich, C.S., and Laubhan, M.K., 2011, Balancing energy development and conservation—A 

Method Utilizing Species Distribution Models: Environmental Management, v. 47, p. 926–936.

Jones, K.B., Kepner, L.P., and Martin, T.E., 1985, Species of reptiles occupying habitat islands in 

western Arizona—A deterministic assemblage: Oecologia, v. 66, p. 595–601.

Jones, L.C., and Lovich, R.E., 2009, Lizards of the American Southwest—A Photographic Field Guide: 

Tucson, Ariz., Rio Nuevo Publishers, p. 567.

Kahle, A.B., 1987, Surface emittance, temperature and thermal inertia derived from thermal infrared 

multispectral scanner (TIMS) data for Death Valley, California: Geophysics, v. 52, p. 858–874.

Klauber, L.M., 1931, A new subspecies of the California boa, with notes on the genus Lichanura: 

Transactions of the San Diego Society of Natural History, v. 6, p. 305–318.

Klauber, L.M., 1951, The shovel-nosed snake, Chionactis, with descriptions of two new subspecies: 

Transactions of the San Diego Society of Natural History, v. 11, p. 141–204.

Krekorian, C.O., 1977, Homing in the desert iguana, Dipsosaurus dorsalis: Herpetologica, v. 33, p. 

123–127.

Linsdale, J.M., 1938, Environmental responses of the vertebrates in the Great Basin: American Midland 

Naturalist, v. 19, p. 1–206.

Liu, C., Berry, P., Dawson, T., and Pearson, R., 2005, Selecting thresholds of occurrence in the 

prediction of species distributions: Ecography, v. 28, p. 385–393.

Lobo, J.M., and Tognelli, M.F., 2011, Exploring the effects of quantity and location of pseudo-absences 

and sampling biases on the performance of distribution models with limited point occurrence data:

Journal for Nature Conservation, v. 19, p. 1–7.

Mateo, R.G., Croat, T.B., Felicísimo, Á.M., and Muñoz, J., 2010, Profile or group discriminative 

techniques?—Generating reliable species distribution models using pseudo-absences and target-group 

absences from natural history collections: Diversity and Distributions, v. 16, p. 84–94.

Mayhew, W.W., 1964, Taxonomic status of California populations of the lizard genus Uma: 

Herpetologica, v. 20, p. 170–183.

McClanahan, L.L., Ruibal, R., and Shoemaker, V.H., 1994, Frogs and toads in deserts: Scientific 

American, v. 270, p. 82–88.

McGuire, J.A., 1996, Phylogenetic systematics of crotaphytid lizards (Reptilia: Iguania: Crotaphytidae):

Bulletin of Carnegie Museum of Natural History, v. 32, p. 1–120.
90

Merriam, C.H., 1890, General results, with special reference to the geographical and vertical 

distribution of species, in North American Fauna No. 3—Results of a biological survey of the San 

Francisco Mountain Region and Desert of the Little Colorado, Arizona: Washington, D.C., U.S. 

Department of Agriculture, Division of Ornithology and Mammalogy, Government Printing Office.

Miller, A.H., and Stebbins, R.C., 1964, The lives of desert animals in Joshua Tree National Monument:

Berkeley, Calif., University of California Press.

Miller, M.R., 1951, Some aspects of the life history of the Yucca night lizard, Xantusia vigilis: Copeia, 

v. 2, p. 114–120.

Montgomery, D.R., and Foufoula-Georgiou, E., 1993, Channel network source representation using 

digital elevation models: Water Resources Research, v. 29, p. 3,925–3,934.

Moore, I.D., Grayson, R.B., and Ladson, A.R., 1991, Digital terrain modeling: A review of 

hydrological, geomorphological, and biological applications, Hydrological Processes, v. 5, p. 3–30.

Nalder, I.A., and Wein, R.W., 1998, Spatial interpolation of climatic Normals—Test of a new method in 

the Canadian boreal forest: Agricultural and Forest Meteorology, v. 92, p. 211–225.

Newbold, T., 2010, Applications and limitations of museum data for conservation and ecology, with 

particular attention to species distribution models: Progress in Physical Geography, v. 34, p. 3–22.

Norris, K.S., 1953, The ecology of the desert iguana Dipsosaurus dorsalis: Ecology, v. 34, p. 265–287.

Parker, W.S., and Pianka, E.R., 1973, Notes on the ecology of the iguanid lizard, Sceloporus magister: 

Herpetologica, v. 29, p. 143–152.

Pearson, R.G. and Dawson, T.P., 2003, Predicting the impacts of climate change on the distribution of 

species: are bioclimatic envelope models useful? Global Ecology and Biogeography, v. 12, p. 

361–371.

Peters, A.J., Eve, M.D., Holt, E.H., and Whitford, W.G., 1997, Analysis of desert plant community 

growth patterns with high temporal resolution satellite spectra: Journal of Applied Ecology, v. 34, p. 

418–432.

Phillips, S.J., Anderson, R.P., and Schapire, R.E., 2006, Maximum entropy modeling of species 

geographic distributions: Ecological Modelling, v. 190, p. 231–259.

Phillips, S.J., and Dudik, M., 2008, Modeling of species distributions with Maxent—New extensions 

and a comprehensive evaluation: Ecography, v. 31, p. 161–175.

Phillips, S.J., Dudik, M., Elith, J., Graham, C.H., Lehmann, A., Leathwick, J., and Ferrier, S., 2009,

Sample selection bias and presence-only distribution models—Implications for background and 

pseudo-absence data: Ecological Applications, v. 19, p. 181–197.

Pianka, E.R., 1966, Convexity, desert lizards, and spatial heterogeneity: Ecology, v. 47, p. 1,055–1,059.

Pianka, E.R., 1967, On lizard species diversity—North American flatland deserts: Ecology, v. 48, p. 

333–351.

Posada, D., and Buckley, T., 2004, Model selection and model averaging in phylogenetics—Advantages 

of Akaike Information Criterion and Bayesian approaches over likelihood ratio tests: Systematic 

Biology, v. 53, p. 793–808.

Pratt, D., 1979, The thermal inertia approach to mapping of soil moisture and geology: Remote Sensing 

of Environment, v. 8, p. 151–158.

Qi, J., Chehbouni, A., Huete, A., Kerr, Y., and Sorooshian, S., 1994, A modified soil adjusted 

vegetation index: Remote Sensing of Environment, v. 48, p. 119–126.

R Core Development Team, 2011, R—A language and environment for statistical computing: Vienna, 

Austria, R Foundation for Statistical Computing, http://www.R-project.org/.

Rodgers, T.L., 1944, A new skink from the Sierra Nevada of California: Copeia, v. 2, p. 101–104.
91

Rosen, P.C., Holm, P.A., and Lowe, C.H., 1996, Ecology and status of shovelnose snakes (Chionactis) 

and leafnose snakes (Phyllorhynchus) at and near Organ Pipe Cactus National Monument, Arizona:

Final Report on the funded project Status and Natural History of U.S. Populations of the Organ Pipe 

Shovelnosed Snake, Arizona Game and Fish Department Heritage Grant IIPAM 192045, p. 68 + 17.

Roth, V.D., 1984, The spider family Homalonychidae (Arachnida, Araneae): American Museum 

Novitates, no. 2790, p. 11.

Rushton, S., Ormerod, S., and Kerby, G., 2004, New paradigms for modelling species distributions?:

Journal of Applied Ecology, v. 41, p. 193–200.

Sawyer, J.O., Keeler-Wolf, T., and Evens, J., 2009, A manual of California vegetation: California 

Native Plant Society, p. 1,300.

Shcheglovitova, M., and Anderson, R.P., 2013, Estimating optimal complexity for ecological niche 

models—A jackknife approach for species with small sample sizes: Ecological Modelling, v. 269, no. 

9–17.

Smallwood, K.S., and Morrison, M.L., 1999, Spatial scaling of pocket gopher (Geomyidae) density: 

Southwestern Naturalist, v. 44, p. 73–82.

Simley, J.D., and Carswell, Jr., W.J., 2009, The National Map—Hydrography: U.S. Geological Survey 

Fact Sheet 2009-3054, 4 p.

Stebbins, R.C., 1944, Some aspects of the ecology of the Iguanid Genus Uma: Ecological Monographs,

v. 14, p. 311–332.

Stebbins, R.C., 1948, New distributional records for Xantusia vigilis with observations on its habitat:

American Midland Naturalist, v. 39, p. 96–101.

Stebbins, R.C., 1951, Amphibians of western North America: Berkeley, Calif., University of California 

Press, 539 p.

Stebbins, R.C., 2003, A field guide to western reptiles and amphibians (3d ed.): Boston, Massachusetts, 

Houghton Mifflin, Boston, 532 p.

Stokland, J.N., Halvorsen, R., and Stoa, B., 2011, Species distribution modelling-Effect of design and 

sample size of pseudo-absence observations: Ecological Modelling, v. 222, p. 1–10.

Suarez-Gracida, C.G., and Alvarez-Castaneda, S.T., 2009, Physical and biological variables related to 

habitat preferences of rodents: Biodiversity and Conservation, v. 18, p. 2,779–2,797.

Sullivan, B.K., 2005, Bufo punctatus, in Lannoo, M.J., ed., Amphibina Declines—The Conservation 

Status of United States Species: Berkeley, Calif., University of California Press, p. 430-432..

Swetnam, T.W., and Lynch, A.M., 1993, Multicentury, regional-scale patterns of western spruce 

budworm outbreaks: Ecological Monographs, v. 63, p. 399–424.

Syfert, M.M., Smith, M.J., and Coomes, D.A., 2013, The effects of sampling bias and model complexity 

on the predictive performance of MaxEnt Species Distribution Models: PLoS ONE, v. 8, no. 2, 

e55158, doi:10.1371/journal.pone.0055158.

Tevis, L., 1966, Unsuccessful breeding by desert toads (Bufo punctatus) at limit of their ecological 

tolerance: Ecology, v. 47, p. 766.

Tracy, C., 1976, A model of the dynamic exchanges of water and energy between a terrestrial 

amphibian and its environment: Ecological Monographs, v. 46, p. 293–326.

Turner, R.M., 1982, Mohave desertscrub, in Brown, D.E., ed., Desert Plants: Biotic communities of the 

American Southwest—United States and Mexico, v. 4, p. 157–168.

Turner, F.B., and Randall, D.C., 1989, Net productivity by shrubs and winter annuals in southern 

Nevada: Journal of Arid Environments, v. 17, p. 23–36.

Turner, F.B., and Wauer, R.H., 1963, A survey of the herpetofauna of the Death Valley Area: Great 

Basin Naturalist, v. 23, p. 119–128.
92

Vandergast, A., Inman, R., Barr, K., Nussear, K., Esque, T., Hathaway, S., Wood, D., Medica, P., 

Breinholt, J., Stephen, C., Gottscho, A., Marks, S., Jennings, W., and Fisher, R., 2013, Evolutionary 

hotspots in the Mojave Desert: Diversity, v. 5, p. 293–319.

van Leeuwen, W.J.D., Davison, J.E., Casady, G.M., and Marsh, S.E., 2010, Phenological 

characterization of Desert Sky Island vegetation communities with remotely sensed and climate time 

series data: Remote Sensing, v. 2, p. 388–415.

Wallace, C.S.A., and Thomas, K.A., 2008, An annual plant growth proxy in the Mojave Desert using 

MODIS-EVI Data: Sensors, v. 8, p. 7,792–7,808.

Wang, P.X., Li, X.W., Sun, W., Li, X.M., Zhang, S.Y., and Liu, A.L., 2004, Soil thermal inertia 

estimation by combining afternoon and morning AVHRR data with a modified diurnal land surface 

temperature change model: in Proceedings of the IEEE International Geoscience and Remote Sensing 

Symposium: IGARSS '04, p. 4,299–4,301.

Weiss, J., Gutzler, D., Coonrod, J., and Dahm, C., 2004, Long-term vegetation monitoring with NDVI 

in a diverse semi-arid setting, central New Mexico, USA: Journal of Arid Environments, v. 58, p. 

249–272.

Wood, S.N., 2006, Generalized additive models—An introduction with R: Boca Raton, Florida, CRC 

Press, p. 392.

Wood, S.N., 2011, Fast stable restricted maximum likelihood and marginal likelihood estimation of 

semiparametric generalized linear models: Journal of the Royal Statistical Society—Series B 

(Statistical Methodology), v. 73, p. 3–36, doi: 10.1111/j.1467-9868.2010.00749.x.

Zeveloff, S.I., 1988, Mammals of the intermountain west: Salt Lake City, University of Utah Press, p. 

365.

Zweifel, R.G., and Lowe, C.H., 1966, The ecology of a population of Xantusia vigilis, the desert night 

lizard: American Museum Novitates, v. 2,247, p. 1–57.
Publishing support provided by the U.S. Geological Survey

Publishing Network, Tacoma Publishing Service Center

For more information concerning the research in this report, contact the

 Director, Western Ecological Research Center

U.S. Geological Survey

3020 State University Drive East

Sacramento, California 95819 

http://werc.usgs.gov/
Inman and others—Mapping Habitat for Multiple Species in the Desert Southwest Open-File Report 2014–1134 —

ISSN 2331-1258 (online)

http://dx.doi.org/10.3133/ofr20141134