Evaluating current and future range limits of an endangered, keystone rodent (Dipodomys ingens)

We use current models of species distribution to predict the future habitat suitability for an endangered, keystone rodent (Dipodomys ingens, giant kangaroo rat). We incorporate the possibility of local adaptation and novel competitive interactions to improve the predictive accuracy of our Maxent models.


| INTRODUC TI ON
Climate is often considered the single most important factor limiting species' ranges (Merriam, 1894;Peterson et al., 2011). Contemporary approaches to modelling species distributions and predicting range shifts under climate change suggest species will relocate to track their climatic niches (Parmesan, 2006). This response is governed by niche conservatism, which occurs when species retain ecological traits related to their niche over time (Wiens et al., 2010). Here we define "niche" as the combination of abiotic and biotic factors that a species experiences and those factors' distribution in geographic space (Hutchinson, 1957). The resulting potential range is the distribution of the bioclimatic envelope, and the realized range is the product of the constraints imposed by biotic interactions (Guisan & Thuiller, 2005). The interaction between abiotic and biotic factors is important in determining whether species follow or deviate from niche conservatism. While coarse-scale studies have detected movements poleward or towards higher elevations, suggesting that species are tracking their climatic niches (Parmesan & Yohe, 2003), long-term and fine-scale studies have revealed more idiosyncratic responses (Gibson-Reinemer & Rahel, 2015;Moritz et al., 2008).
Three potential problems with expecting species to exclusively track their contemporary, realized, abiotic niche are as follows: (a) most studies focus on temperature without regard to precipitation (Rapacciulo et al., 2014); (b) populations may have adapted to local climate regimes (Gibson-Reinemer & Rahel, 2015); and (c) species may be constrained by biotic factors that are only currently correlated with abiotic proxies (Guisan & Thuiller, 2005).
Even with the improved accuracy from multiple climatic variables, broad correlations between species distributions and climate may obscure local adaptations. Some species' ranges may shift to track a fluctuating niche, while other species' ranges remain in place, exposing them to changes in climate (Holt, 1990). Non-contiguous populations within a species could experience dissimilar climatic conditions, and over time, the population-level niche could vary based on local conditions (Gibson-Reinemer & Rahel, 2015). Modelling genetically unique subpopulations should improve the predictive accuracy of distribution models (Gonzalez, 2011).
Considering multiple aspects of climate may make predictions more accurate but neglects the potentially powerful influence of biotic interactions (Post, 2013). Small mammal communities are structured by internal mechanisms to avoid competition (e.g. microhabitat partitioning, resource selection, temporal separation) (Price, 1978). While these processes may reduce the effects of competition within communities, the boundaries between communities could serve as intense areas of competition, which can influence species distributions (Boulangeat, Gravel, & Thuiller, 2012;Gaston, 2003). Climate change is expected to affect community assembly, with species expanding across community boundaries at different rates, inducing novel competitive interactions (Montoya & Raffaelli, 2010). Most species distribution models-whether explicitly stated or not-incorporate biotic interactions through indirect abiotic proxies, such as attempting to capture competitive species boundaries with differences in temperature (Guisan & Thuiller, 2005).
However, because species respond to climate change at various temporal scales and via different mechanisms, the relationship between climate and competition may be decoupled when projecting into the future (Elith & Leathwick, 2009;Guisan & Thuiller, 2005;Lurgi, Lopez, & Montoya, 2012). Consideration of biotic limitations directly may better explain the mechanisms limiting current distributions and, therefore, more accurately forecast species distributions under future climates.
The giant kangaroo rat (Dipodomys ingens) is an ideal species to test the inclusion of local adaptation and biotic interactions in species distribution modelling. D. ingens are a state and federally listed endangered species endemic to California's Great Central Valley, a global biodiversity hotspot (Myers, Mittermeier, Mittermeier, Fonseca, & Kent, 2000;USFWS, 1987). They currently only reside in a few remnant and isolated populations, the largest of which are found in the Carrizo Plain National Monument ("Carrizo") and the Ciervo-Panoche Natural Area ("Panoche"). These two populations are geographically isolated by ~150 km, and recent work suggests deep genetic isolation dating back at least 6,000 years, making local adaptation more likely (Statham, Bean, Alexander, Westphal, & Sacks, 2019). Although these populations may experience different local conditions, the species exists within a narrow range of habitat characteristics .
These habitat specialists thrive in desert-grasslands with low annual precipitation (<30 cm), sandy loam soils, and flat or very low-grade slopes (Grinnell, 1932;Hawbecker, 1951). Their burrow structures are used by other vertebrate species, including some classified as endangered, and their seed caching may serve as a food source for invertebrate inhabitants (Prugh & Brashares, 2012). They are thus considered ecosystem engineers and a keystone species (Prugh & Brashares, 2012). Understanding the limitations to D. ingens range will aid in crafting more effective conservation strategies and help protect an endangered ecosystem.
Precipitation appears to play a key role in limiting D. ingens distributions . Persistence in areas of low annual precipitation is limited by food resources, particularly after consecutive years of low rainfall. However, the specific mechanisms by which precipitation limits D. ingens in the wetter parts of their range are undetermined (Bean, 2012). D. ingens could be limited by precipitation in several ways, leading to three non-exclusive mechanistic hypotheses. First, as a burrowing, seed-caching species, high precipitation could directly affect their ability to maintain burrows, or cause seed spoilage, depleting seasonal food stores ("Precipitation" hypothesis, Valone, Brown, & Jacobi, 1995). Second, higher precipitation causes increased growth of dense, non-native vegetation (e.g. Bromus madritensis spp. rubens) that impedes D. ingens movement, decreasing foraging efficiency and increasing the risk of predation ("Vegetation" hypothesis, Germano, Rathbun, & Saslaw, 2012). Climatic water deficit combines the effects of precipitation, temperature and radiation to provide a metric of available water (Flint & Flint, 2014). A lower deficit would result in increased soil moisture, which could improve growth conditions for dense nonnative vegetation, inhibiting saltatorial movements (Stephenson, 1998). Finally, while D. ingens is considered competitively dominant within its community, wetter areas of their range are suitable to potential novel competitors ("Competition" hypothesis, Grinnell, 1932;Prugh & Brashares, 2012;Bean, 2012). Populations in areas of higher precipitation have a greater chance of interacting with California ground squirrels (Otospermophilus beecheyi) which currently occupy the fluctuating wetter edges of D. ingens habitat. Direct contact with this aggressive species could result in antagonistic interactions (Trulio, 1996). In fact, Williams, Germano, and Tordoff (1993) observed O. beecheyi raiding D. ingens seed caches and report anecdotal evidence of a direct mortality or eviction event of D. ingens due to burrow invasion. O. beecheyi is a diet and habitat generalist that has shown little response to climate change over the last century and could pose a threat to D. ingens expansion into wetter territory (Eastman, Morelli, Rowe, Conroy, & Moritz, 2012;Grinnell & Dixon, 1918;Hubbart, 2012;Moritz et al., 2008). We investigate this species particularly because of its life history and current apparent distribution that create its potential for competition. As the largest of the kangaroo rats, D. ingens dominates its current community, thus we suggest that it would require a larger, more aggressive species to interrupt or displace this hierarchy. While D. ingens may encounter a suite of novel interactions with other species, O. beecheyi, according to current conditions, appears to be the most likely.
The Carrizo and Panoche populations experience different climatic regimes-the Panoche receives 100 mm (~45%) greater average annual precipitation than the Carrizo. It is therefore possible that the two populations of D. ingens adapted to their respective local conditions. This could reduce niche overlap between populations, but predict less range contraction, given the populations are adapted to a wider range of climatic conditions.
To test three hypotheses of D. ingens range limitation, we identified important variables in species distribution models. By incorporating more direct mechanisms-that is local adaptation and biotic interactions-we attempted to improve the predictive accuracy of habitat suitability from species distribution models for D. ingens. We created a suite of distribution models using Maxent including: (1) rangewide D. ingens, (2) population-specific D. ingens in the Panoche, (3) population-specific D. ingens in the Carrizo and the same set of models for O. beecheyi (4-6). To test the efficacy of the models including local adaptation, we compared the rangewide model (1) to the population-level niche models (2 and 3) to determine differences of population-level models. Then, we used the Panoche and Carrizo models to project the population-level estimates of habitat suitability rangewide, for both current and future climate. Finally, we estimated niche overlap between D. ingens and O. beecheyi to assess the possibility of a competitive interaction limiting D. ingens range expansion. To understand how different approaches to modelling contemporary range limits might impact future climate projections and to inform future management, we projected the top performing models into the future.

