The study area encompasses 132,441 km², including the Level III Mojave Desert ecoregion and some adjacent ecoregions (Omernik and Griffith, 2014) in Utah, Arizona, Nevada, and California, USA ( Figure 1 ). The study area perimeter was determined using the probability threshold of 0.3 or higher from a previous species distribution model for Joshua trees (Godsoe et al., 2009) and based on field observation. Baseline survey elevations were from 400 m to 2,200 m to encompass the range of Joshua trees but exclude extraneous search areas because of the size of the study area. The elevational range across the study region was −86 m in Death Valley, California, to 3,632 m at the summit of Charleston Peak, Clark County, Nevada.

Publicly available Google Earth satellite imagery was used to identify presence or absence of adult Joshua trees of both species (Esque et al., 2020a). We created a structured data set for both species with a repeatable protocol (Isaac et al., 2020). Rather than a stratified sampling design, we aimed to sample 100% of the study area at 500 m resolution. Using the Fishnet toolbox in ArcMap 10.5 we created a grid of 811,900 500 m × 500 m cells in the USA Contiguous Albers Equal Area Conic USGS coordinate system (SR-ORG:7301). Observers inspected 672,043 of the survey cells during our initial mapping phase, encompassing 82.7% of the gridded study area. The remainder of the gridded study area, or 139,857 cells, occurred in south-central Nevada near US Department of Energy and Defense facilities (i.e., National Nuclear Testing Facility and Nellis Air Force Base) where we found obscured imagery presumably due to national security concerns. Thus, presence/absence determinations were precluded in this area ( Figure 2 ). As an alternative, we modeled Joshua tree habitat in this region using SDM algorithms based on the true presence and absences for the rest of both species’ ranges (see below).

Visual scans of adult Joshua tree presence and absence were conducted at a standardized eye elevation of ~250 m (i.e., the altitude above the land surface) with a target search time of roughly 45 s per cell. Observers either placed waypoints within each cell at the location of a prominent Joshua tree or designated absence. Surveys of satellite imagery were mostly limited to adult Joshua trees (usually branched trees >2 m height) because initial field surveys indicated that smaller trees are usually not distinguishable, and those <1 m tall are undetectable (Esque et al., 2021).

The quality of Google Earth imagery was variable, but usable, across the gridded study area except for the region of obscured imagery described above (see Completing coverage for Joshua tree distributions in obscured area using SDMs). The most recent imagery (2021) was evaluated first, but if presence/absence was not readily assigned because of image quality, the eye altitude or time frame within the historical satellite imagery (2003 to 2021) was varied to try and detect Joshua trees. Remaining questionable cells were re-evaluated using secondary satellite surveys (see below).

We further refined the habitat map by correcting errors using field validation, secondary satellite searches, empirical point data from co-authors’ unpublished datasets, points provided by staff at Nellis Air Force Base, filtered research grade iNaturalist observations (GBIF, 2023), and some SEInet observations (Esque et al., 2020a; SEINet, 2023). Field validation included 15,001 km of driving along both paved and dirt roads during twenty site visits throughout the grid system (29,050 cells; Figure 2 ). Where site visits were not possible, the most experienced observers re-evaluated cells using secondary satellite searches (from any year available, with the best imagery) in the following three scenarios: 1) cells determined to be without Joshua trees but adjacent to cells having Joshua trees present; 2) cells having Joshua trees present but surrounded by cells without Joshua trees; and 3) cells in areas known to cause confusion because of other issues (e.g., plants, rocky outcrops, fire scars, the urban/wildland interface). Data from iNaturalist and SEInet that were inconsistent with our database were also field validated.

