Functional connectivity and home range inferred at a microgeographic landscape genetics scale in a desert‐dwelling rodent

Abstract Gene flow in animals is limited or facilitated by different features within the landscape matrix they inhabit. The landscape representation in landscape genetics (LG) is traditionally modeled as resistance surfaces (RS), where novel optimization approaches are needed for assigning resistance values that adequately avoid subjectivity. Also, desert ecosystems and mammals are scarcely represented in LG studies. We addressed these issues by evaluating, at a microgeographic scale, the effect of landscape features on functional connectivity of the desert‐dwelling Dipodomys merriami. We characterized genetic diversity and structure with microsatellites loci, estimated home ranges and movement of individuals using telemetry—one of the first with rodents, generated a set of individual and composite environmental surfaces based on hypotheses of variables influencing movement, and assessed how these variables relate to individual‐based gene flow. Genetic diversity and structure results evidenced a family‐induced pattern driven by first‐order‐related individuals, notably determining landscape genetic inferences. The vegetation cover and soil resistance optimized surface (NDVI) were the best‐supported model and a significant predictor of individual genetic distance, followed by humidity and NDVI+humidity. Based on an accurate definition of thematic resolution, we also showed that vegetation is better represented as continuously (vs. categorically) distributed. Hence, with a nonsubjective optimization framework for RS and telemetry, we were able to describe that vegetation cover, soil texture, and climatic variables influence D. merriami's functional connectivity at a microgeographic scale, patterns we could further explain based on the home range, habitat use, and activity observed between sexes. We describe the relationship between environmental features and some aspects of D. merriami‘s behavior and physiology.