| Study area
The San Joaquin Valley, the southern portion of California's Great Central Valley, is a desert characterized by mild winters with low rainfall, and hot, dry summers (Germano et al., 2011). Much of the land includes agriculture and oil development that has replaced D. ingens habitat and is also at risk for solar energy development (USFWS, 2010;Williams et al., 1998). In California, the average annual temperature has increased 1.7°C over the last century, 70% higher than the U.S. national average (Moser, Franco, Pittiglio, Chou, & Cayan, 2009). During this time, the state also experienced extended periods of drought (National Drought Mitigation Center, 2016). The maximum temperature in the Panoche is expected to increase between 2.5 and 4.1°C by the years 2070-2099. Annual precipitation is expected to remain stable or to increase by up to 28 mm; the temperature in the Carrizo is predicted to increase by 1.8 to 3.4°C by the years 2070-2099, and precipitation is expected to decrease between 30 and 38 mm (Cal-Adapt, 2016).

| Live trapping
Occurrence data came from trapping across the range in 2010-2017.
Target-based trapping consisted of using extra-long Sherman live traps to identify the presence of D. ingens at sites with occupancy signs, such as burrows mounds or scat. Trapping protocol followed Prugh and Brashares (2010) and Alexander (2016). For modelling purposes, presence was recorded at the site level to reduce autocorrelation and to be more appropriate for the scale of the environmental variables. All work was performed under Humboldt State Animal Care Protocols 13-14.W.109-A and 16/17.W.96-A and followed American Society of Mammalogists guidelines (Sikes et al., 2011).
Our trapping locations included 197 presences in the Panoche, 101 in the Carrizo, and 334 rangewide, covering much of the range of D. ingens. However, we were unable to acquire observed presences in the Lokern and Elk Hills regions, areas of known GKR occupancy.
We therefore supplemented our data with 33 occurrences from the California Natural Diversity Database (CDFW, 2018).