We derived 18 environmental variables to serve as covariates in the species distribution models (SDMs), which together characterize climate, topography, vegetation (e.g., NDVI variables), and soil surface properties for the study region ( Table 1 ). Precipitation and temperature layers were created using ClimateNA v. 7.3 (Wang et al., 2016), which downscales PRISM data (Daly et al., 2008) and corrects for elevational variation. Our contemporary climate analyses included data from the 30-year period between 1980–2010. Satellite metrics incorporating plant canopy and soil surface data from the moderate-resolution imaging spectroradiometer (MODIS) satellite were averaged across 17 y (2003–2020) to represent a norm for the study region (NDVI amplitude and maximum – USGS eMODIS Remote Sensing Phenology, https://doi.org//10.5066/F7PC30G1). A layer representing soil surface texture was downloaded and mosaicked from the SoilGrids 2.0 web portal (Poggio et al., 2021). All topographic metrics were calculated by aggregating a 30 m digital elevation model to the 500 m × 500 m resolution used for modeling (National Elevation Dataset, http://ned.usgs.gov/).

Climatographs were generated for Y. brevifolia and Y. jaegeriana from the final corrected presence grid cells using gaussian kernel density estimates from cell values for each environmental variable. The climatographs were used to compare the climate between adult Joshua tree species’ occupied habitat and can be compared with partial response curves resulting from SDMs.

We used an ensemble modeling approach to create SDMs for Y. brevifolia, Y. jaegeriana, and a rangewide model with both species’ datasets (hereafter, rangewide model). Presence and absence points were assigned to each species’ range based on Rowlands (1978); Lenz (2007), and author’s unpublished genetic data (CIS). The species are non-overlapping except in the hybrid zone which was obscured by poor imagery; in that area the data for both species were used in combination for modeling (i.e., species specific points were not designated). Next, we fit SDMs to the individual species occurrences separately, along with a rangewide SDM featuring all points (i.e., from both species). We used a custom script in R version 4.1.3 (R Core Team, 2022) to implement model cross-validation, model averaging, parallel processing, and partial response curves for model terms. Our ensemble modeling approach included two algorithms: generalized additive models (GAM; R package “mgcv” version 1.8-22; Wood, 2017) and random forests (RF; R package “ranger” version 0.12.1; Wright and Ziegler, 2017). Both algorithms have been consistently strong performers among SDM algorithms (Franklin, 2010). However, due to consistently better cross-validation performance, we only fit RF models for the rangewide dataset, whereas both GAM and RF were compared for the individual species models. All GAM models were fit with restricted maximum likelihood (REML) and an extra penalty allowing smooth terms to be penalized to zero (“gam” option select=TRUE in “mgcv” package) to aid model selection. Random forest models were fit with 1,500 random trees and “ranger” package defaults for the regression model (Wright and Ziegler, 2017).

For each algorithm and dataset, we considered eight candidate models that included 10 uncorrelated terms and another 8 terms used such that we avoided multicollinearity as environmental predictors ( Table S1 ). All models included the same 10 uncorrelated terms (NDVI amplitude – AMPn, NDVI maximum – MAXn, HCL, precipitation seasonality – PCV, precipitation ratio – Pratio, Sand, Slope, topographic position – TPI, temperature seasonality – TSD, and temperature range – Trange; Table 1 ), but varied with respect to eight additional terms to avoid multicollinearity (Climate moisture deficit – CMD, Annual heat/moisture index – AHM, mean annual precipitation – MAP, mean annual temperature – MAT, summer precipitation – SP, summer maximum temperature – STMX, winter precipitation – WP, and winter temperature minimum – WTMN). This set of models allowed us to contrast seasonal and annual climate variables which, though correlated, may relate to different critical life stages or phenological processes for the different Joshua tree species.

To account for potential bias due to spatial aggregation (Veloz, 2009), and for computational efficiency, we rasterized presence and absence points to the modeling resolution (500 m × 500 m) and applied a spatial thinning procedure in which a maximum of 50 points could be randomly sampled from any 5 km² area (Fourcade et al., 2014). This procedure reduced data density without changing the spatial arrangement of points. Following this procedure, the models included 128,438 points for Y. brevifolia (including presences and absences), 114,705 points for Y. jaegeriana, and 198,305 points for the rangewide model (note: rangewide points are not additive because of potential species overlap in a large area). This procedure was repeated five times for each species and for the rangewide dataset. Next, for each set of thinned points, we applied a 5-fold cross-validation procedure to split the data into training folds (used for model fitting) and testing folds (used for model evaluation). Overall, this process resulted in 25 total cross-validation runs used to evaluate models (5 sets of randomly thinned points × 5-fold cross-validation). It has been shown that conventional random cross-validation may underestimate model error and/or transferability, particularly when applied to spatially structured data (Bahn and McGill, 2012; Wenger and Olden, 2012). For this reason, we used a spatial cross-validation approach developed by Valavi et al. (2019) and implemented in the R package “blockCV” v2.1.4. Using the packages’ “spatial blocking feature”, we split data into random training and testing blocks that were both geographically separated and optimized to contain relatively equal numbers of presences and absences. With this procedure, we intended to reduce overfitting to the training data while identifying models that were generalizable and performed well across the study extent.

To measure model performance of the cross-validated models, we considered several metrics of model prediction accuracy including AUC (i.e., the area under the receiver operating characteristic; Fielding and Bell, 1997), the Boyce Index (Hirzel et al., 2006), and the True Skill Statistic (TSS; Allouche et al., 2006). For GAM, we also calculated each model’s average AIC (with each model being fit to the same subsets of data during cross validation) to help identify well-performing, parsimonious models. These metrics were averaged across cross-validation runs for each model to obtain an overall estimate of performance. Moreover, for each cross-validation run, we also generated model predictions as raster grids for the study extent constituting the predicted habitat suitability probabilities. We then created a final prediction raster for each model as the weighted average of the 25 individual cross-validation predictions based on the TSS scores (such that better performing model predictions from the cross-validation runs featured more heavily in the final raster prediction surface for each model). Similarly, we created overall algorithm ensemble predictions as the weighted average of raster predictions from individual models, again based on TSS scores. For each ensemble model, we also calculated model calibration curves by binning the predicted habitat probalities into 5 probability classes of width 0.2 (e.g., 0–0.2, 0.2–0.4, etc.) and comparing the midpoints of each bin with the observed frequency of presences within that subset of data. Calibration curves are useful in illustrating how well habitat probabilities correspond to the actual frequency of presence observations.

To aid model interpretation, we derived relative importance values for each predictor present in the candidate models for each algorithm. Relative importance for predictors in random forest models were based on the mean decrease in accuracy from permutations leaving out each term, while for GAM, relative performance was based on the parameter’s chi-square statistics. We also derived partial variable response curves for each of the top nine predictors present in the candidate models for each species, as well as the rangewide dataset. These curves indicate the shape and direction of relationships between predictors and habitat probability values. For GAM, response curve functions for predictors were also averaged across all the models in which each predictor occurred to create a model-averaged response curve for each predictor, which we overlay on the individual curves from candidate models. For random forest models, we calculated partial dependence curves using the R package “pdp” (Greenwell, 2017).

We define the modeled habitat for Joshua trees as areas where SDM-predicted habitat probability values were above the threshold that maximized the sum of model sensitivities (true positive rate) and specificities (true negative rate; Liu et al., 2005). SDM ensemble habitat layers were thresholded separately for each species and for the combined dataset. Potential habitat for Joshua trees includes locations with modeled habitat where empirical Joshua tree presences were not detected (e.g., where there were no presence points from surveys). We quantified the area of habitat resulting from empirical distribution data versus the modeled and potential habitat generated by SDMs.

To compare the modeled habitats of Y. brevifolia and Y. jaegeriana, we calculated Schoener’s D statistic and the I statistic defined in Warren et al. (2009), both measures of environmental overlap which range from 0 (no overlap in covariates) to 1 (complete overlap). As a second measure, we calculated AUC values for two scenarios: (1) Y. brevifolia’s habitat probabilities (from surveyed presences) within the SDM predicted for Y. jaegeriana and projected across the full range of both species; and (2) Y. jaegeriana’s habitat probabilities (from surveyed presences) within the SDM predicted for Y. brevifolia and projected across the full range. These calculations illustrate how well each species’ individual SDMs predict occupied habitat for the other species, as well as whether any such relationships are spatially symmetric. Finally, to contrast the environmental tolerances of each species, we overlaid the response curves for environmental predictors that were among the top nine predictors in the individual SDMs for either species (overlaid separately for GAM and random forest). These response curve overlays contrast the relationships between each species’ modeled habitat by each environmental covariate.

Field validation of 29,050 grid cells resulted in status changes for 3,703 of them (12.8%) (i.e., changes from absence to presence, or vice versa; Table S2 ). Experienced observers re-evaluated 76,578 cells using secondary satellite searches, changed the status of 36,857 cells (48.1% of re-surveyed cells, or 5.5% of all cells surveyed). The majority of these cells were changed from absence to presence (72.4% of changed cells; Table S2 ).

Based on the gridded image surveys and final corrections, Yucca jaegeriana currently occupies more than 16,683 km² of habitat (excluding the area of obscured imagery) in the eastern Mojave Basin and Range ecoregion, and conterminous ecotones of Sonoran Basin and Range, Arizona/New Mexico Mountains, and Southern Basin and Range ecoregions ( Figure 2 ; Omernik and Griffith, 2014). Nevada has the greatest amount of occupied Y. jaegeriana habitat (7,137 km²), followed by Arizona, California, and Utah ( Table 2 ). Yucca jaegeriana generally occurs in 11 large and homogeneous populations, along with a few smaller peripheral populations. Unoccupied cells within these populations are rare ( Figure 2 ). Occupied cells outside the well-defined populations are infrequent and mostly within 0.4 to 2.0 km of population edges, with only two individual cells isolated by >3 km. While uncommon, such isolates occur throughout the range of Y. jaegeriana (77 out of 102,274, 0.25 km² grid cells, Figure 2 ). The elevational limits of Y. jaegeriana are 392 m near Alamo State Park, AZ and 2,319 m in the Sheep Range, NV. The latitudinal limits are near Aguila, AZ in the south and Dry Lake Valley, near Caliente, NV in the north. Longitudinal limits are near Hawkins, AZ in the east, and the Avawats Mountains on Ft. Irwin Military Base, CA in the west.

Yucca brevifolia occupies 15,955 km² of habitat (excluding the area of obscured imagery) in the western Mojave Basin and Range ecoregion and adjacent Sierra Nevada and Southern California Mountain ecoregions (Omernik and Griffith, 2014; Figure 2 ). California has the most occupied habitat for Y. brevifolia followed by Nevada. Native populations do not occur in Utah or Arizona ( Table 2 ). Yucca brevifolia has an extensive population on the southwestern edge of its range with many narrow pinch points throughout. Besides another large stand in the northwest periphery of its range that is contiguous with the region of obscured imagery, there are several smaller isolated populations of Y. brevifolia southward and eastward throughout California ( Figure 2 ). Occupied cells that are isolated outside of large Joshua tree stands are infrequent with 170 cells isolated by <3 km and only 2 occupied cells isolated by >3 km from larger populations. In the southwestern range of Y. brevifolia, unoccupied clusters of cells are more frequently interspersed within occupied populations, creating a more diffuse pattern of occupancy in comparison with Y. jaegeriana ( Figure 2 ). The elevational limit of Y. brevifolia is 600 m near Cantil, CA and the upper is 2,605 m on Maturango Peak, CA on the Naval Air Weapons Station – China Lake. The latitudinal limits are from near Indio, CA in the south and just south of Tonopah, NV in the north, and longitudinal limits are near Twentynine Palms in the east and the western limit is currently at the junction of Orwin Way Road and the Quail Canyon Motocross Road in Los Angeles, Co, CA.

Climatographs were used to display the frequency distributions of environmental variables within occupied Joshua tree habitat ( Figure 3 ). While annual heat/moisture index (AHM) and winter precipitation (WP) had similar distributions between Y. brevifolia and Y. jaegeriana habitats, other variables – particularly summer precipitation (SP), precipitation ratio (Pratio), temperature standard deviation (TSD), and precipitation coefficient of variation (PCV) – showed marked discrepancies ( Figure 3 ). Yucca jaegeriana experiences a greater range of temperature variation (higher TSD) than Y. brevifolia throughout the year, but a more even precipitation distribution (i.e., lower PCV). Interestingly, due to differences in the precipitation regime, Y. jaegeriana and Y. brevifolia share a similar profile of annual aridity (AHM), despite Y. jaegeriana experiencing higher mean annual temperatures (MAT). For Y. brevifolia, annual precipitation is largely restricted to the winter months due to its western position in the Mojave precipitation gradient. However, Y. jaegeriana receives more bi-modal precipitation benefitting from tropical monsoonal storms in its easterly position (Hereford et al., 2006; Figure 3 ), and consequently higher annual precipitation. However, the climatic moisture deficit model (based on Dobrowski et al., 2013), which incorporates additional variables beyond MAP and MAT (e.g., downward shortwave radiation and wind velocity), suggested that Y. jaegeriana has a greater overall moisture deficit than Y. brevifolia in portions of its range.

SDMs predicted the Joshua tree habitats with high model performance ( Tables 3 , S3 ). Most AUC values from spatial cross-validation were greater than 0.8, and scores for the algorithm ensemble models were higher for RF than for GAM ( Table 3 ). The RF model for the rangewide map performed similarly compared to the separate species models, with an AUC = 0.868 and R ² = 0.809 ( Table 3 , Figure 4 ). Given that spatial cross-validation is a more stringent approach for model evaluation than traditional random cross-validation (i.e., models are scored based on their ability to predict into distinct geographic areas), these metrics suggest that the SDMs are generalizable rather than overfit and accurate in their predictions across the full distribution. Inter-model standard deviations were generally low (<0.2 in habitat probabilities) but more pronounced in areas with fewer data points, such as the northern part of the range where empirical coverage is low. Standard deviations for the individual species models were also higher where we extrapolated predictions into the range of the other species to facilitate habitat comparisons ( Figures 5 , 6 ), but these areas were typically of lower predicted habitat suitability, and the rangewide ensemble SDM is not affected by this issue. Model calibration curves for the ensemble models indicated that, collectively, the models tended to predict somewhat higher probabilities than observed presences at the low end of values (<0.4), and somewhat lower probabilities than expected at the higher end ( Figure S1 ). However, the estimates were close enough to the observed values that they would not result in mispredictions when the model threshold is applied to determine a suitability cutoff.

Climate variables were typically more important predictors in the SDMs than the vegetation or topography variables for each species. For Y. jaegeriana, the GAM model with the best (lowest) AIC was model 3, while for Y. brevifolia, the GAM with the lowest AIC was model 7 ( Supplementary Table S3 ), and models 3 and 7 were the top two models for both Joshua tree species on this metric ( Table S3 ). The environmental variables for these models only differed by including the climate moisture deficit (CMD for model 3), and summer maximum temperature (STMX for model 7). Despite the large differences in AIC, GAM models were largely consistent in predictive performance based on other metrics, with AUC ranging from 0.796 to 0.798 for Y. jaegeriana, and from 0.852 to 0.859 for Y. brevifolia ( Table S3 ) among the eight models. No single GAM model scored highest across all performance metrics.

Among RF models, we observed less variation in predictive performance across models than for GAM ( Table S3 ). For Y. jaegeriana, models 2 and 1 had the highest AUC scores, while models 2 and 5 had the highest TSS scores ( Table S3 ). Both models 1 and 2 contained a measure of annual aridity (AHM); similarly, model 5 reflected annual precipitation and temperature variables (MAP, MAT) but not their seasonal components. However, the models varied in AUC by only 0.004 (0.859 versus 0.855). For Y. brevifolia, models 4 and 3 had the highest AUC and TSS scores respectively ( Table S3 ). These models only differed by including either seasonal precipitation (SP, WP, Model 3) or annual precipitation (MAP, Model 4). However, as with Y. jaegeriana, AUC differed by only 0.004 between the highest and lowest performing RF models for this species (0.878 versus 0.874), suggesting similar predictive performance. The rangewide, multi-species RF models also showed consistently high predictive performance, with Models 7 and 8 showing the highest AUC scores ( Table S3 ). These models both contained seasonal precipitation variables but varied in representing seasonal (Model 7) vs. annual (Model 8) temperature averages.

For Y. jaegeriana, AHM and MAP were the strongest predictors of Joshua tree habitat averaged across algorithms, followed by SP, MAT, Pratio, PCV, and TSD ( Table 4 ). For Y. brevifolia, MAT was the strongest predictor averaged across GAM and RF followed by TSD, PCV, SP, AHM, WTMN, and Pratio ( Table 4 ). The best eight environmental predictors averaged between RF and GAM differed in rank but were largely shared between the two species, differing only with respect to STMX (for Y. jaegeriana) and MAP (for Y. brevifolia). For the rangewide RF model, AHM was the strongest predictor, followed by MAP, MAT, SP, WP, STMX, Pratio, PCV, TSD, and WTMN ( Table 4 ).

We used partial response curves to illustrate relationships between Joshua tree modeled habitat and individual environmental variables across gradients in the landscape. Both Y. brevifolia and Y. jaegeriana have similar peaks in probability of occurrence for the annual heat/moisture index (AHM), climatic moisture deficit (CMD), and mean annual precipitation (MAP) for the GAM and RF models ( Figure 7 ). In contrast, Y. jaegeriana are predicted to occur at somewhat higher temperatures than Y. brevifolia for both RF and GAM algorithms ( Figure 7 ). Summer precipitation (SP) increased habitat probabilities at higher values for Y. jaegeriana in the GAM and RF models, and this relationship was also reflected in higher probability values with increasing ratio of summer to winter precipitation (Pratio). Both the RF and GAM response curves indicated that habitat probabilities for Y. brevifolia may also increase in response to a higher ratio of summer precipitation, which occurs in Y. brevifolia habitat in the vicinity of the hybrid zone and in the most northerly parts of the species range in NV. Similarly, while Y. jaegeriana typically occurs in areas with more bimodal precipitation and hence lower precipitation seasonality (PCV), response curves suggested that higher PCV values (i.e., more winter-dominated precipitation) increased habitat probability in the GAM model for Y. jaegeriana, likely reflecting the large amount of potential habitat occurring in the western Mojave that is often congruent with Y. brevifolia habitat ( Figure 5 ). This is also illustrated in the Fort Irwin area, where the GAM model predicted higher habitat suitability values than the RF model (see Figure 5 ). Higher levels of temperature variability were more influential in the SDMs for Y. jaegeriana for both RF and GAM ( Figure 7 ). Full partial response curves from GAM and RF models for each species are available in Supplementary Figures S2–S6 .

Niche similarity statistics used on SDMs suggested considerable overlap between the two species, with D = 0.629 and I = 0.882 for the multi-algorithm ensemble predictions. However, AUC comparisons suggested that the niche similarity was asymmetric: while the Y. jaegeriana ensemble SDM predicted Y. brevifolia presences with an AUC = 0.764, the Y. brevifolia ensemble SDM predicted Y. jaegeriana presences with an AUC = 0.667. We note that SDMs extrapolated beyond the known range of a species are prone to uncertainty: we projected each species SDM into the range of the other species to illustrate where the habitat may be similar, but this does not imply occupancy of both species in these areas, and these observations are corroborated by the habitat mapping reported here.

Spatial comparisons of the ensemble SDMs for each species further illustrate that potential habitat (modeled habitat where presences were not detected in surveys) for Y. jaegeriana extends farther into the habitat of Y. brevifolia than vice versa ( Figure 8 ). In fact, potential habitat for Y. jaegeriana extends throughout the entire range of Y. brevifolia, while potential habitat of Y. brevifolia only occurs in small patches beyond its extant range ( Figures 8B, C ), in addition to the area of obscured imagery where both species share potential habitat. Overall, potential habitat derived from the ensemble SDM for Y. brevifolia covers approximately 16,400 km² and is of similar size to the area of occupied habitat ( Figure 8B ). However, potential habitat had a median probability value of 0.624, lower than the median of 0.726 for grid cells in occupied habitat (e.g., those overlapping with surveyed presences). Potential habitat for Y. brevifolia also overlaps with approximately 1,834 km² of occupied Y. jaegariana habitat. Potential habitat for Y. jaegeriana covers approximately 27,638 km² and is substantially larger than the occupied habitat for this species ( Figure 8C ). As observed for Y. brevifolia, potential habitat for Y. jaegeriana may be of lower current suitability, with a median habitat probability of 0.61 versus a median of 0.73 for occupied habitat. Approximately 5,410 km² of Y. jaegeriana’s potential habitat is currently occupied by Y. brevifolia, mirroring the asymmetric niche similarity measures noted above. Overall, the two species share approximately 9,669 km² of potential habitat that is not currently occupied by either species. Much of this area of shared potential habitat falls within the area of obscured imagery ( Figure 8 ); hence, a caveat to the calculations presented here is that the amount of occupied habitat in the obscured imagery area is unknown.

Potential habitat predicted for both species tended to occur at higher elevations than currently occupied habitat, suggesting that SDMs did not typically extrapolate into warmer areas. For example, potential habitat for Y. brevifolia was at significantly higher elevation than occupied habitat (mean of 1459 m versus 1325 m, respectively; t = 71.07, P < 0.0001), received more annual precipitation (MAP of 206 mm versus 189 mm; t = 57.52, P < 0,0001), and had lower annual temperatures (MAT of 13.7°C versus 14.4°C; t = 61.62; P < 0.001). Similarly, potential habitat for Y. jaegeriana was higher in elevation (average of 1408 m versus 1078 m; t = 230.99, P < 0.0001) and cooler (MAT of 14.1°C versus 17.2°C; t = −295, P < 0.0001) than occupied habitat. However, precipitation was lower in Y. jaegeriana’s potential habitat (MAP of 190 mm versus 210 mm in occupied habitat; t = −86.52, P < 0.0001), reflecting a westward shift into areas receiving less summer precipitation that are more typical of Y. brevifolia habitat.

While higher elevation and lower temperatures in potential habitat suggest these areas could serve as refugia as the climate warms, accessibility to current populations could be a major constraint. For Y. jaegeriana, potential habitat that is contiguous with actual habitat is rarely more than a few kilometers deep around Y. jaegeriana stands ( Figure 8 ). Outlying patches and diffuse cells of potential habitat are likely inaccessible for colonization under present climatic conditions because the unoccupied intervening areas do not support Joshua trees and Joshua trees seem to be dispersal limited (Vander Wall et al., 2006). While the pattern is similar for most of Y. brevifolia’s range, there are some large potential habitat patches in the southwest of the range that are largely contiguous with occupied habitat.

Our species distribution models are based on the most comprehensive presence and absence data yet available for Joshua trees, enabling us to evaluate the environmental variables underlying each species’ distribution. Thus, we present the first highly accurate and nearly complete empirically derived range maps for Yucca brevifolia and Y. jaegeriana. By coupling Google Earth imagery to detect adult Joshua trees with extensive field validation, we identified presences on 0.25 km² grids across the ranges of both species. While presence data were lacking for one area of obscured imagery where national defense is a priority for the Department of Energy and the Department of Defense, we used our high-resolution presence and absence data to fit SDMs across the full range, resulting in completed rangewide distribution maps.

Identifying large, multi-stem Joshua trees using remotely sensed Google Earth data is relatively straightforward, especially in the absence of other large non-target succulent species that are similar in physiognomy. However, Joshua trees can be challenging to distinguish remotely from other large species including Mojave yucca (Y. schidigera Roezl ex Ortgies), soaptree yucca (Y. elata Engelm.), banana yucca (Y. baccata Torr.), giant saguaro cactus (Carnegiea gigantea (Engelm.) Britton & Rose), desert almond (Prunus fasciculata (Torr.) A. Gray); sugar sumac (Rhus ovata S. Watson), creosote bush (Larrea tridentata Coville), or trees such as pinyon pines (e.g., Pinus monophyla Torr. & Frém.) and junipers (e.g., Juniperus osteosperma (Torr.) Little). Clonal Joshua tree forms can also be difficult to observe remotely where they are prevalent, e.g., on the western margin of Y. brevifolia distributions. Areas scarred by wildfires can have exceptionally low densities of Joshua trees, which also makes them difficult to detect. Joshua trees that were within natural distributions, but recognizably in cultivation (i.e., planted specimens within fences, urban settings) were not included as presences. Despite these difficulties, extensive field validations of imagery-based presence and absence classifications indicated the method to be highly effective, with 87% of over 29,000 field surveyed grid cells requiring no change ( Table S2 ). Secondary imagery-based validations of over 76,000 grid cells by experienced observers further refined our presence and absence data, resulting in reclassifications of 5.5% of grid cells and likely increased accuracy in the most challenging areas of imagery ( Table S2 ).

Yucca jaegeriana habitat occurs in discrete well-defined and homogeneously occupied populations of the eastern Mojave Desert and neighboring ecoregions, while Y. brevifolia populations have a sprawling distribution of considerably less homogeneous patches along the western boundary of the Mojave Desert. Furthermore, empirical demographic measurements at the leading edges of Joshua tree distributions indicate that small founder trees occurring there extend less than a kilometer from the edges of established Joshua tree stands. This is consistent with observations during demographic transect surveys, in which young Joshua trees are well represented within some Joshua tree stands and taper off within a few 100 m near the stand edges (Smith, et al., unpublished data). These distributional patterns are generally consistent with previous Joshua tree range maps (Rowlands, 1978; Godsoe et al., 2009; Cole et al., 2011; Smith et al., 2011; Wilkening et al., 2020); however, the habitat we defined has 23.8% less areal coverage than the distribution identified in the most recently published distribution map (see Supplementary Table S4 ; Wilkening et al., 2020). Within populations we found absences of Y. jaegeriana only in small, widely scattered locations; in contrast, Y. brevifolia populations were more diffuse owing to a greater number of absences dispersed throughout the populations, especially in the southwest of the range. Based on satellite imagery some of those absences within Y. brevifolia populations appear to be the result from urban development, fire, and other cumulative disturbances.

The upper and lower elevational, latitudinal, and longitudinal limits were revised for both species where appropriate (see Figure 2 ). Differences we noted from previous work mostly represent the higher resolution of imagery now available and the benefits of widespread ground searches, rather than recent changes in species range limits. The northern limit for Y. jaegeriana range was increased by a few kilometers but the eastern limit that we interpreted from Rowlands (1978) was reduced by a few kilometers (P. Rowlands – pers. comm.), while the southern and western limits remain unchanged (Rowlands, 1978). The northern limit of Y. brevifolia was increased by 35 km, and the western limit was not previously well-defined, and misrepresented by a herbarium record with an erroneous locality, but is currently at the junction of Orwin Way road and the Quail Canyon Motocross road in Los Angeles, Co., CA. The eastern limit remains in Tikaboo Valley, NV (verified genetically near the hybrid zone – Starr et al., 2013) and the southern limit also remains the same as identified by Rowlands (1978).

Applying SDMs to our Joshua tree presence and absence data facilitates understanding the ecological correlates of Joshua tree distributions. Clarifying the roles of abiotic and biotic ecological factors are key to understanding species distributions (Lexer and Fay, 2005), and how Joshua trees may interact with future disturbances such as climate change. In this study, climate variables were the most important environmental correlates of Joshua tree distributions compared with remotely sensed or topographic variables; however, we did not include some important biotic variables in our analyses such as pollinator biology or intraspecific genomic variation also important to Joshua tree distributional patterns (Godsoe et al., 2008; Smith et al., 2008; Royer et al., 2016; Royer et al., 2020).

Several variables highlight differences between Y. jaegeriana and Y. brevifolia in the climatographs generated from survey-based presence data, and functional response curves generated from SDMs. Variables with the greatest contrast between species were the precipitation coefficient of variation, the precipitation ratio, and temperature standard deviation ( Figure 3 ). Climatographs showed that the overall frequency distributions of both species were similar for mean annual precipitation and winter precipitation ( Figure 3 ), however, only Y. jaegeriana occupies regions receiving regular summer precipitation. This influence is also correlated with the distribution of the species in relation to the precipitation coefficient of variation. Yucca jaegeriana also experiences more variable temperature patterns than Y. brevifolia ( Figure 3 ). For example, mean annual temperature had a similar amplitude between species but was clearly skewed toward higher temperatures for Y. jaegeriana and lower temperatures for Y. brevifolia which is similar to functional response curves for SDMs ( Figure 7 ). The standard deviation for monthly mean temperatures (TSD) was also higher for Y. jaegeriana, indicating the species experiences both warmer and more variable temperatures.

Higher temperatures experienced by Y. jaegeriana are ameliorated by greater summer precipitation (SP), such that the overall aridity profile (e.g., AHM, Figure 3 ) is similar between species. Patterns in seasonality and amount of precipitation experienced by Y. jaegeriana contrast with those of Y. brevifolia, particularly for the most easterly populations. This contrast is the result of a regional precipitation pattern of greater winter precipitation in the west, a more even precipitation pattern central to the range and near the hybrid zone, and bi-modal precipitation in the east (Hereford et al., 2006). Central and northern populations of Y. jaegeriana share a more similar climate with Y. brevifolia populations, and SDMs suggest that Y. jaegeriana species could extend further into the western Mojave absent limits on pollination, dispersal and other unaccounted factors. However, the Mojave Desert and surrounding regions are predicted to experience rapidly changing climate over the next several decades (Dai, 2013), and local adaptation within the range of each species could lead to differences in the population’s response to future changes in climate, particularly if climate change alters the seasonality of precipitation. Recruitment of Joshua tree populations follows narrow seasonal precipitation and temperature cues (Reynolds et al., 2012). Like other large and long-lived desert plants (Steenbergh and Lowe, 1969; Steenbergh and Lowe, 1977; Jordan and Nobel, 1979), seedling and juvenile Joshua trees are more vulnerable to drought-related mortality than adults (Esque et al., 2015). Rapid shifts in the climate regime could limit a population’s ability to recruit while adult Joshua trees continue to survive in undisturbed habitat, leading to lag effects (Svenning and Sandel, 2013) that may not be detected using the methods presented here. Lags in population change for long-lived plants are already observed among the large tree-like quiver tree (Aloidendron dichotomum) of South Africa and Namibia (Foden et al., 2007; Brodie et al., 2021).

Although the sister species of Joshua tree currently exist in allopatry, except for the narrow hybrid zone (Godsoe et al., 2009) – a pattern not uncommon to many sister species (Barton and Hewitt, 1985) – our SDM projections of potential habitat suggest this pattern may reflect ecological relationships other than habitat suitability, as previously hypothesized (Godsoe et al., 2009). Potential habitat for Y. brevifolia is closely associated with occupied habitat and does not extend far into the range of Y. jaegeriana east of the hybrid zone with a few scattered exceptions ( Figure 8 ). Conversely, for Y. jaegeriana west of the hybrid zone (where this species does not currently exist), extensive potential habitat is predicted and closely mimics the habitat of Y. brevifolia with few exceptions. Considering this pattern, as well as the tendency for introgression of Y. jaegeriana genetic material westward into Y. brevifolia populations (Starr et al., 2013; Royer et al., 2020), one wonders why Y. jaegeriana is not the dominant species in a front moving westward across the Mojave Desert. We suggest this may be largely dependent on biological factors that we did not account for in this analysis. Such factors include but are not limited to pollination interference by mis-matched pollinators, lack of vigor in hybrids for a variety of reasons, low capability of pollinator dispersal, reproductive isolation, and low seed dispersal distances (Vander Wall et al., 2006; Smith et al., 2008; Godsoe et al., 2009; Waitman et al., 2012; Royer et al., 2016; Royer et al., 2020). As Godsoe et al. (2009) previously suggested, the natural experiment between Joshua tree species in the hybrid zone is not sufficient to understand all the dynamics driving the patterns observed there. Further experimentation using known matrilines in combination with growth chambers and common gardens to tease out differences in species performance by genetics, physiology, and phenology would be a good start toward unraveling these patterns.