| INTRODUC TI ON
Since the term was coined in 2003, landscape genetics (LG) has been a constantly growing field that combines population genetics, landscape ecology, and spatial analytical techniques to quantify the effects that the landscape has on microevolutionary processes (Balkenhol, Cushman, Storfer, & Waits, 2015;Manel, Schwartz, Luikart, & Taberlet, 2003;Storfer, Murphy, Spear, Holderegger, & Waits, 2010). One such process is the movement of genes among populations (i.e., gene flow) that, in animal species, is primarily determined by the dispersal of individuals, at different scales, across the landscape (Reding, Cushman, Gosselink, & Clark, 2013). A variety of features within the landscape matrix may limit or facilitate the movement and gene flow of individuals (structural and functional connectivity). Landscape representation in LG has been traditionally modeled as resistance surfaces (RS), which can be defined as a spatial layer that assigns a value to each of a set of landscape features, and where values denote the degree to which a feature limits or facilitates connectivity across the landscape. Resistance surfaces can be viewed as hypotheses of the relationship between landscape variables and gene flow (Spear, Balkenhol, Fortin, McRae, & Scribner, 2010;Spear, Cushman, & McRae, 2015;Zeller, McGarigal, & Whiteley, 2012).
A common method for developing RS is the parameterization of resistance values based on expert opinion (Murray et al., 2009) to test alternative costs ratios (e.g., 2:1, 10:1, 100:1), in which one variable is hypothesized to always facilitate (e.g., habitat value = 1) and the other to restrict movement at different orders of magnitude (nonhabitat value = 2, 10, 100). This approach has been amply applied in natural systems (Howell, Delgado, & Scribner, 2017;Spear & Storfer, 2008. Simulation studies have demonstrated that the value of nonhabitat relative to habitat is key for detecting an effect on gene flow, in which a higher contrast in the costs ratios will make this relationship more detectable (Cushman, Shirk, & Landguth, 2013;Jaquiéry, Broquet, Hirzel, Yearsley, & Perrin, 2011). More recently, LG studies using this parameterization approach have also implemented mirror-like cost ratios, by assigning high cost values not only to nonhabitat but also to habitat patches (i.e., 10:1, 2:1, 1:2, 1:10), which eliminates the potential bias of only assigning increasing cost values to variables considered to restrict gene flow (see Hohnen et al., 2016). Notwithstanding, expert opinion RS development is based on arbitrary costs with no consensus about the values assigned for the cost ratios, frequently assuming a linear relationship between continuous variables and genetic distances, which may not always be the case . Because accurate inferences about these relationships are dependent on both, the correct representation of the landscape and of the values of landscape variables, we require methods for assigning resistance values that adequately avoid subjectivity (Richardson, Brady, Wang, & Spear, 2016;. A recent proposal is to simultaneously evaluate multiple surfaces and a wide range of cost values, without making any assumptions about their relationship with genetic distances. Specifically, Peterman, Connette, Semlitsch, and Eggert (2014) and Peterman (2018) used Ricker and monomolecular equations to transform data, in combination with linear mixed-effects models and genetic algorithms. This RS optimization framework does not make a priori assumptions regarding the scale and direction of the resistance relationship, allowing to perform both the optimization of categorical surfaces and the simultaneous optimization of multiple RS.
Landscape genetic studies have focused mostly on terrestrial animals and temperate forest ecosystems (Dyer, 2015;, whereas tropical and specially desert ecosystems are significantly lacking. Deserts represent one of the most extended ecosystems on Earth-nearly a third of the globe-exhibiting environmental characteristics like extreme temperatures and low precipitation regimes, which result in low net primary productivity and scarce vegetation cover (Whitford, 2002). For such reasons, desert ecosystems may seem to harbor little heterogeneity and, consequently, no significant role of landscape features in determining evolutionary processes like gene flow; the latter is likely associated with their underrepresentation in LG literature .
Rodents are considered a key ecological component of deserts due to the fundamental role they play for the structure and dynamics of these ecosystems, especially by the dispersal of seeds and soil removal (Brown & Heske, 1990a, 1990b. In addition, rodents have been proposed as ideal for conducting LG research (Waits, Cushman, & Spear, 2015) because of their small body size, short generation times, and limited dispersal abilities. Accordingly, LG studies have focused on different rodent species living in a wide range of environments, including tropical dry forests (spiny pocket mouse, Liomys pictus; Garrido-Garduño, Téllez-Valdés, Manel, & Vázquez-Domínguez, 2015), savannas of South Africa (Natal multimammate mouse, Mastomys natalensis; Russo, Sole, Barbato, von Bramann, & Bruford, 2016), subantarctic forests and Patagonian steppes (Longtailed pygmy rice rat, Oligoryzomys longicaudatus; Ortiz et al., 2017), and even in urbanized areas (White-footed mouse, Peromyscus leucopus, Munshi-South, 2012;and Norway rats, Rattus norvegicus, Gardner-Santana et al., 2009). Interestingly, vegetation has been underlined as a variable facilitating gene flow in some of those systems, including the canopy cover across New York City for P. leucopus (Munshi-South, 2012), tropical dry forest corridors in L. pictus (Garrido-Garduño et al., 2015), and forest patches in chipmunks (Tamias striatus) inhabiting agroecosystems (Kierepka, Anderson,  One of the most conspicuous desert rodents is kangaroo rats (genus Dipodomys, family Heteromyidae). In particular, the Merriam's kangaroo rat Dipodomys merriami is one of the smallest species in the genus, with a distribution that encompasses the desert region of southwestern United States and northern Mexico. It is a burrowdwelling, nocturnal rodent characterized by a long tail with a tuft of hair at the tip, large rear legs for bipedal locomotion, and with sexual dimorphism where males are larger than females. It is considered solitary and territorial, characterized by dispersing adult males (mean of 60 m) and phylopatric females that make parental investment (Behrends, Daly, & Wilson, 1986;Randall, 1993;Zeng & Brown, 1987). Their mating system is polygynandrous, where matings occur mainly between close neighbors (Randall, 1993) during February and July. They show a mean survival of 3.5 years, with two litters per year and an average of three young per litter (Zeng & Brown, 1987), feeding commonly on mesquite seeds (Prosopis sp.) and always building their burrows under bushes like mesquite and creosote (Larrea tridentata) (Murrieta-Galindo & Cuatle-García, 2016;Reynolds, 1958).
As a desert organism, D. merriami is challenged by the low productivity, extreme temperatures, and low precipitation regimes of its habitat. Indeed, D. merriami has shown to be highly sensitive to food shortage under controlled conditions (Banta, 2003), in agreement with a desertic environment where food is scarce, and thus rendering mesquite (and other desert plants) key for its survival (Reynolds, 1958). Moreover, D. merriami's daily activity (e.g., foraging, dispersal, breeding) is affected by ambient temperature (29-34ºC;French, 1993;Banta, 2003) and relative humidity (40%; Reynolds, 1958;Frank, 1988;Walsberg, 2000). Although some information about the natural history and ecology of D. merriami exists (e.g., Behrends et al., 1986;Zeng & Brown, 1987;Randall, 1993), the effect of landscape and environmental features on genetic structure and connectivity has not been evaluated.
The space where an animal species performs all its activities (e.g., feeding, reproduction, dispersal, avoiding predators) defines what is known as home range (Burt, 1943); a home range represents an interplay between the environment and an animal's understanding of that environment (i.e., its cognitive map) (Powell & Mitchell, 2012).
The selection of a place to live among different alternatives available (habitat selection), tightly linked to the home range, is a hierarchical process that involves both innate and learned behaviors (Partridge, 1978). Recent advances in telemetry, geographical positioning systems, and home range estimators have facilitated the description of habitat selection patterns and home ranges in wild species, at diverse micro-and macroenvironmental scales; allowing, in turn, the understanding of organismal movement, how and why animals use specific resources, and other elements that determine fitness and survival (Kernohan, Gitzen, & Millspaugh, 2001;Demšar et al., 2015). Importantly, combining knowledge on habitat use, behavior, and population genetics can significantly contribute to the interpretation of functional connectivity in landscape genetics (Portanier et al., 2018).
We here apply the optimization framework developed by Peterman et al. (2014) and Peterman (2018) to evaluate, at a microgeographic scale, the effect of landscape features on gene flow of Dipodomys merriami in northern Mexico. For that purpose, we first characterized genetic diversity levels and population genetic structure; next, we generated a set of environmental surfaces based on hypotheses of variables influencing the species movement, to finally assess how these variables relate to individual-based genetic differentiation. In order to obtain estimates of the species habitat use and home range that will inform our genetic results, we evaluated movement of individuals with telemetry. Given that D. merriami is tightly linked to shrub cover and occurrence of fine-sized soils for protection against predators, feeding, and burrow construction, we expected that (a) females and males will mostly overlap home ranges; (b) the presence of vegetation associated with sandy and gravel soils would be a key predictor of gene flow, while uncovered habitat (i.e., bare soil and/or rocky landscapes) would restrict movement; (c) environmental variables (e.g., temperature, humidity) will significantly influence patterns of functional connectivity in this desert-dwelling rodent.

| Study site and sampling
The Chihuahuan Desert is one of the largest deserts in North America, known as the Bolsón de Mapimí in Mexico, where our study was conducted at the Mapimí Biosphere Reserve (MBR), a protected natural area and a Man and Biosphere (MAB-UNESCO) site (Montaña, 1988) (Figure 1). The MBR extends 342,387 ha, representative of the arid ecosystems from northern Mexico that harbors a high endemism and species richness (Conanp, 2006). The climate is dry (mean annual precipitation: 271 mm, mean annual temperature: 20.8°C, rainy season from July to October and elevation range: 1,100-1,650 m; Montaña & Breimer, 1988;Conanp, 2006). The vegetation is dominated by shrubby species like mesquite (Prosopis glandulosa) and creosote bush (Larrea tridentata). Despite the zone is relatively flat, the landscape has been described as a toposequence, on the base of landforms, soils and vegetation units organized along an elevation gradient (Grünberger, 2004;Montaña & Breimer, 1988) and vegetation units: "magueyal" (predominating Agave sp. and L. tridentata), "nopalera" (Opuntia sp.), and "pastizal" (Prosopis sp.) (Martínez & Morello, 1977) (Figure 1). Sampling was performed from 31 May to 11 June 2015, using a modification of the trapping web sampling method (Anderson, Burnham, White, & Otis, 1983) in order to include the most environmental heterogeneity across the scale of the study; importantly, this trapping web has proven to be effective for capturing D. merriami in the MBR (Aragón, Castillo, & Garza, 2002;Hernández et al., 2005). Briefly, the web design consists of  Figure S1 in the Appendix S1). Trapped individuals were sexed and measured; a tissue sample (earpunches) was taken for genetic analysis and stored in labeled Eppendorf tubes with 96% ethanol. All individuals were released at the sampling site. Procedures were conducted according to the American Society of Mammalogists guidelines for use of wild mammal species (Sikes, 2016) and with the corresponding collecting permits.

| Population genetic analyses
We extracted DNA from tissue samples to test 14 fluorescently labeled microsatellite primers developed for D. spectabilis, from which we successfully amplified nine with D. merriami (see the Supporting Information Appendix S1 and Table S1 for more information).
We tested for deviation from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) across the entire population Genepop v.6 (Rousset, 2008). The value was adjusted for multiple comparisons using a Bonferroni correction (Rice, 1989). We checked for the presence of null alleles and stuttering with Microchecker (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004). Genetic diversity levels were assessed by estimating number of alleles, number of effective alleles, unbiased expected and observed heterozygosity, and F IS with gstudio in R (Dyer, 2014), while relatedness between individuals was assessed with ML-Relate (Kalinowski, Wagner, & Taper, 2006).
We inferred population structure with two approaches, the Bayesian clustering method implemented in Structure v.2.3.4 (Pritchard, Stephens, & Donnelly, 2000) and the discriminant analysis of principal components (DAPC; Jombart, Devillard, & Balloux, 2010). For Structure, we tested from K = 1 to K = 13 genetically distinct clusters with 20 replicates for each K, based on an initial burn-in of 2 x 10 5 followed by 5 x 10 5 Monte Carlo Markov chain iterations and an admixture model with correlated allele frequencies. In addition, we incorporated the sampling information (locprior option) to the model by considering each trapping web as a different location. The posterior probability lnP(K) (Pritchard et al., 2000) and ∆K (Evanno, Regnaut, & Goudet, 2005), implemented in Structure Harvester (Earl & vonHold, 2012), were used for identifying the most likely number of clusters. Although both methods are suitable for analyzing spatially continuous data, the Structure algorithm considers population genetics models (Pritchard et al., 2000). Because some natural systems can violate some model-assumptions which in turn may lead to errors in the assignment of individuals to populations, we performed an additional multivariate, nonmodel-based approach (DAPC), using adegenet 2.1.0 in R (Jombart et al., 2010), following a two-stage procedure: First, a principal component analysis (PCA) is performed with the genetic data; then, the principal components (PCs) of the PCA are processed with a linear discriminant analysis (LDA). Since DAPC relies on the determination of predefined groups that are often unknown, they must be identified a priori. Thus, we first inferred genetic clusters by running the K-means clustering algorithm, from K = 1 to K = 13, with find.clusters. For this step, we retained all 76 PCs based on the recommendation of keeping all (or at least 80%) of the information (Jombart & Collins, 2015). Next, based on the lowest Bayesian information criterion (BIC) value, we determined the optimal number of clusters, with which we ran the dapc function. Unlike the K-means clustering algorithm, DAPC benefits from not using too many PCs. Finally, to determine the number of PCs and also avoid overfitting during discrimination, we used the optim.a.score and xvalDapc functions; the scatter function was used to build an ordination plot with the results. Finally, we also estimated pairwise F ST values between genetic clusters inferred by each method with hierfstat 0.04-22 in R (Goudet, 2015).
Considering that evidence of population structure can be found when family members are included in a sample, as would be the case for D. merriami, even when such structure is absent (Anderson & Dunham, 2008), we identified first-order relatives (full siblings, FS; and parent-offspring, PO) based on the previous relatedness analysis (Kalinowski et al., 2006). Next, we removed one individual of each dyad and ran both analyses again with the same parameters for this unrelated dataset. For DAPC analysis, we retained 59 PCs in the find.clusters function with this dataset (see Results). This procedure has shown to improve accuracy on estimates of genetic structure (Rodríguez-Ramilo & Wang, 2012), particularly when conducting landscape genetics studies (Peterman et al., 2014;Ruiz-Lopez et al., 2016).
Finally, genetic dissimilarity at the individual level (i.e., between all pairs of individuals) was estimated as the proportion of shared alleles (D PS ) with adegenet in R (Jombart et al., 2016). We chose D PS because it makes no biological assumptions and can be used for populations at any level of ploidy or inbreeding; it has also proved to be an adequate metric for performing individual-based genetic distance estimates (Shirk, Landguth, & Cushman, 2017). To test for isolation by distance, we performed a simple Mantel test between genetic (D PS ) and Euclidean distances. Additionally, a Mantel correlogram was estimated based on 50 m classes considering the mean dispersal distance reported for D. merriami (Zeng & Brown, 1987) and our own home range results. Spearman correlation significance for both analyses was based on 10,000 permutations. Euclidean distances were estimated with gstudio (Dyer, 2014), while the Mantel test and correlogram were calculated with ecodist (Goslee & Urban, 2007), both in R.

| Landscape data
Since most of the freely available environmental data is at coarse scales (e.g., WorldClim ca. 1 km 2 resolution; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), we generated a set of environmental surfaces at our study spatial scale which were hypothesized to affect survival or movement of D. merriami: humidity (as a proxy for precipitation at a very local scale), temperature, elevation, plant cover (hereafter vegetation), and soil. Humidity and temperature data were collected using HOBO Data Loggers (Pro v2 and UX100-003, ONSET Computer Corporation), and elevation was obtained with a GPS (Garmin). These three variables were measured at each trap location and at the center of the trapping web (i.e., 21 points per trapping web and 105 by transect), for a total of 315 points; the values for the entire study area were obtained by the krigging interpolation method (Holdaway, 1996)  Using lower-resolution imagery to characterize land cover can lead to incorrect or misleading evaluations of connectivity if not verified on the field (Zeller, Nijhawan, Salom-Pérez, Potosme, & Hines, 2011), and on the other hand, some satellite bands are not available in Google Earth imagery, thus lacking information about certain landscape features (Boyle et al., 2014). Accordingly, due to the fine scale of this study, we chose to use both Landsat 8 (lower-resolution) and Google Earth (higher-resolution) data to more adequately discern vegetation and soil variables. All surfaces were processed with R v.3.3.2 (R Core Team, 2016) and resampled to a resolution of 5 m for landscape genetics analysis.

| Landscape genetics analyses
We followed the optimization framework developed by Peterman et al. (2014) to determine the resistance values of our surfaces.
Briefly, this approach uses monomolecular and Ricker functions (Bolker, 2008) to transform continuous resistance surfaces; it relies on a genetic algorithm (GA; Scrucca, 2013) that adaptively explores the parameter space, seeking to maximize the relationship between pairwise landscape distances (least-cost or resistance) and pairwise genetic distances, while making no a priori assumptions about their relationships. The monomolecular [y = r (1-exp -bx )] and Ricker [y = r exp -bx ] are two exponential-based functions used for ecological modeling, differing in the curve shape of the relationship they are modeling. This curve shape is mainly determined by shape (x) and magnitude (b) parameters, which produce a saturating exponential (growth or decay) curve for the monomolecular function, and a hump-shaped curve (skewed to right or left) for the Ricker function (Bolker, 2008). During the optimization process, the genetic algorithm searches all possible combinations of these parameters for transforming resistance surfaces, denoted by "r" in the monomolecular and Ricker equations (Peterman, 2018;Peterman et al., 2014).

The optimization framework was performed with
ResistanceGA in R (Peterman, 2018; https://github.com/wpeter man/ResistanceGA) following two steps: First, all surfaces were independently optimized from pairwise resistance distances estimated with gdistance in R (Van Etten, 2017), with the commu-teDistance function, by exploring resistance values up to 2,500.
Previous studies using this optimization approach (e.g., Peterman et al., 2014, Ruiz-Lopez et al., 2016, Khimoun et al., 2017 have measured resistance distance using circuitScape (McRae, 2006); however, it is known that commuteDistance is functionally equivalent to circuitScape, with the advantage that it can be run in parallel (Kivimäki, Shimbo, & Saerens, 2014;Peterman, 2018). All these processes were performed using an eight-neighbor connection scheme for assessing connectivity. We conducted three independent optimization runs for each surface to confirm convergence and parameter estimates (see the Appendix S1 for the R optimization scripts).
We used AIC as our objective function during optimization, which was determined from linear mixed-effects models (lmem). The lmem were fitted by the maximum-likelihood population effects (MLPE) parameterization, to account for the nonindependence of values within pairwise distance matrices (Clarke, Rothery, & Raybould, 2002;Van-Strien, Keller, & Holderegger, 2012). The dependent and predictor variables were pairwise genetic distance (D PS ) and pairwise scaled and centered resistance distance, respectively. MLPE parameterization was done with lme4 (Bates, Maechler, Bolker, & Walker, 2014) in R, and support of the optimized resistance surfaces was assessed using the AICc (Akaike's information criterion corrected for small/finite sample size; Akaike, 1974). To evaluate the robustness of our model selection and optimization given different combinations of samples, we performed a bootstrap resampling of the data (Peterman et al., 2014;Ruiz-Lopez et al., 2016). Next, to control for potential bias in our results, 75% of the samples were randomly selected without replacement and each surface was then fit to the subset of individuals; the average rank, average model weight, and the percentage that a surface was selected as the best (top ranked) model following 10 000 iterations were estimated.
Before the second optimization step, we did a Spearman coefficient correlation test (rho =ρ) with R between the surfaces that showed a greater selection percentage than distance alone, and selected, based on Cohen (1992), the surfaces that showed a small to medium correlation (ρ < 0.29), in order to avoid correlated variables in the multisurface model. Finally, we performed a multisurface optimization for selected surfaces to generate composite surfaces.
Bootstrap model selection was performed again (75% of samples and 10,000 iterations) to obtain the average rank, average model weight, and the top ranked model of individual (i.e., univariate) and composite surfaces.

| Radiotracking
We selected 17 individuals from three different webs along one transect (Figure 1), where trapping success was highest; each individual was equipped with a TXB-003G radiotransmitter (Telenax, Mexico) attached to the suprascapular area with a drop of instant-dry glue.
The radiotransmitters weighed ca. 0.6 g, representing 1.3%-2% of the animals' body weight, which is below the 5% maximum proposed by White and Garrot (1990). Radiolocations were taken one day after the individuals had been released (Springer, 2003), using a handheld three-element Yagi directional antenna and an RX-TLNX receiver (Telenax), between 21:00 and 01:30 hr, with intervals of 30 min and never recording the same individual consecutively, assuring data independence (Kernohan et al., 2001). Radiolocations were taken every night until reaching 10 per individual; georeferenced data, time, associated vegetation, and if the individual was actually observed were recorded for each one.
We used a probability-based statistical estimator, the Kernel estimator, to calculate home range size (kernel estimating functions) with adehabitatHR (Calenge, 2015) in R v.3.3.0, using the smoothing parameters h 1CSV and h ref ; area estimation did not differ between them thus we report the latter because it exhibited a better resolution for the isopleths that delimit regions with different probability within home ranges. This method estimates the intensity of area use as a two-dimensional relative frequency distribution of an animal's location over time (Worton, 1989), avoiding biases due to its low sensitivity for extreme data (Rodgers & Carr, 1998) that other methods have (e.g., the minimum convex polygon). We tested for potential differences between female and male home range sizes using nonparametric Mann-Whitney U tests and estimated their overlap percentage, with R. We also estimated the Euclidian distance between radiolocations tracked consecutively for each individual, and based on the time of each radiolocation, we evaluated activity peaks considering the maximum distances the individuals moved during each radiotracking night. Finally, we recorded an indirect estimate of habitat use as the presence (percentage) of a vegetation type per radiolocation; differences between sexes were analyzed with Mann-Whitney U tests.

| Population genetics
A total of 76 individuals were captured over the study area, 18 individuals in transect 1 (T1), 37 in T2, and 21 in T3. Eight microsatellite loci were polymorphic across all samples (Supporting Information Table S1 in the Appendix S1), whereas locus DS109 was monomorphic and was excluded from the analyses. Five loci (Ds1, Ds3, Ds19, Ds46, and DS98) deviated significantly from HWE after Bonferroni correction, while there was no evidence of LD (Bonferroni corrected p < 0.05); these same loci showed evidence of null alleles but not stuttering errors were detected. Because D. merriami's social structure encompasses groups of individuals with different degrees of relatedness (see relatedness results below), some evidence of deviation from HWE or null alleles is expected (i.e., not resulting from a systematic nonamplification of alleles; Bergl & Vigilant, 2007; Mapelli, Mora, Mirol, & Kittlein, 2012). Hence, all eight loci were included in the following analyses. Notably, we amplified a locus (Ds46) previously reported as unsuccessful in D. merriami, and we found that Ds19 is not a X-linked locus in this species (Davis et al., 2000).
Genetic diversity results showed that the number of alleles per locus ranged from 8 to 30 (mean=15.1) and the number of effective alleles from 3.2 to 17.7 (mean=8.16). The mean observed heterozygosity and the unbiased expected heterozygosity across loci were 0.65 and 0.84, respectively; overall F IS was 0.22 (Supporting   Information Table S2 in the Appendix S1). Regarding relatedness, we detected 27 first-order relationships (full siblings and parentoffspring), 11 occurred in different web/transect, six occurred in different webs along the same transect, and 10 occurred in the same web/transect. After removing one individual from each dyad, we obtained a dataset with 59 individuals (hereafter, unrelated dataset), used to test the effect of close relatives in population genetic structure.
Results about population genetic structure using Structure and DAPC with both datasets were equivocal. Five clusters were detected by Structure for the full dataset (n = 76), with both statistics (mean lnP(K) = −2,628.435 and ∆K = 12.331) (Supporting Information Figure S2a in the Appendix S1). However, the genetic clusters showed pairwise F ST values of 0.027 to 0.239 and had no congruent geographic spatial pattern (Supporting Information Figure S2b in the Appendix S1). On the other hand, the curve of Bayesian information criterion (BIC) values for the DAPC results suggested the presence of K = 1 to K = 3 genetic clusters; notably, BIC values decreased rapidly, reaching their lowest value at K = 2 before rising again (Supporting Information Figure S3 in the Appendix S1). Because multiple possible K values have been described as a characteristic scenario for continuously distributed species using this method (Jombart, 2008), we performed DAPC with each the K = 2 and K = 3 assigned clusters by retaining 20 PCs, which comprised 71.7% of the total variance (Supporting Information Figure S3a in the Appendix S1). In both cases, groups were discriminated but in which clusters are not spatially structured; also, F ST values between clusters were low (0.04 for K = 2; 0.042-0.059 for K = 3).
Regarding the unrelated dataset (n = 59), one genetic cluster was obtained with Structure (K = 1) by the mean lnP(K) = −2,195.545, whereas 11 clusters were detected with ∆K (20.781) (Supporting Information Figure S2c). For the DAPC, we retained all 59 PCs for the first step in the function find.clusters (Supporting Information Figure S4a in the Appendix S1). Results with DAPC were the same as those obtained with the full dataset; hence, we also applied DAPC to K = 2 and K = 3 (retaining 12 PCs and comprising 56.5% of the total variance), obtaining differentiation between groups but with clusters overlapping spatially (Supporting Information Figure S4b,c in the Appendix S1); F ST values were again markedly low (0.043 for K = 2, 0.05-0.073 for K = 3).
Differences between Structure and DAPC for detecting genetic structure in our study system can be explained by their assumptions.
DAPC requires predefined groups, making this decision crucial for downstream interpretation of genetic data. In this context, clusters can be visualized as tools to summarize and understand the data, but recognizing that complex systems are not always subject to this clear-cut representation (Jombart & Collins, 2015). As mentioned above, obtaining multiple values of K using the K-means algorithm has been related to continuously distributed species, but only on a wide geographic extent (Guerrero et al., 2018), hence performing better for island-based models than for continuous models (Jombart & Collins, 2015). Moreover, this algorithm uses a simple measure of group differentiation and is likely to fail to identify the correct number of clusters in complex population models (Jombart et al., 2010).
On the other hand, discrepancies between the Structure results for the full and the unrelated datasets, based on ∆K, agree with a family-induced population structure (Anderson & Dunham, 2008); in addition, ∆K does not perform adequately when K = 1 (Evanno et al., 2005). Consequently, we consider results support one genetic cluster (no structuring) based on the above analytical issues, the low F ST values between suggested clusters, the fine-spatial scale of the study area, and importantly, the social structure of D. merriami.

| Landscape genetics and connectivity
Genetic distance (D PS ) between individuals for LG analyses was estimated based on the unrelated dataset (Supporting Information Figure S5 in the Appendix S1). The Mantel test showed a nonsignificant negative correlation between the genetic and Euclidean distances (r = −0.022; p > 0.3), while the Mantel correlogram exhibited only one small positive and significant value at 1,700 m (r = 0.073; p < 0.05) (Supporting Information Figure S6 in the Appendix S1). Optimization and model selection results showed that the vegetation cover and soil resistance surface obtained with the NDVI were the best-supported model (54.8% of the times based on 10,000 bootstrap replicates; Table 1), with an inverse monomolecular function (Figure 2a). Furthermore, NDVI was a significant predictor of genetic distance on the generalized linear mixed-effects model (Supporting Information Table S3 in the Appendix S1). The NDVI optimized surface assigned high resistance to areas of the landscape with predominantly bare soil (<0.1), with a fast decrease in resistance where vegetation is present (> 0.1) (Figure 2a).
The following best-supported functional form was humidity (inverse-reverse Ricker), with different functional forms for temperature (inverse Ricker) and elevation (inverse-reverse monomolecular), in which resistance values were lowest around 44% humidity, 30ºC temperature, and 1,152 m (Figure 2b,c,d).
Additionally, humidity, elevation, and temperature resistance surfaces explained 30%, 4.2%, and 4.9%, respectively, of the variation in the pairwise genetic data than distance alone (Table 1) Table S4 in the Appendix S1). Accordingly, we constructed composite surfaces with these three comparisons, and also one including all four layers. Results of the bootstrap model selection showed that NDVI was the best-supported model (54.9%), followed by humidity (19.7%) and NDVI+humidity (10.7%) ( Table 2; the contribution of each variable to each multisurface model is shown in Table 3). The surface including all four layers (Combination 4) had no support and performed poorly in the generalized linear mixed-effects model (Supporting Information Table S5 in the Appendix S1).

| Home range
Of the 17 individuals with radiotransmitter, we lost radio signal for eight before obtaining enough data for analyses; hence, we report the results for nine individuals (two females, six males, one juvenile, sex undetermined) (Table 4). We obtained 91 independent radiolocations, and six individuals were directly observed while moving.
The majority of the radiolocations used to calculate home range size per individual are within the estimated area, which shows these are regular activity zones. The home range size estimated was on average 0.6294±0.264 ha, 0.6957±0.3770 ha for males and 0.2453 ha for females (Figure 3, Table 4), with significant differences between sexes (U = 45, p = 0.003). The largest home range was a male's (R8), completely overlapping with that of another male (R11) and a female (R13). A 77.4% home range overlap was observed between a male and a female (radios R7 and R15), while among R8, R10, R11, and R13 varied from 10.3% to 100% (Figure 3, Table 4). Only two of the nine individuals showed a second-order relationship (half-siblings), the R7 male and the R15 female.
Two activity cycles were recorded for the nine D. merriami individuals with radiolocations, during 8 days, one from 21:30 to 00:55, and the other from 23:00 to 01:10 ( Figure 4)

| Genetic diversity and microgeographic familyinduced structure
The Considering that trapping was conducted during the breeding season, our results may reflect movements of individuals across their home range area (interpopulation dispersal), generating a nonequilibrium pattern evidenced by heterozygote deficiency, a phenomenon observed for other Dipodomys species like the banner-tailed kangaroo rat D. spectabilis (Busch et al., 2007). In this context, D.
merriami exhibits a mean distance dispersal of around 60 m, while it has been suggested it can move longer distances during the breeding season with no distinction between sexes (Behrends et al., 1986;Zeng & Brown, 1987). Our results are in agreement, considering that the longest dispersal between radiolocations we observed was 57-61 m, performed by males. On the other hand, significant home range size differences between sexes were found (0.695 ha for males and 0.245 ha for females), despite being the breeding season; the latter can be related to less female movement compared with males due to gestation and offspring care (Murrieta-Galindo & Cuatle-García, 2016;Nader, 1978). Also, our findings show that females overlap home ranges with males but not between them; only two radiotracked individuals were related (half-siblings), a male and a female that greatly overlapped their home ranges (70%). We acknowledge our low sample size and that results need be taken with caution; nonetheless, results are congruent with the behavior of D.
merriami and support our prediction. Particularly, this species is considered a solitary and territorial rodent, characterized by adult male dispersal and phylopatric females that perform parental investment (Behrends et al., 1986;Randall, 1993).
Our results showed both no effect of geographic distance on genetic differentiation of D. merriami and no genetic structuring within the scale studied. Additionally, D. merriami exhibits frequent burrow shifts to avoid predators (Behrends et al., 1986), which has been related to a long-life span (up to four years; Zeng & Brown, 1987), and consequently, a long-term stability of populations that can result in a social phenomenon of tolerance by individual familiarity and mate  (Behrends et al., 1986), leading to mating events between close neighbors (Randall, 1993). Hence, our overall F IS results and lack of population genetic structure agree with a family-induced pattern driven by first-order-related individuals, a biologically meaningful aspect of D. merriami described by ecological and behavioral data (Randall, 1993;Zeng & Brown, 1987), and now supported by the genetic component. Our study evidences the implications that a family-induced structure at a microgeographic scale can have on landscape genetic inferences, specifically for selecting the unit level (i.e., populations or individuals) and genetic differentiation metrics for analyses (Shirk et al., 2017). Indeed, a correct assessment of population structure should always be conducted keeping in mind the more biologically relevant patterns (Anderson & Dunham, 2008;Bergl & Vigilant, 2007;Rodríguez-Ramilo & Wang, 2012;Ruiz-Lopez et al., 2016).

| Environmental features and genetic connectivity
Our findings that gene flow in D. merriami is best explained by the normalized difference vegetation index (NDVI) enabled us to suggest a link between the observed genetic pattern and the ecological processes underlying it (e.g., the species' dispersal, foraging, physiology). Finding these mechanistic links in landscape genetics research is one of the goals of the optimization framework we implemented (Peterman, 2018). The simultaneous optimization of multiple resistance surfaces (RS) and our creating a vegetation cover and soil texture surface based on the NDVI allowed us to explain the variation in genetic data better than any other individual or composite surface, suggesting these features are biologically relevant for D. merriami (Peterman et al., 2014;Ruiz-Lopez et al., 2016;. Although it has been suggested that multiple RS should be used for  , results do vary; for instance, Ruiz-Lopez et al. (2016) found that a composite surface comprised of fire density and the distance to the nearest village describes the variation in genetic data of the red colobus monkey (Procolobus gordonorum), suggesting a strong influence of anthropogenic activities on the species movement. On the other hand, Khimoun et al. (2017) showed that individual optimization of land cover RS has a higher support compared to composite surfaces in the insular tropical bird Plumbeous warbler (Setophaga plumbea).
According to our expectations, we identified that vegetation tems (Kierepka et al., 2016). Moreover, D. merriami activity peaks exhibit a behavior tightly associated with predator avoidance (Daly, Behrends, Wilson, & Jacobs, 1992;Soltz-Herman & Valone, 2000), where movement is limited not only on open areas but also by lunar light (Daly et al., 1992;Fuentes-Montemayor et al., 2009). Indeed, as our activity cycle results show, D. merriami is more active at crepuscular hours during full moon, while as the moonlight decreases its activity increases at midnight.
Interestingly, despite studies indicate that D. merriami feeds preferentially on seeds of Prosopis glandulosa and builds its burrows under this mesquite (Cox, De Alba-Avila, Rice, & Cox, 1993;Reynolds, 1958), the dominant vegetation type on our study region is creosote (Larrea tridentata). Indeed, D. merriami exhibits a wider habitat use that includes preferentially L. tridentata, but also Opuntia sp., J. dioica, and P. glandulosa, highlighting the key role of this generalist species on desert ecosystem dynamics (Brown & Heske, 1990b;Murrieta-Galindo & Cuatle-García, 2016). Soil type, in this case sandy and gravel soils, is also a key factor associated with the plant species present, burrow construction, and individual movement within D. merriami's home range.
Most landscape genetics inference studies have assumed either a positive or negative linear relationship between landscape features and cost surfaces (Garroway, Bowman, & Wilson, 2011;Koen, Bowman, & Walpole, 2012), including examples with rodent populations (Chiappero et al., 2016;Howell et al., 2017;Mora et al., 2017;Ortiz et al., 2017), despite that nonlinear responses are expected to be more common (Marrotte & Bowman, 2017;. Here, we found nonlinear relationships relative to landscape resistance for both humidity and temperature, variables that influence genetic connectivity in different species like the northern quoll (Dasyurus hallucatus; Hohnen et al., 2016) and Liomys pictus (Garrido-Garduño et al., 2015). Temperature has been explicitly proposed as determinant for genetic connectivity in climate sensitive species, for example, the American pika (Ochotona princeps), a heat intolerant, and cool microclimate restricted species for which an increase in temperature adversely affects gene flow (Castillo, Epps, Davis, & Cushman, 2014). Given the microgeographic scale of our study compared to the above mentioned, the nonlinear relationships we find may be associated with some aspects of D. merriami's physiology, directly related to connectivity across the landscape. The thermoneutral zone (TNZ) refers to the temperature gradient where an organism's metabolism is minimized but also leads to higher rates of water loss, which ranges between 29º and 34ºC for this species (French, 1993). As a desert-dwelling mammal, D. merriami has evolved certain physiological traits, for instance during its nocturnal active phase, it selects cooler TNZ temperatures (30.3-31.5ºC) for water conservation (Banta, 2003). Also, it is known that 40% humidity favors its movement on the surface (i.e., not across their burrows underground), whereas values outside this range may be detrimental (Frank, 1988;Reynolds, 1958;Walsberg, 2000

| Thematic resolution for detecting landscape genetic patterns
Scale has been recognized as a central question in ecology (Levin, 1992 and genetic data, as highlighted by Cushman and Landguth (2010).
In addition, the two approaches we used to represent the same feature at different scale attributes partially allowed us to exhibit their effects at a microgeographic scale, where vegetation was better represented as continuously distributed, thus evidencing the need to characterize that feature according to its "more-real" nature. Finally, there are recent advances aimed to make categorical data more ecologically relevant in LG research (Peterman, 2018), although they still need to be evaluated empirically. Our study is an empirical example of how the pixel size and thematic resolution need be considered, particularly for small body size species and at microgeographic areas, when developing RS that are relevant to the study system.

| CON CLUS IONS
Evaluating the effects of landscape features (landscape matrix) on individual movement and gene flow in natural populations has been challenged by the need to avoid subjectivity when assigning resistance values and landscape scales (Khimoun et al., 2017;Richardson et al., 2016;. Our study is novel in diverse aspects, where we present a landscape genetics study with a desert-dwelling rodent species-an ecosystem rarely investigated under this genetics approach-using a nonsubjective optimization framework for resistance surfaces and an accurate definition of thematic resolution (Khimoun et al., 2017;Peterman, 2018;Peterman et al., 2014 Based on our overall results, we describe patterns of the relationship between environmental features and some aspects of the behavior and physiology of D. merriami.

ACK N OWLED G M ENTS
We are grateful with L. León-Paniagua and J. Golubov for discussions throughout the entire project, and with Jane Remfert and

CO N FLI C T S O F I NTE R E S T
The authors declare that they have no conflicts of interest.