Local niche differences predict genotype associations in sister taxa of desert tortoise

To investigate spatial congruence between ecological niches and genotype in two allopatric species of desert tortoise that are species of conservation concern.


| INTRODUC TI ON
The relationships between the distributions of species and their ecological properties have long been central to biogeographic inquiry (Grinnell, 1917;MacArthur, 1972). More recently, quantitative methods to define ecological niches have become essential toolsets for investigating how species are distributed in environmental and geographic space (Araújo & Guisan, 2006;Guisan & Zimmermann, 2000). Species distribution modelling (SDM) is one such toolset that relates locations of species observations to explanatory variables hypothesized to influence or define a species' Hutchinsonian niche (Franklin, 2010). SDM quantifies the relationships between environmental conditions at locations where a species has been observed to those locations where it has not in order to predict how likely it is to occur at unobserved locations. These relationships, hereafter "species-environment relationships," are represented by model coefficients and can range from simple linear parameters in the most basic form of generalized linear modelling (GLM), to complex combinations of basis functions (e.g., polynomials, splines) in generalized additive modelling (GAM) or machine learning frameworks (Franklin, 2010).
For example, under conditions of allopatric speciation, SDM can be used to develop niche models for sister species in order to quantify niche differences using ecological overlap metrics (e.g., Broennimann et al., 2011;Godsoe, 2013;Rödder & Engler, 2011) and statistical tests (e.g., Nunes & Pearson, 2017;Warren et al., 2008;Warren, Glor, & Turelli, 2010). Substantial similarity between niches may suggest niche conservation between the two allopatric species (Warren et al., 2008;Wiens & Graham, 2005), such that the species-environment relationships are maintained through time and across taxa, even in the presence of environmental change or speciation (Wiens & Graham, 2005). However, the assumptions normally imposed in SDM may affect its utility for niche comparisons, including (a) imposing a single scale for all species-environment relationships and (b) fixing the spatial scale of those species-environment relationships across space. The former is well understood, with a consensus that no single scale is most appropriate for studying ecological phenomena (Levin, 1992;Wiens, 1989) and that understanding landscape structure and ecosystem processes may require multi-scale approaches (Rahbek & Graves, 2001;Seo, Thorne, Hannah, & Thuiller, 2009;Willis & Whittaker, 2002;Wu, 2004). This is evident in SDM where climate may dominate distributions at global to regional scales, yet topography and surface characteristics may influence species at regional to local scales (Beever, Swihart, & Bestelmeyer, 2006;Mackey & Lindenmayer, 2001).
The second constraint assumes that species-environment relationships (and the model coefficients that define them) do not vary across the geographic range of each species. Mixture zones, especially those with hybridization, may represent areas where a species' niche differs from the rest of its range due to local adaptation in different habitats (Lenormand, 2012). In these areas, species-environment relationships defined from the entirety of a species' range may not adequately characterize local differences, such that a single set of model coefficients may not represent the variation in species-environment relationships across space (Foody, 2004;Miller, 2012;Osborne, Foody, & Suárez-Seoane, 2007) because model coefficients may covary with location (Atkinson, 2001;Fotheringham, 1997). SDM practitioners have developed tools for incorporating potential spatial variation in species-environment relationships by calibrating separate models on subregions of a species' distribution (e.g., Peterson & Holt, 2003) or through the use of additive or multiplicative indicator variables. However, these solutions require a priori knowledge of the configuration of any hypothesized subregions and therefore limit exploration of spatial patterns in model coefficients.
Identifying spatially varying relationships across ecological boundaries (such as between two adjacent species) can highlight differences in data quality spanning those transition zones (e.g., differences in sampling regimes; Cheng & Fotheringham, 2013) and can provide information about the nature of the boundary, such as potential secondary contact between vicariant populations (Endler, 1977;Gay, Crochet, Bell, & Lenormand, 2008;Jiggins & Mallet, 2000). Through a conservation biogeography lens, boundaries and the gradients spanning them are key concepts that give context to phylogenetic diversity and the conditions that influence speciation (Hoffmann & Blows, 1994).
Tools such as geographically weighted regression (GWR; Fotheringham, Brunsdon, & Charlton, 2002) have been used to explore locally varying processes and their spatial scale across boundaries (Cheng & Fotheringham, 2013) and have also been used to investigate locally varying patterns in species richness (Bickford & Laffan, 2006;Foody, 2004;Holloway & Miller, 2015) and species-environment relationships (Kupfer & Farris, 2006;Miller, 2012;Miller, Franklin, & Aspinall, 2007;Osborne et al., 2007). In GWR, the spatial scale of species-environment relationships is represented with bandwidth parameters that determine the degree to which nearby observations are given higher weights than more distant ones (Fotheringham et al., 2002). Large bandwidths approximate classical regression with a single set of coefficients to represent species-environment relationships, while small bandwidths result in highly local models with species-environment relationships that can vary across the landscape. However, previous implementations of GWR have required a single bandwidth for all explanatory variables (Fotheringham et al., 2002), thus precluding a multi-scale approach.
A recent development has enabled estimation of separate spatial scales for each explanatory variable by optimizing multiple bandwidth parameters-multiscale geographically weighted regression (MGWR; Fotheringham, Yang, & Kang, 2017), which allows a multiscale approach to exploring species-environment relationships.
Here we develop a case study of two sister taxa to explore geographic patterns of niche differences between them, with focus on their differing conservation status. These two species, Gopherus agassizii (Agassiz's desert tortoise) and Gopherus morafkai (Morafka's desert tortoise) diverged approximately 6 Ma due to geographic isolation by the Bouse embayment, a putative marine transgression of the ancestral Gulf of California along the lower Colorado River, which has resulted in allopatric speciation (Murphy et al., 2011). These two cryptic species were only recently distinguished phylogenetically and taxonomically due to differences in genetics, reproductive ecology and seasonal activity (McLuckie, Lamb, Schwalbe, & McCord, 1999;Murphy et al., 2011), but are not readily distinguished morphologically. Prior to the taxonomic split, a distinct population segment (DPS) defined as the Mojave population (Figure 1; tortoises west and north of the Colorado River) was listed as threatened with extinction and given legal protection under the U.S. Endangered Species Act (ESA; Department of the Interior, 1990) and has received extensive monitoring yielding a wealth of georeferenced observations. The remaining "Sonoran Population," later elevated as the distinct species, G. morafkai, does not have the same legal protection or monitoring effort (Murphy et al., 2011;Service, 2015).
While the Colorado River defines the geographic division between the species, recent genetic work has identified a secondary contact zone where G. agassizii (the western species) occurs in a small population east of the Colorado River (Edwards et al., 2015;McLuckie et al., 1999). This secondary contact zone likely emerged only 2.5 ka as a result of avulsion in the Colorado River, but now G. agassizii in this zone are isolated from individuals occurring west of the Colorado River. This small population faces threats from increasing development in the region and is not legally protected.
This situation is further complicated by evidence of natural hybridization between G. agassizii and G. morafkai occurring in this secondary contact zone, and by the lack of a clear definition of habitat for this population of G. agassizii east of the river. Recent work has also suggested that this population occupies habitat with intermediate characteristics to that of the pure lineages (Edwards et al., 2015).
Habitat for G. agassizii and G. morafkai outside the contact zone is better defined, with known habitat ranging from valley bottoms and alkaline areas surrounding playas in the Mojave, to bajadas and alluvial fans, arroyos, rocky slopes and ridges in the upland regions in the Sonoran desert (Nussear & Tuberville, 2014). Differences in habitat characteristics between the two species span physiography (Nussear & Tuberville, 2014), geology (Burge, 1978), vegetation (Bury, Esque, DeFalco, & Medica, 1994) and climate (Nussear, Esque, Inman, Gass, & Thomas, 2009;Tracy et al., 2004). However, no formal tests of niche similarity have been conducted to date, further complicating delineations of these two iconic species.
We use SDM and MGWR in a coupled modelling approach to identify differences in the ecological niches of G. agassizii and G. morafkai, and explore spatially varying species-environment relationships in the recent secondary contact zone. We (a) formally test for differences in their ecological niches, (b) identify boundaries represented by differences in their niches and (c) determine which of three spatial F I G U R E 1 Study area used to create pooled-taxa species distribution models (light grey) and region of habitat for the two species of desert tortoise, Gopherus agassizii (Agassiz's tortoise; light orange) and Gopherus morafkai (Morafkai's tortoise; purple). The focal study area (thin black line) encompassing the contact zone was used to reduce computation time for local models and genotype assessment. The Colorado River (blue) separates California and Arizona and creates the division between the two species of desert tortoise. A Distinct Population Segment defined as the Mojave population includes individuals located west of the Colorado River. USA Contiguous Albers Equal Area Conic projection (SR-ORG:7301) delineations better describes landscape patterns of genotypic variation. These delineations include (a) the Colorado River (the current geographic boundary defining each species), (b) the Mojave and Sonoran Basin and Range ecotone, and (c) geographic patterns in local niche differences identified in this study. The results of this study will inform conservation planning across the transition zone of these two species.