| Environmental variables
We acquired climatic data from the Basic Characterization Model (Flint & Flint, 2014). We used a digital elevation model to calculate slope (USGS, 2013) and used estimates of soil texture to incorporate important aspects of burrowing requirements (NRCS, 2003). All environmental variables were resampled to the coarsest resolution for modelling, about 900 m (Table 1). While this scale is a product of data availability, it is also appropriate for modelling this species because individuals have a home range of about 200-1,000 m 2 (Cooper & Randall, 2007). This means that one cell could represent the environmental variables selected by an individual, or a group of individuals, which we classify as a site.

| Driving surveys
In 2017, we conducted driving surveys to record sightings of O. beecheyi. Due to logistics of travelling between the two populations, we were restricted to strategically surveying only one population extent. We considered the Panoche population of D. ingens to be at the greatest risk of competition with O. beecheyi, and thus, we drove our surveys within a 50 km radius of the Panoche (Figure 1).
With an additional observer, we drove five survey routes ranging between 130 and 190 km. In order to survey a range of environmental conditions, the five routes were stratified into five equally  (Flint & Flint, 2014). Surveys occurred every 4 km along the survey route. Two observers used binoculars to search the 360° visual field surrounding the survey point for O. beecheyi for 2 min (Downey, 2003). We estimated that, due to topography and land use, even the greatest distance from the survey point to a squirrel detection was within a 900 m pixel. If a squirrel was detected before the 2 min were complete, the point was coded as a presence and the survey continued to the next location. Opportunistic ground squirrel sightings during the surveys and during D. ingens trapping were also recorded.

