Global loss of biodiversity is a well-known concern for conservationists and managers, but detailed spatial maps of local biodiversity for use by local managers are often lacking. We used a suite of existing species distribution models to calculate spatial variation in low-elevation species richness across Clark County, Nevada, USA, comprising much of the eastern Mojave Desert. We then used a macroecological model to estimate true latent low-elevation biodiversity across the county, correcting for potential taxonomic bias in the estimates of species richness. We found that species richness and biodiversity tended to be higher along the Muddy and Virgin rivers and in the Las Vegas valley. Biodiversity was positively associated with flat, rocky landforms, low elevation, late seasonal greenup, and lower differences between winter and summer temperatures. We present a brief example for local managers to apply the new publicly available low-elevation species richness and biodiversity spatial layers.

Study area

The Mojave Desert is the smallest North American desert and is generally arid. Precipitation patterns exist along a gradient, with the western Mojave drier and experiencing winter rains and the eastern Mojave wetter with both winter rain and snow and summer monsoon rains (Germano et al. 1994, Keeler-Wolf 2003, Pietrasiak et al. 2014). At low elevations in the eastern Mojave Desert, precipitation can average 11 cm per year (Abella et al. 2009).  Mojave Desert landforms are typically alluvial fans and basins, dominated by the creosote bush (Larrea tridentata) and white bursage (Ambrosia dumosa) shrub alliance at lower elevations, blackbrush (Coleogyne ramosissima) shrub communities at middle elevations, and pinyon-juniper (Pinus-Juniperus) forests at higher elevations (Keeler-Wolf 2003, Abella et al. 2009).

Clark County is the southernmost county in the U.S. state of Nevada and encompasses ~20,800 km², covering much of the eastern Mojave Desert (Figure 1). Approximately 88.9% is federally owned and managed with the remainder being state or private land. Elevation ranges from 137 m to 3,634 m (Abella et al. 2009). Riparian shrub communities (e.g., Salix spp., invasive Tamarix ramosissima) occur along the Virgin and Muddy rivers in the northeastern portion of the county. Approximately 12,800 km² are designated under some minimal type of conservation management, such as wilderness areas, federally protected lands, military lands, state parks, etc. (Clark County Desert Conservation Program, unpublished data). Included within the conservation acreage is the Boulder City Conservation Easement (BCCE), a 353 km² area managed by Clark County to protect habitat for Mojave desert tortoises and other species and which we used as an example of applying the species richness and biodiversity layers derived here (Figure 1). Anthropogenic activities on public lands can be extensive, including off-road vehicles, shooting, and illegal dumping. As one example of quantifying anthropogenic activities, from 1995–2006 an average of nine million people annually visited the Lake Mead National Recreation Area along the eastern border of Clark County (Figure 1; Abella et al. 2009).

Conceptual treatment of SDMs

As of 2018, the Clark County Desert Conservation Program had either commissioned or had access to SDMs for 55 plant, reptile, mammal, avian, and invertebrate species. The list of species chosen for deriving SDMs was not comprehensive of all species present within Clark County, but was reflective of species occurring at low elevations that were most subject to degradation from private land development and were covered under the multiple species habitat conservation plan. Therefore, the resultant estimates of richness should be interpreted as richness of low elevation species, hereafter ‘species richness’. Importantly, we then used the species richness layer and landscape data to generate a macroecological model generalized to all low-elevation species, and hereafter refer to this predictive layer as ‘biodiversity’ (Calabrese et al. 2013).

Data processing

Our first step in data processing was to classify all SDMs into one of three categories: suitable for use as provided (i.e., a continuous estimate of the relative probability of occurrence, scaled between zero and one), suitable for use after processing (e.g., species density estimates that can be scaled between zero and one), or unsuitable for this task (e.g., vector layers with qualitative classification of species’ occurrence or a probability raster that was only generated for a small portion of Clark County). For SDMs with density estimates that were all less than 1, the direct density estimate was retained as approximately equivalent to the probability of occurrence in that habitat cover type (i.e., low average observed densities = low probability of occurrence in a given plot). For density estimates that exceeded 1, densities were rescaled to a zero-to-one scale by dividing the observed value by the maximum value rounded up to the nearest tenth place to ensure all values were less than 1. For ‘suitable after processing’ SDMs we assumed that within-species variation in density was a proxy variable for variation in probability of occurrence.