| Study area
Our study area included the known range of G. agassizii and G. morafkai across 68,323 km 2 in the Southwestern United States, encompassing parts of California, Arizona, Nevada and Utah (Figure 1). This region is characterized as the Mojave Basin and Range Level III Ecoregion and Sonoran Basin and Range Level III Ecoregion (Wiken, Nava, & Griffith, 2011), hereafter the Mojave Desert and Sonoran Desert, respectively. The subregion encompassing the genetic sampling locations used by Edwards et al. (2015), hereafter referred to as the focal study area (Figure 1), offered an opportunity to explore spatial patterns in species-environment relationships across the ecotone between the Mojave and Sonoran deserts and in the secondary contact zone between G. agassizii and G. morafkai.

| Modelling overview
We developed a two-step modelling approach drawing on the strengths of both SDM and MGWR ( Figure 2) to explore spatial F I G U R E 2 Modelling Overview. Two-step modelling approach using species distribution modelling (SDM) and multiscale geographically weighted regression (MGWR) to explore spatial patterns in species-environment relationships of Gopherus agassizii and Gopherus morafkai.
In step 1, we use SDM to develop range-wide ecological niche models for each species separately and test hypotheses of niche equivalency. We also pool both species to develop a single model of their combined ecological niche and use the mapped residuals from this pooled model as a measure of local deviation. In step 2, we use MGWR to explore spatial patterns in the relationships between these residuals and explanatory variables that may enumerate differences between the two species and their hybrids within the focal study area TA B L E 1 Names, abbreviations and general description of 13 explanatory variables considered for modelling the pooled distribution of Gopherus agassizii (Agassiz's tortoise) and Gopherus morafkai (Morafkai's tortoise)