| Species distribution modelling
We used Maxent to create species distribution models for D. ingens and O. beecheyi (Phillips, Anderson, & Schapire, 2006). Maxent is a machine-learning algorithm that calculates the maximum entropy probability distribution of occurrence points under a set of environmental constraints. It assumes that the distribution of a given species is driven by some combination of this set of environmental variables and can explore complex relationships between them.
Maxent is useful for studying rare and declining species because it requires the input of only presence locations and performs well with small data sets. Maxent thins presence locations by interpreting a single point per pixel, according to the resolution of the input variables, reducing the risk of autocorrelation. Maxent samples background locations (default = 10,000) to compare to these presence locations. No inference is made as to whether background locations are presence or absence sites, hence they are considered "available." Once a model is estimated, Maxent projects it to the entire study area, resulting in a map of values ranging from zero to approximately one. These values have multiple interpretations, but not all are appropriate for every species. For the purposes of this study, Maxent values are referred to as habitat suitability. Bean, Prugh, et al. (2014) found that Maxent accurately predicts long-term habitat suitability for D. ingens at a coarse temporal scale, which is appropriate for our species and environmental data.
We first created a suite of models of current D. ingens distribution and then included future projections for temperature and precipitation, according to the Community Climate System Model's emissions scenario under the most severe representative concentration pathway (rcp8.5) (Gent et al., 2011). We chose to model the future predictions under this pathway to capture the greatest potential change.
Lower emissions scenarios could underestimate the potential impact of climate change and rcp8.5 allows us to compare the current model to the riskiest scenario, which may be more likely to occur depending on current and future policy changes and political climates.
For current climatic conditions, we created species distribution models for three subsets of the D. ingens range: Panoche, Carrizo and a rangewide model that contained both populations. After selecting biologically relevant predictor variables, we decoupled predictor variables with |r| greater than 0.7, meaning the variables remained as possible environmental predictors, but were not paired in a single model (Pearson, 1920). Mean annual temperature was highly correlated with all but the soil characteristics and was thus removed as a predicator. Using a jackknife test of variable importance and each variable's per cent contribution to the models, we finalized the set of variables available for candidate models (Table 1). Aspect was removed from the final candidate model set for lack of contribution in preliminary model evaluation. We adjusted the beta value, a regularization parameter designed to optimize data fitting, on the preliminary models (β = 0.5, 1 and 2) and used Akaike's information criterion corrected for small sample size (AICc) to assess model fit.
The top models were all created using β = 1, and thus, we determined that this value was the most appropriate (Morales, Fernández, & Baca-González, 2017). We used the default auto-features setting to select feature types and selected the complementary log-log TA B L E 1 Environmental variable layers that were included in the candidate model sets for species distribution models for Dipodomys ingens and Otospermophilus beecheyi  Information Table S1.1 in Appendix S1).
The best model for each geographic subset was chosen using AICc, and model discrimination ability was assessed using area under the curve (AUC) and the Continuous Boyce Index (CBI; Hirzel, Lay, Helfer, Randin, & Guisan, 2006). CBI is an improvement of the Boyce Index that uses a moving window rather than discrete bins to estimate the Spearman rank correlation of the ratio of predicted to expected presence locations (Boyce, Vernier, Nielson, & Schmieglow, 2002). We calculated CBI using ten bins with a moving window between the maximum and minimum suitability values. We used 75% of the data to train the model and 25% to test using CBI. The result was a value between negative one and one, where a positive one represents a model that perfectly predicted presence, zero means the model is no better than random chance, and negative values predicted presence in the areas of the lowest suitability (Hirzel et al., 2006). We used only GBIF locations to create the California-wide models, to avoid over-sampling our study areas. For the Panoche and Carrizo models, we included survey locations and anecdotal sightings, but thinned the locations using the geogThin function in R package "enmSdm" so that clusters of squirrels within one pixel were represented by a single location (Muscarella et al., 2014;Smith, 2018).
To visualize contemporary and future distributions and inform management, we selected threshold values to convert continuous habitat suitability models to maps excluding the lowest suitability values. For all extents, the threshold for "not suitable" was selected based on the 5th percentile suitability value of all occurrence points within the respective current model. We subdivided suitability into two categories, low and high. We considered "low" suitability to be between the 5th and 50th percentile, and all greater values were considered "high" suitability. We used the same values for thresholding future models. These threshold values were only used for mapping purposes; un-thresholded models were used for further analyses.