Simple spatial re-projection, raster snapping, or raster resampling were applied where necessary to standardize all of the SDMs to the same resolution (i.e., 250m x 250m), raster origin coordinates, and coordinate system (i.e., NAD83 UTM Z11N). All NoData cells were assigned a value of 0 prior to stacking. Some SDMs were resampled to a higher resolution (250m) than that of the original SDM (1km). In practice, this meant that a grid cell that was previously 1km x 1km with a single raster value now was comprised of four grid cells, each 250m x 250m, and each with the same original raster value. In terms of spatial accuracy, this falsely implies a higher resolution at the level of the individual SDM, but we decided it was acceptable to do at the level of the SSDM because it allowed them to accommodate the resolution of the 250m SDMs. This approach maximized inclusion of all information from the highest-resolution SDMs. Some of the avian SDMs were generated at an even higher resolution than 250m x 250m. These were resampled to a 250m x 250m resolution taking the maximum value within the window to reflect the highest probability that a species was observed within the larger window. After processing, all SDMs were a continuous, probabilistic probability of species occurrence, or p-SDM (Calabrese et al. 2013).

Stacking p-SDMs and macroecological model

Next we ‘stacked’ all of the p-SDMs together via summation over overlapping grid cells, yielding a continuous index of species richness (a p-SSDM; Figure 2), a method that has been found to be robust and accurate when compared with validation data (Calabrese et al. 2013, Zellmer et al. 2019, Zurrell et al. 2020). The stacked p-SSDM is useful in that it retains information on the input species, but there is a potential for bias in the p-SSDM with regards to true biodiversity depending on the level of latent bias in which species were chosen for creation of the individual SDMs (Distler et al. 2015). For example, if only reptile species were used for the species richness model, that model may not represent biodiversity of plants. To resolve this we conducted a second analysis, whereby we built a macroecological model that sought general rules in the eastern Mojave Desert driving species assemblages (Calabrese et al. 2013, Distler et al. 2015). We used linear regression to model the p-SSDM species index as a function of environmental predictors. The environmental predictor variables were chosen from the six most important variables identified by Inman et al. (2014) in a study on SDMs for 15 species across the Mojave Ecoregion. The five predictor variables used in this macroecological model were elevation, seasonal thermal difference (i.e., the difference between winter and summer surface temperatures), surface texture (i.e., an index of sandy versus rocky surface), topographic position (i.e., high values equal valley bottom while low values equal mountain or ridge top), and seasonal greenness timing (i.e., the average date of maximum vegetation greenness; Table 1). Winter precipitation was a sixth environmental predictor variable initially included. However, winter precipitation was highly correlated with elevation (r = 0.863) and was excluded from further analysis because previous work has shown elevation to be a primary driver of regional variation in biodiversity (Mateo et al. 2012).

We then generated a systematic grid of sample points 2.25 km apart from each other to sample the p-SSDM species richness raster and the environmental predictor variable rasters while ensuring even coverage across the county and minimizing spatial autocorrelation. After removing points from the sample grid that fell within ‘disturbed’ areas (i.e., graded or paved construction as determined from aerial imagery analysis, L. Bice unpublished data) there were 3,827 sample points. We then used linear regression in Program R v.3.5 (R Core Team 2018) to derive the macroecological statistical model. Finally, we used the coefficients from the statistical model in Raster Calculator in ArcGIS 10.4 (ESRI, Inc.) to create the continuous biodiversity predictive surface. All environmental raster layers were available at a resolution of 1 km, so the output macroecological layer is at this resolution as well.