Citation
Gopherus agassizii (Agassiz's tortoise) patterns in species-environment relationships of G. agassizii and G. morafkai across this secondary contact zone. In the first step, we use SDM to develop range-wide ecological niche models for each species separately and test hypotheses that the niches of these two species are more different than would be expected by chance. We then pool location data for both species and develop a single model of their combined ecological niche and use the mapped residuals from this pooled model as a measure of local deviation. We assume that if the two species exhibit different ecological niches, residuals from a pooled model will represent how poorly their combined niche predicts the probability of presence at a given location. In the second step, we use MGWR to explore spatial patterns in the relationships between these residuals and hypothesized explanatory variables that may enumerate differences between the two species and their hybrids within the focal study area.

| Species distribution modelling
We observations for the two species available for model calibration. To reduce the effects of spatial sampling bias, we implement a background weight correction with the FactorBiasOut algorithm (Dudik, Phillips, & Schapire, 2005) and use a bias grid as an estimate of the sampling bias by creating a kernel density raster of observations for each species. The bandwidth for each kernel was estimated using cross-validation to minimize mean square error (Baddeley, Rubak, & Turner, 2015) and was linearly rescaled to 1-20 range to give greater background selection probability to areas with higher densities of observations (Elith et al., 2011).
In order to compare niches of G. agassizii and G. morafkai, a common set of explanatory variables is needed; we therefore considered a set of 13 explanatory variables (  (Elith et al., 2011). We stopped removing variables when a decrease of 0.05 in the area under the receiver operating characteristic curve (AUC; Fielding & Bell, 1997) was observed with 20% withheld test data.
The selected set of explanatory variables was used to calibrate models for each species separately using a bootstrap framework with 100 iterations and also to create a pooled model by treating the two species as a single taxon and pooling observations. We report the test AUC for each model, as well as relative contributions for each explanatory variable on training gain.

| Niche comparisons
We hypothesized that the niches of the two species would show similarities due to relatedness and niche conservation, but that differences would also be apparent due to geographic isolation over the past 6 million years. We therefore compared their niches using be used as a measure of variable importance (Phillips & Dudik, 2008).
We asked if the explanatory variables had the same importance to each species, such that high correlation of these contribution scores across species would indicate niche overlap and provide additional evidence that these two species share niche properties, whereas low correlation would suggest that some explanatory variables are more important to one species than the other. other explanatory variables at their mean value (Phillips & Dudik, 2008). Differences in the shape of species-environment response curves suggest niche differences between the two species. We hypothesized that differences due to prolonged geographic isolation over the past 6 million years would be limited to physiographic variables as a result of the substantial differences in terrain between occupied habitats. For example, G. agassizii occurs more often in valley bottoms and gentle slopes with smaller sediment size, while G. morafkai occurs more often in rockier slopes and bajadas (Nussear & Tuberville, 2014). We also assumed that differences in their responses to temperature would be minimal because these species have evolved in comparable climates.
Our third method to compare the species' niches used ran- For each of these randomization tests (equivalency and asymmetrical similarity), we define niche similarity with the expected fraction of shared presences overlap metric (ESP; Godsoe, 2013), which measures the degree to which two probability of presence surfaces agree.

| Local niche models and spatial scale
In order to further investigate differences in species-environment relationships between G. agassizii and G. morafkai, we calibrated local species-environment relationships within our focal study area around the secondary contact zone to estimate local variation that may exist across this region. We hypothesized that if differences in species-environment relationships were evident between the two species, then spatial gradients in those relationships may also be evident in the ecotone between them. We expected that these species-environment relationships would be expressed at different spatial scales due to differences between climatic and topographic constraints on distributions, where climate may affect distributions at regional scales, while topography may influence local scale patterns (Beever et al., 2006;Mackey & Lindenmayer, 2001 Tables S1.1-S1.4 in Appendix S1) in GRASS GIS with the "i.pca" module using normalization. We explore eigenvalues and their weight We therefore thinned the calibration data to 3 per 10 km 2 to reduce computation time, resulting in a dataset with 2,156 records. This thinning also created a uniform density of observations across the study area to minimize bias towards the more heavily sampled species, G. agassizii. MGWR can use adaptive bandwidths as estimates of spatial scale, and therefore, we use a Gaussian spatial kernel for each explanatory variable to allow each variable to converge on a separate bandwidth using AIC with small sample correction to avoid overfitting (Fotheringham et al., 2017). Non-linear regression coefficients were not considered as they have not been implemented in MGWR, and because at very local scales, species-environment relationships are expected to approximate linear responses due to the limits in the local range of each explanatory variable (Fotheringham et al., 2002). Bandwidths and their approximate spatial scale are reported for each explanatory variable along with local parameter estimates and model R 2 . We report the spatial scale of each explanatory variable as the product of the average distance between locations in our calibration dataset and the bandwidth, and use inverse distance weighted interpolation to create regression coefficient maps from local parameter estimates for areas that were thinned prior to running MGWR.

| Habitat-genotype association
We hypothesized that landscape patterns in the interpolated MGWR coefficient maps (representing spatially varying species-environment relationships) would be congruent with previously reported phylogenetic differences found in the secondary contact zone identified by Edwards et al. (2015). We represented the phylogenetic structure of sampled populations using the admixture proportion (Q) of a pure G. agassizii genotype from STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). This index, hereafter genotype association index, represents the probability that an individual contained G. agassizii genotypes (Edwards et al., 2015) and was interpolated across our study area using inverse distance weighting to create a map for the Mojave genotype. We used Kendall's rank correlation coefficient (tau) for paired samples to assess correlations between each explanatory variable's local regression coefficients and the genotype association index. A nonparametric test was chosen because our genotype association index did not meet assumptions of normality.
We then asked if natural divisions in the local species-environment relationships exist, and if present, do they coincide spatially with the genotype association index. We identified divisions with K-medoids optimal partitioning in multivariate space of the local regression coefficient maps for all explanatory variables using the package "cluster" (Maechler, Rousseeuw, Struyf, Hubert, & Hornik, 2016) in R (R Core Team, 2016). The optimal number of clusters was estimated by minimizing within-cluster variance (Hennig & Liao, 2013). Cluster assignments were mapped back to geographic space and compared to the genotype association index. Here we used spatial autoregressive lag models (SAR lag; Anselin, 2001) with the package "spdep" (Bivand & Piras, 2015) in R (R Core Team, 2016) to determine which of three delineations best explained the gen-

| Species distribution modelling
The single set of variables selected to describe the ecological niche of each species included the following: precipitation of the driest month (Ppt_dry), precipitation seasonality (Ppt_CV), surface texture (Surf_Text), soil moisture (S_moist), temperature evenness (Isotherm), photosynthetic activity (Veg_Amp), topographic position index (Topo_ Index) and surface material (Surf_Mat; see Table 1 for descriptions), and resulted in models for G. agassizii and G. morafkai with test AUC scores of 0.733 and 0.875, respectively ( Figure 3). The pooled model showed reduced performance, with a test AUC score of 0.697, suggesting that the ecological niches of each species were different from one another and not well represented with a single, pooled model.

| Local niche models and spatial scale
The reduced set of nine principal components identified to investigate spatial patterns in local species-environment relationships in-  Local regression coefficient maps are provided in Figure S1.1.

| Habitat-genotype association
Kendall's tau values representing the degree to which local regression coefficient maps from MGWR were correlated with our genotype association index, ranged from −0.43 to 0.40 (Table 4) and indicated a modest overall agreement between any given species-environment relationship and genotype. However, when considered together, we identified two multivariate clusters in these local regression coefficient maps, which, when mapped back to geographic space, were largely coincident with the boundary separating G. agassizii and G. morafkai ( Figure 5). However, this division between the two clusters did not exactly coincide with the Colorado River, instead suggesting a boundary approximately 40 km to the east of the Colorado River in the northern portion of the focal study area ( Figure 5). The mean genotype association index for the two clusters was 0.98 and 0.15 for the Mojave and Sonoran clusters, respectively, indicating that the Mojave cluster was most strongly associated with the Mojave genotype and the Sonoran cluster was not. The mapped clusters of local species-environment relationships were better able to predict the genotype association index than either the ecoregions or the geographic delineation between the species, with a ∆AIC score of >2 between the next best SAR lag models (Table 5). Overlap between the mapped clusters and the Mojave and Sonoran ecoregions suggested that the Mojave cluster was more closely aligned with the Mojave ecoregion than with the current geographic delineation of the Mojave population of G. agassizii (Table 6). In contrast, the Sonoran cluster most closely aligned with the current geographic delineation of G. morafkai, indicating that the current geographic delineation of G. morafkai is a better representation of Sonoran habitat in the focal study area than the Sonoran ecoregion alone (Table 6).

| D ISCUSS I ON
We explore ecological niche differences between two allopatric species of conservation concern, G. agassizii and G. morafkai, and find that while these two species occupy broadly similar ecological niches in their respective ecoregions, they differ subtly in their selection of habitat. Moreover, spatial scale differed among key species-environment relationships, confirming expectations that climate may dominate species' distributions at coarse scales, while responses to topography and land surface characteristics may be more apparent at fine scales (Pearson & Dawson, 2003). Habitat barriers such as water and developed surfaces (i.e., lakes, road and cities) had negative effects on habitat at coarse scales independent of location, while soil conditions, vegetation and physiographic characteristics exhibited local effects that varied across the region encompassing the recent secondary contact zone. We also find that local variation in species-environment relationships provided greater support for the phylogenetic differences observed among individuals than does the current geographic delineation between the two species. Our results contribute additional evidence that G. agassizii and hybridized individuals east of the Colorado River and west of Kingman AZ ( Figure 5) exhibit ecological niches that are more similar to G. agassizii in the rest of its range than to their proximal sister taxa, G. morafkai.

| In support of phylogenetic boundaries: Local species-environment relationships
We found evidence for two, but not three, multivariate clusters in the local regression coefficient maps. A third category, if coinciding with regions containing hybrids, might suggest that hybridized individuals select habitat in ways that are locally different than either of the two pure genotypes. Previous work has shown that these hybrids do occupy habitats with characteristics that span those of both their Mojave and Sonoran parental lineages in terms of topographic, surface textural and vegetation characteristics (Edwards et al., 2015). Our delineation of two categories does not counter these findings because here we explore local niche differences, that is local differences in species-environment relationships, rather than differences in occupied habitat as was explored by Edwards et al. (2015). In this simple example, these two groups show similar positive local species-environment relationships even though they occupy different regions of this hypothetical gradient (i.e., values of 10 vs. 5).
Using a coupled approach with species distribution modelling and multiscale geographically weighted regression, we find that individuals in the secondary contact zone exhibit habitat preferences that are more akin to G. agassizii than G. morafkai even though some of them occupy habitats that are only marginally different from either parental lineage. This is where coupling SDM with local modelling methods such as MGWR departs from traditional habitat assessments-differences in local habitat selection are uncovered rather than differences in local occupied habitat. Local habitat selection is evaluated in context of nearby environments, while occupied habitat is a measure of differences between occupied areas.
In the case of G. agassizii and G. morafkai, regional differences in occupied habitat are clearly evident. Differences span climate, vegetation, physiography and geology (Nussear & Tuberville, 2014) and are consistent with the subtle differences we found in their ecological niches when quantified with SDM based on range-wide species-environment relationships. For example, when we modelled distribution of each species separately, we found differences in species-environment relationships for seasonal vegetation amplitude (Veg_Amp; Figure S1.2), an explanatory variable describing the seasonal vegetation green-up potential (Meier & Brown, 2014). We also found subtle differences in the range-wide relationships for surface texture (Surf_Text; Figure S1.2) and topographic position index describing physiographic relief (Topo_Index; Figure S1.2). G. agassizii tend to occupy regions with finer surface textures such as alluvial soils, while G. morafkai occur more often in rocky soils and bajadas (Van Devender, 2006). In contrast, both species share similar rangewide unimodal species-environment relationships for summer soil moisture content (S_moist; Figure S1.2), thereby avoiding very dry and very wet soils. Similarly, both G. agassizii and G. morafkai appear to have range limits defined by cold winter temperatures, as each can tolerate high summer temperatures through behavioural aestivation (Nussear & Tuberville, 2014). This suggests that while the two species occupy different habitats, they exhibit similar selection for certain environmental conditions.
In contrast, we found substantial differences in species-environment relationships at local scales, where soils (SOIL2, SOIL3), precipitation (CLIM1) and vegetation (VEG1, VEG3) variables were optimized with short bandwidths. This suggests that local variations in the species-environment relationships of these explanatory variables contribute to overall niche differences between the two species. For example, the SOIL2 PCA showed a spatial scale of ~50 km and was most influenced by summer and winter potential evapotranspiration. The local regression coefficient map for this variable tended to show positive species-environment relationships west of the Colorado River, and negative relationships east of the River, suggesting that individuals west of the Colorado River tended to select habitat with higher potential evapotranspiration given locally available conditions ( Figure S1.1). Differences in these local species-environment relationships, such as precipitation (e.g., summer and winter; CLIM1), terrain (e.g., slope and rockiness; PHYS1) and vegetation (e.g., phenology and canopy growth; Here, we find that the most parsimonious spatial model explaining the landscape pattern of genotype association was the two multivariate clusters of local species-environment relationships rather than the Mojave and Sonoran ecoregions or the current geographic delineation of the two species and their protection status.

| Importance for conservation
Efforts to preserve biodiversity have placed new emphasis on measures of biodiversity beyond taxonomic diversity.
Understanding landscape patterns in phylogenetic diversity is especially important to conservation goals aimed at maximizing the resilience of biodiversity in the face of rapid global change (Flynn et al., 2011;Legendre et al., 2005), and for identifying conditions where recent lineage divergence has contributed to local niche differences that may aid in adapting to changing environments  Moritz, 2002). Identifying spatially structured variation in habitat selection, coupled with an understanding of genotypic structure, is therefore important for predicting potential outcomes of spatial conservation decisions (Ferrier & Drielsma, 2010;Whittaker et al., 2005). Often, conservation prioritization focuses on hotspots (Myers et al., 2000;Naeem et al., 2012;Winter et al., 2013) delineated on the basis of taxonomic diversity (Ferrier et al., 2004;Myers et al., 2000), phylogenetic diversity (Crozier, 1997;Helmus et al., 2007;Scoble & Lowe, 2010;Vandergast et al., 2013;Wood et al., 2013) or measures of evolutionary potential, such as sequence diversity (Tamura & Nei, 1993) or divergence (Nei & Li, 1979). However, the ability to compare landscape measures of genetic diversity to measures of local niche differences and habitat selection (e.g., species-environment relationships) presents new opportunities to investigate the confluence of genetics and ecology in context of conservation biogeography.
Conservation managers tasked as stewards of healthy and sustainable ecosystems often request spatially explicit information that supports management objectives. Lake Mead National Recreation Area, a unit within the National Park Service, is the unit responsible for stewardship of ~6 million ha of land in southern Nevada and northwest Arizona. Park managers seek information on tortoise distributions and lineages in order to prioritize protection and restoration of tortoise habitat impacted by invasive weeds, fire, road disturbance, recreation and development (Brooks & Esque, 2002;Esque et al., 2010;Lovich et al., 2011).

| A novel coupled modelling approach
The use of local regression to explore spatial variation in species-environment relationships is not new to SDM but has been difficult to apply given the widespread reliance on binary (presence-absence or presence-background) calibration data necessitating logistic regression. Local models using a logistic regression framework can suffer from complete separation of response classes at fine spatial scales (Fotheringham et al., 2002), thereby forcing models to use large bandwidths approximating range-wide models (Miller, 2012). This is especially problematic when calibration datasets exhibit extreme sampling bias. We mitigate against this problem by calibrating a local MGWR model on the residuals of a pooled niche model from both taxa to explore local deviation in species-environment relationships.
Modelling residuals enable the use of local Gaussian models, and multiscale local regression methods allow more flexible application of the method to ecological data where relationships between response and explanatory variables may be apparent at varying scales.
Cases of allopatric speciation are especially well suited to these coupled methods because gradients in ecological niche conservation across two taxa can be explored at a local level within hypothesized mixture zones. Moreover, delimiting small regions of interest is necessary when computationally intensive MGWR models require extreme processing times due to their use of iterative back-fitting algorithms to fit optimal bandwidth vectors (Fotheringham et al., 2017). In contrast, presence-background SDM methods assume that the entirety of a species' range is sampled (Elith & Leathwick, 2009;Franklin, 2010) and require large study areas. This coupling approach allows each method to use an appropriate spatial domain for its respective model assumptions and computation limitations. In this way, MGWR can be used in a subregion of the species' distribution to explore local variation in species-environment relationships as expressed in deviations from these predictions, that is residuals.

ACK N OWLED G EM ENTS
We thank J. Meuller and M. Sappington for enumerating Lake Mead National Recreation Area research needs, and A. Vandergast and K. Drake for providing helpful comments and suggestions. We are also grateful for 3 anonymous referees for their careful reading and thoughtful contributions to this work. Any use of trade, product or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.

DATA AVA I L A B I L I T Y
All relevant data are hosted by ScienceBase, a public repository hosted by the U.S. Geological Survey. Data can be accessed at the following DOI: https ://doi.org/10.5066/P91V2S8C.