| Niche overlap
The top D. ingens population models were compared to evaluate niche overlap as evidence for separate adaptation. Subsequently, overlap between the same models was calculated based on future climate projections. Low niche overlap would indicate that populations experience different environmental regimes and could suggest adaptation. We used Warrens I (Warren, Glor, & Turelli, 2008) to calculate the degree of niche overlap, a commonly used strategy for evaluating niche conservatism using habitat suitability models (Wiens et al., 2010). This similarity statistic assumes that habitat suitability is an appropriate measure of the niche, which is reasonable for our coarse temporal scale with distinct differences in each population's environmental variables (see Supporting Information

| Dipodomys ingens
The top models for two of three geographic extents were unambiguous; all other models had ΔAICc greater than seven (see Supporting Information Table S1.1 in Appendix S1). For the Carrizo model, there were two models within seven AICc, but we chose the model with the lowest AICc value. The Panoche model included climatic water deficit, minimum temperature, maximum temperature, per cent clay and slope (Table 2, Figure 2). The Carrizo model included precipitation, maximum temperature, per cent clay and slope (Table 2, Figure 2). The rangewide model included precipitation, minimum temperature, maximum temperature, clay and slope (Table 2, Figure 3). The response curves indicated selection for low slope, high minimum temperature, low precipitation and a peak in suitability around 23% clay in the soil (see Supporting Information Figure S3.3 in Appendix S3). All top models had AUC scores above 0.9 (see Supporting Information Figure S4.4 in Appendix S4   and future (2070-2099) species distribution models for giant kangaroo rats (Dipodomys ingens). Models were trained locally and projected rangewide, San Joaquin Valley, California. Higher suitability is indicated by dark shading. The model trained in the Panoche is shown in blue, and the model trained in the Carrizo is shown in red. Purple areas indicate overlap between the two models. The grey shading indicates the historical distribution presented by Williams et al. (1992), and the crosshatching shows areas of agriculture or residential use. The future map (right) was projected using CCSM4 rcp8.5 values were only used for map construction, making viewing of three separate models with the same extent possible in a static image. However, all analyses were performed using the full models.
When populations were modelled together in the rangewide model, the prediction of future habitat suitability was low and covered an extremely small area (Figure 3). However, when the two populations were modelled separately and each projected into the future (2070-2099), less range contraction was predicted ( Figure 2). The future projection of the historical model also predicted far less range contraction than the rangewide model ( Figure 4).

| Otospermophilus beecheyi
The top models for all three extents-Panoche, Carrizo and California-wide-differed but were unambiguous (Table 2,

| Niche overlap
All

| D ISCUSS I ON
Our models of local, contemporary D. ingens habitat suitability more accurately reflect the true distribution represented in the San Joaquin Valley. We believe this translates into a useful prediction of future habitat suitability. When populations were modelled together, a single, rangewide model predicted extreme future range contraction, with almost complete extirpation from their current range.
However, population-specific modelling provided a more optimistic forecast. The future prediction of the historical model, which represented wider niche breadth, was the most optimistic prediction.
Used in combination, these models highlight areas of consistently suitable habitat, taking into account historical niche, as well as local adaptation. The areas of overlap likely represent patches appropriate for conservation prioritization.
Species distribution models of D. ingens in the past have focused on rangewide or Carrizo-specific models, leaving the Panoche population relatively underrepresented .
Habitat For D. ingens, the environmental factors included in the top local models were similar for both populations and consistent with previous literature Grinnell, 1932;Williams, 1992), with the exception of the precipitation variable. The limitation imposed by precipitation could be defined by either excessive or insufficient rainfall. In the Carrizo model, precipitation was an influential variable, contributing 45% to the model (Table 2). Excessive precipitation may cause direct effects of water infiltration and burrow collapse (Germano, Rathbun, & Saslaw, 2001). Alternatively, precipitation could better define the areas limited by too little water, hindering primary productivity. The Panoche model included climatic water deficit rather than mean annual precipitation. This is consistent with the Vegetation Hypothesis, suggesting D. ingens in the Panoche are more limited by consequences of dense vegetation due to increased water availability than they are by direct effects of precipitation. Tingley et al. (2009) compared niche shifts of 53 bird species from historical (1911)(1912)(1913)(1914)(1915)(1916)(1917)(1918)(1919)(1920)(1921)(1922)(1923)(1924)(1925)(1926)(1927)(1928)(1929) to contemporary (2003)(2004)(2005)(2006)(2007)(2008) ranges, showing that species' responses to changing climate could not be predicted solely from temperature or precipitation, but estimates were far more accurate when the two were used in tandem.
The local models predicted far less range contraction than the rangewide model. However, the historical model, which modelled the species' historical environmental niche, predicted even less range contraction than the local models. This could be because it  (Good, Williams, Ralls, & Fleischer, 1997;Loew, Williams, Ralls, Pilgrim, & Fleischer, 2005) found relatively low genetic differentiation in neutral microsatellites between the Panoche and Carrizo,   indicating recent or no divergence, more recent work relying on unique haplotypes in cytochrome b suggested a much earlier (~6,800 years) splitting time (Statham et al., 2019). In two divergent populations with large effective population size, selection for specific climatic regimes may occur relatively rapidly (Hoffman & Sgrò, 2011). Additional work is needed to determine whether these populations-or others throughout the range-have in fact adapted to local climatic conditions or maintained niche conservatism.
The potential for competition was approximated using niche overlap between species. The results presented were estimated from the California-wide niche of O. beecheyi, which captures much of the range of the environmental associations of this generalist species ( Figure 5). This comparison suggested that niche overlap was greatest in the Panoche population, and the local models highlight the potential for interaction on a finer scale. Here we focus on the Panoche populations due to a greater threat of competition, as well as the more robust O. beecheyi model, as a result of survey efforts.
The areas of niche overlap between the local species models in the Panoche, the population of greatest concern for competition, are within confirmed current D. ingens habitat, including areas where the two species have been observed in close proximity ( Figure 5) The current recovery plan for D. ingens highlights the importance of acquiring and conserving specific locations of confirmed occupancy (Williams et al., 1998). An updated management plan for D. ingens should consider local adaptation, biotic interactions and historical, current and future habitat suitability. Suitable areas of overlap between the historical, local and rangewide models should be of particular concern and importance, but of these three potential futures, the overlap between historical and local models should be a priority (Figure 7). This combination retains historical niche characteristics while incorporating the effects of local adaptation over time. Within the priority areas, those occupied by O. beecheyi populations should be considered less than ideal habitat.
F I G U R E 7 Four separate future projections of giant kangaroo rat (Dipodomys ingens) habitat suitability combined. Two models were trained locally, within the Panoche and Carrizo populations. One was trained rangewide, including both populations. The last is the model of historical distribution based on aerial imagery (A.L. Rutrough, personal communication). The four models were divided into low and high suitability. Low suitability was assigned a score of one, high suitability a score of two, and the models were all added. Darker shading indicated higher suitability and more model agreement. All models are projected into 2070-2099 using CCSM4 rcp8.5. The grey shading indicates the historical distribution presented by Williams et al. (1992), and the crosshatching shows areas of agriculture or residential use. San Joaquin Valley, California Of particular interest is Cuyama Valley, just south of the Carrizo. Future suitability in the Cuyama Valley is high according to the historical, rangewide and local models (Figure 7). The Valley is currently used for ranching, agriculture and oil production, but is predicted to be suitable for  (Germano, 2010;Williams et al., 1993).
Conservation prioritization based on future suitability would improve the chance of survival for D. ingens populations. This may include a combination of future distribution predictions, connectivity modelling (i.e., Alexander, 2016) and strategic land conservation.
The protection of this keystone species and ecosystem engineer will help maintain the landscape for other taxa and help ensure the health of the community in the future.

ACK N OWLED G EM ENTS
We thank The Bureau of Land Management, Bureau of Reclamation

DATA ACC E S S I B I L I T Y
As the raw information represented in the figures pertains to a listed endangered species, it is not publicly available. However, raster layers of the models are available to managers and the scientific community upon request. She is interested in spatial ecology, climate change and endangered species with an emphasis on mammalian conservation.

William T. Bean is an Assistant Professor in the Department of
Wildlife at Humboldt State University. His research examines habitat needs across spatial scales in the past, present and future to inform conservation.