Fine-scale analysis reveals cryptic landscape genetic structure in desert tortoises.

Characterizing the effects of landscape features on genetic variation is essential for understanding how landscapes shape patterns of gene flow and spatial genetic structure of populations. Most landscape genetics studies have focused on patterns of gene flow at a regional scale. However, the genetic structure of populations at a local scale may be influenced by a unique suite of landscape variables that have little bearing on connectivity patterns observed at broader spatial scales. We investigated fine-scale spatial patterns of genetic variation and gene flow in relation to features of the landscape in desert tortoise (Gopherus agassizii), using 859 tortoises genotyped at 16 microsatellite loci with associated data on geographic location, sex, elevation, slope, and soil type, and spatial relationship to putative barriers (power lines, roads). We used spatially explicit and non-explicit Bayesian clustering algorithms to partition the sample into discrete clusters, and characterize the relationships between genetic distance and ecological variables to identify factors with the greatest influence on gene flow at a local scale. Desert tortoises exhibit weak genetic structure at a local scale, and we identified two subpopulations across the study area. Although genetic differentiation between the subpopulations was low, our landscape genetic analysis identified both natural (slope) and anthropogenic (roads) landscape variables that have significantly influenced gene flow within this local population. We show that desert tortoise movements at a local scale are influenced by features of the landscape, and that these features are different than those that influence gene flow at larger scales. Our findings are important for desert tortoise conservation and management, particularly in light of recent translocation efforts in the region. More generally, our results indicate that recent landscape changes can affect gene flow at a local scale and that their effects can be detected almost immediately.

Introduction

Features of the landscape can either restrict or promote movements of individuals in natural populations and consequently influence the degree of genetic connectivity among populations. Characterizing the effects of landscape features on genetic variation is essential for understanding how landscapes shape patterns of gene flow and spatial genetic structure of populations, and this objective forms the basis of landscape genetics [1], [2]. Novel approaches in landscape genetics have yielded new insights about how geographic and environmental features structure genetic variation, especially at fine spatial scales, insights that have had implications for ecology, evolution, and conservation biology.

Detecting cryptic population genetic structure and identifying features of the landscape that influence gene flow are central to understanding local- and regional-scale movements and barriers. Decomposing a genetic mixture of samples into their component, panmictic subpopulations (with or without spatial location information) helps to identify genetic discontinuities that could restrict gene flow among populations (e.g. [3], [4], [5], [6]). Both natural (e.g., rivers, mountains; [7], [8], [9], [10] and anthropogenic (e.g., roads, urban areas; [8], [11], [12] barriers that correlate with genetic discontinuities have been identified within landscapes, and recent studies have distinguished among alternative barrier scenarios [13], [14]. Understanding patterns of movement throughout a landscape and environmental features that influence those movements has been particularly useful in population conservation and management, to identify and prioritize populations for action and to design effective management strategies.

Most landscape genetics studies focus on patterns of movements of individuals at a regional scale. For example, some of the best-known studies in landscape genetics of vertebrates (Columbia spotted frogs, [7]; wolverines, [15]; roe deer, [8]; and black bears, [13]) have characterized patterns of gene flow and genetic structure over several thousand to several hundred thousand square kilometers. These studies have provided unique insight into large-scale landscape features and environmental conditions that influence gene flow at a regional scale. However, organisms can also respond to the landscape at a smaller scale, perceiving slight changes in the local landscape and making movement decisions based on these perceptions. Although it is unclear whether factors influencing the movement of individuals at fine spatial scales may be important predictors of genetic structure at broader scales, recent research suggests that they may not [16]. Thus, the genetic structure of populations at a local scale may be influenced by a unique suite of landscape variables that have little bearing on connectivity patterns observed at broader spatial scales. These local effects may be particularly pronounced in harsh environments, where habitat suitability varies considerably over small distances [17], or for organisms that see habitats as coarse- versus fine-grained.

We carried out a landscape genetics study on the desert tortoise (Gopherus agassizii) to investigate fine-scale spatial patterns of genetic variation and gene flow in relation to features of the landscape. The desert tortoise occupies portions of three major North American deserts (Sonoran, Mojave, and Colorado; though the Sonoran desert tortoise has recently been proposed as a new species, Gopherus morafkai, [18]), and is federally listed as threatened in the northwestern one-third of its geographic range (the ‘Mojave population’, including all tortoises north and west of the Colorado River; [19]). At a broad scale, desert tortoises are separated by the Colorado River, and long-term restricted gene flow between these regions has led to substantial morphological and genetic differences between Sonoran and Mojave desert tortoises [20], [21], [22], [23], [24].

Within the threatened Mojave population, gene flow among populations appears to be influenced primarily by geographic distance [24], [25], [26], [27], [28]. Indeed, the desert tortoise has been characterized as “perhaps the ideal organism for the IBD [isolation by distance] model; one that is distributed across the landscape in isolated patches and for which the difficulty of dispersal is a function of distance” [26]. This may be true at a regional scale, but makes little sense at a local scale, where desert tortoises make daily movement decisions based in part on landscape features such as elevation, slope, brush cover, and water availability. Desert tortoises are familiar with the locations of resources within their home ranges, including cover sites, mates, forage, mineral licks, and drinking sites, and can travel to these resources for several hundred meters [29]. In addition, the microhabitat within a burrow is significantly cooler and more humid than the surrounding environment [30], perhaps encouraging aggregation of tortoises around optimal burrow sites. These characteristics indicate that desert tortoises may be locally aggregated in accordance with landscape variables at a finer scale than has previously been investigated.

Local landscape effects on desert tortoise gene flow would be obscured if sampling and analyses were performed at a regional scale, where the main force driving patterns of genetic structure is probably geographic distance. Thus, fine-scale characterization of genetic structure is the most relevant for understanding the influence of the landscape on gene flow in the desert tortoise. This study had two main objectives. The main aim of this study was to assess the pattern of gene flow at a fine scale, where several potential barriers to local gene flow exist. In general, we expected to find low levels of population substructure within a local population due to the fairly high population density in the study area and high dispersal potential for desert tortoises. Our second goal was to characterize the extent to which landscape variables, such as elevation, slope, roads, habitat discontinuities, and power lines, would explain desert tortoise gene flow compared to a standard isolation-by-distance model. The integration of molecular genetic techniques with spatial analyses will provide unique insight into landscape genetics at a fine-scale, which in turn will permit more efficient conservation and management of tortoises in this region.

Results

Clustering analysis revealed that desert tortoises were structured into two genetically distinct subpopulations. The most likely number of subpopulations in the dataset was two for both geneland and structure, although K = 1 was only slightly less likely using structure (Figures 1, 2, 3). This result was not unexpected, as structure does not always detect distinct subpopulations in cases where genetic differentiation between them is weak [31]. This limitation was somewhat alleviated by the high power afforded by our microsatellite marker set (unbiased probability of identity (PID) = 1.9×10−10; [32]). There was a high level of agreement between the two clustering algorithms in terms of assignment of individual tortoises to the two subpopulations, despite the fact that the programs utilize different statistics to describe population membership (structure calculates q, the proportion of an individual's sampled genome that is characteristic of each subpopulation, and thus allows for individuals to be of mixed ancestry, whereas geneland calculates the probability that an individual belongs to each subpopulation; Figures 2 and 3).

Figure 1

Likelihood values and ΔK calculations for STRUCTURE runs.

Likelihood values obtained for 5 replicate runs of structure at K = 1 to K = 10. Mean and standard deviation across the five runs is provided. ΔK was calculated using the method outlined in [65].

Figure 2

Results of GENELAND analysis, for a) 50% threshold and b) 90% threshold.

Black circles indicate individuals assigned with ≥50% (or ≥90%) probability of membership to the ‘West’ subpopulation, and white circles indicate individuals assigned with ≥50% (or ≥90%) probability of membership to the ‘East’ subpopulation.

Figure 3

Results of STRUCTURE analysis, for a) 50% threshold and b) 90% threshold.

Black circles indicate individuals assigned with ≥50% (or ≥90%) probability of membership to the ‘West’ subpopulation, and white circles indicate individuals assigned with ≥50% (or ≥90%) probability of membership to the ‘East’ subpopulation.

The total population exhibited high levels of genetic variation. The number of alleles per locus ranged from 4 to 22, and of the 16 loci, 13 had at least 9 alleles and 9 had at least 18 alleles. The two subpopulations, denoted ‘East’ and ‘West’, exhibited weak but significant genetic differentiation (FST = 0.0046, p = 0.002), even when adjusted for the high heterozygosity of the microsatellites used (G′ST = 0.00931). This low level of genetic differentiation between subpopulations indicates either a recent separation between these two subpopulations or a high level of ongoing gene flow between them. As expected for weakly differentiated subpopulations, standard genetic diversity measures for each subpopulation were similar (Table 1). Individuals that were not assigned to either subpopulation using the more stringent assignment threshold (90%) indicate that they are of mixed or ambiguous ancestry. These individuals were scattered throughout the landscape, rather than being concentrated along a zone of contact, indicating weak genetic differentiation across the landscape (Figures 2b and 3b).

Table 1

Genetic diversity estimates for inferred subpopulations and the total population.

Population	N	A	AR	He	FIS
Total	859	14.4	13.5	0.738	0.043 (p<0.0001)
East	572	13.8	13.2	0.735	0.044 (p<0.0001)
West	287	13.6	13.5	0.738	0.041 (p<0.0001)

N, sample size; A, mean number of alleles per locus; AR, allelic richness; He, expected heterozygosity; FIS, inbreeding coefficient. P values are given in parentheses.

Genetic distance between individuals (a) was correlated with elevation (r = 0.008, p = 0.009), roads (Manix Trail r = 0.014, p = 0.00009; Fort Irwin Road r = 0.0237, p = 0.00001) and slope (r = 0.0394, p = 0.00001; Table 2). Given that elevation, roads and slope all were significantly related to genetic distance, we used partial Mantel tests to remove the influence of each variable while correlating the other with genetic distance. There was no significant relationship between genetic distance and elevation once the road- or slope-based correlation was removed. This suggests that elevation was a nuisance variable and not the main factor influencing gene flow. In contrast, the relationship between genetic distance and roads was always significant when slope was partialled out of the relationship. Likewise, the relationship between genetic distance and slope remained significant when the effect of roads was removed (Table 2). Mantel tests suggested no correlation between genetic distance and geographic distance (r = 0.004, p = 0.144), power lines (r = 0.007, p = 0.203), or soil type (r = 0.006, p = 0.013; Table 2).

Table 2

Mantel test results comparing matrices of genetic distance, topographic distance, elevation, slope, and position relative to roads and power lines among all pairs of individuals.

Correlation	Partialled out	Mantel r	P
ar 1 x Distance	—	0.0044	0.14400
ar x Elevation	—	0.0075	0.00898
ar x Manix Trail	—	0.0141	0.00009
ar x Fort Irwin Rd	—	0.0237	0.00001
ar x Power Lines	—	0.0070	0.20300
ar x Slope	—	0.0394	0.00001
ar x Soils	—	0.0056	0.01289
ar x Elevation	Manix Trail	0.0053	0.05150
ar x Elevation	Fort Irwin Rd	0.0003	0.46300
ar x Elevation	Slope	0.0050	0.05395
ar x Manix Trail	Elevation	0.0131	0.00030
ar x Manix Trail	Slope	0.0115	0.00102
ar x Fort Irwin Rd	Elevation	0.0225	0.00001
ar x Fort Irwin Rd	Slope	0.0203	0.00001
ar x Slope	Elevation	0.0390	0.00001
ar x Slope	Manix Trail	0.0386	0.00001
ar x Slope	Fort Irwin Rd	0.0375	0.00001

The measure of genetic distance used was Rousset's a inter-individual distance (Rousset 2000).

The spatial genetic autocorrelation analysis revealed small, positive genetic correlation values for only the first distance class, and values not significantly different from zero at larger distance classes. We did not see a clear shift toward negative r-values at larger distance class, but instead observed a slight decreasing trend of r-values with increasing distance classes. Although it is possible that a large effective population size could be obscuring an isolation by distance pattern [33], our results suggest an overall lack of isolation by distance in this population.

Discussion

Desert tortoises in our study area exhibit weak genetic structure at a local scale. We identified two genetically differentiated subpopulations within an apparently contiguous population. A small but significant level of genetic differentiation among the two subpopulations and similar levels of genetic diversity within subpopulations suggest that gene flow is restricted only slightly, or that the subpopulations were separated recently. Given the extremely long generation time for desert tortoises (20–25 years; [34]), ‘recent’ separation of subpopulations could have occurred dozens of years ago. Similarly, gene flow sufficient to maintain a low level of differentiation among subpopulations could be much less than one migrant per year or even one migrant every few decades in this species. This low level of population differentiation is consistent with other studies of desert tortoise in geographic areas of similar size [26].

Regardless of whether the separation between subpopulations is temporally or spatially slight, it is likely that some environmental feature (or suite of features) is restricting local movements of the desert tortoise. Barriers to gene flow that are relevant over ecological time can be difficult to detect because genetic divergence may be low [35], [36], [37]. In our study, thorough sampling (of both individuals and loci) at a local scale permitted detection of biologically relevant barriers to gene flow in spite of the fact that these barriers have resulted in only a very small departure from panmixia. These barriers warrant additional investigation to determine their continuing impact on desert tortoise movements. However, we urge caution in the interpretation of these results for conservation and management to avoid over-emphasizing the ecological significance of such small effects.

Our landscape genetic analysis identified both natural (slope) and anthropogenic (roads) landscape variables that have significantly influenced gene flow within this local population. The higher correlation of genetic distance with slope compared to roads indicates that gene flow is influenced predominantly by the slope of the landscape. Slope has been found as a barrier to gene flow in a limited number of other species, mostly pond-breeding amphibians [7], [38], [39], [40], [41], and indicates that gene flow is impeded by steep slopes. Accidental falls have been identified as a potentially important source of mortality [42] and could be serving as an effective barrier to gene flow over steep slopes. Alternatively, habitat preferences may influence desert tortoise movements and gene flow with respect to slope. In desert habitats, plant communities often differ drastically with slope due to variation in temperature, precipitation, evapotranspiration, and radiation at a micro scale. Desert tortoises are highly selective while foraging, selecting high quality plants and plant parts (those with a high Potassium Excretion Potential Index; [43]). These preferences may be even more pronounced in the western Mojave, where tortoises rely heavily on a peak of winter annual production for foraging.

Roads also influence desert tortoise movements. Desert tortoise pairs from the same side of a road exhibited significantly less genetic differentiation than tortoise pairs from opposite sides of a road. This pattern was found for both Manix Trail and Fort Irwin Road, despite differences in traffic volumes and construction materials. Both roads were constructed relatively recently (expanded to current size ∼1970 s), and it may seem that road construction is too recent to have impacted desert tortoise population structure. However, our results confirm simulation studies that suggest that new landscape features can have rapid effects on genetic structure that are detectable almost immediately [44]. The authors quantified the lag time to detect barriers in the landscape, and found a lag time of as little as one generation when using individual-based analyses, indicating that contemporary spatial genetic patterns can be caused by current landscapes [44]. This is particularly true for species such as desert tortoises, with high dispersal ability relative to the extent of the landscape and density of the population.

There are several possible ways in which roads could impact desert tortoises. First, roads may influence gene flow in desert tortoises directly via increased mortality during road crossings [45]. Second, roads may impact tortoises indirectly by providing increased human access to desert habitat, thereby magnifying anthropogenic impacts such as poaching and predation [46], [47], [48], [49]. Both of these mechanisms would serve to limit tortoise dispersal across roads, which would lead to an increase in genetic differentiation across the road as we observed. Both mechanisms also would be expected to yield a higher effect for Fort Irwin Road than for Manix Trail, a pattern consistent with our data. The third way in which roads could impact desert tortoises is by influencing burrowing behavior. Desert tortoises spend nearly 90% of their time in burrows, indicating that burrow placement is of particular importance in this species [50], [51]. Desert tortoises prefer steeply eroded banks of washes and areas with sufficient plant cover for burrow sites [52]. It is conceivable that roads provide habitat that is functionally similar to the washes and offer preferred burrow sites on either side of a road. Alternatively, plant productivity in the desert is often greater along roadsides (“edge-enhancement”; [53], [54], potentially offering preferred burrow sites [55]. The correlation between genetic distance and roads may therefore reflect an association with one or the other side of a road, and a propensity for tortoises to move along, rather than across, roads. Although these data may seem to be in direct contrast to the wealth of data documenting the negative impacts of roads on desert tortoises [29], [45], [56], other recent studies have also found desert tortoises in association with roads more often than expected [55], [57].

Our study indicates that both natural and anthropogenic features of the landscape influence patterns of genetic structure and gene flow in desert tortoises at a local scale. This is in contrast to patterns observed at a regional scale, where geographic distance plays a major role in shaping patterns of genetic structure [24], [25], [26], [27], [28]. In a recent landscape genetic study of desert tortoises, [28] identified elevation as a potential barrier to gene flow at a regional scale, and found no influence of slope (as part of a landscape resistance model). Although apparently in direct contrast to our results, it demonstrates that different landscape variables influence movements at different temporal and spatial scales [58], [59].

Although genetic structure in this population is currently affected by slope more than by roads, the extremely short time lag between the emergence of roads as a barrier and detection of a genetic effect indicates that roads may become increasingly important in shaping the evolutionary trajectory of desert tortoise populations. In addition, our results empirically support recent research suggesting that landscape effects on population structure may be limited in scale [16], but detectable almost immediately [44]. This is encouraging for those interested in characterizing the influences of recent landscape changes on fine-scale genetic structure and gene flow. The ability to detect barriers to gene flow almost immediately is also important for agencies and individuals working on endangered species management methods such as translocation, captive breeding and release, and reserve and corridor design. However, we urge caution in the interpretation of such data, as our results also demonstrate that barriers can be detected even when their effect on gene flow is relatively weak.

Methods

Ethics Statement

Permits to collect blood samples and transfer them to the laboratory included US Fish and Wildlife Service permit TE-102235-3 and California Scientific Collecting Permit SCP 801179-015.

Sample Collection

Blood samples were collected from desert tortoises in the Southern Expansion Area (SEA) of the National Training Center in Fort Irwin, CA, and within the Southern Translocation Area (STA) to the south of Fort Irwin. The SEA is an area of approximately 9,700 ha that was systematically searched for desert tortoises during a clearance operation to completely remove desert tortoises in this area for subsequent translocation to the STA. All desert tortoises found (the entire population) were sampled. In the STA, samples were collected opportunistically. Desert tortoises in this region are part of the western Mojave Recovery Unit [60]. Geographic location, sex, elevation, soil type, and slope (based on a 30-meter digital elevation model around each point, with 0 being flat and 90 vertical) were collected for each sample. The location of samples relative to three putative barriers was also noted. Barriers included two roads that run roughly north-south through the SEA and STA, and are flat relative to the surrounding landscape (i.e., not sloped or elevated): Fort Irwin Road, a paved road with relatively high traffic volume and Manix Trail, an unpaved dirt road with low traffic volume. The third barrier was a large power line corridor running southwest-northeast through the STA.

Laboratory Methods

DNA was extracted from blood samples using the BioSprint 96 robotic workstation (Qiagen) and the manufacturer's protocol. Desert tortoises were genotyped at 16 microsatellite loci developed for desert and gopher (G. polyphemus) tortoises [61], [62], [63]. All 16 loci were amplified in 4 multiplex reactions (Multiplex A: GOA 1, 11, 22, 23; Multiplex B: GOA 3, 4, 6, 14; Multiplex C: GOA2, 8, 12, 13; Multiplex D: Goag 4, Goag7, GP30, GP61). Each reaction contained ∼20 ng genomic DNA, 3–12 pmol primer, 2× PCR buffer, 2 mM each dNTP, and 1 Unit of Taq DNA polymerase. Thermocycler conditions were: 2-min initial denaturation step at 95°C; 30 cycles of 30 s at 95°C, 30 s at the annealing temperature (multiplex A = 59°C, B = 58°C, C = 56°C, D = 54°C), and 2 min at 72°C; followed by a 45 min soak at 60°C. PCR products were separated using an ABI 3130 automated DNA sequencer and the data were analyzed using genemapper version 3.7. Repeated genotyping for quality control was performed for 9% of the samples in the dataset (n = 1,291 genotypes); 7 repeated genotypes (0.5%) exhibited dropout of the larger allele and were reconciled in the final dataset. The final dataset consisted of 859 desert tortoises with a complete suite of data: multilocus microsatellite genotype, geographic location, sex, elevation, slope, and soil type.

Data Analysis

To infer the number of subpopulations in the dataset and to assign individual samples to these subpopulations, we employed two Bayesian clustering techniques, geneland (version 0.3; [64]) and structure (version 2.1; [3], [65]). In geneland, we first varied the number of populations from 1 to 6, using a matrix of genotypes, spatial coordinates for each individual, and 950,000 stored MCMC iterations (10,000,000 iterations, thinning = 10, burnin 50,000). Allele frequencies were drawn from independent Dirichlet distributions [3], as this model has been shown to perform better than the alternative model (F-model) [5]. We allowed for correlated allele frequencies, allowed for null alleles (recommended whether or not the dataset contains null alleles; [66], and an uncertainty of 10 m on spatial coordinates. This model was run 5 times to assess precision. Alternate values for uncertainty (1 m, 100 m, 1000 m), and number of stored iterations (450,000 and 1,950,000) did not significantly alter the results. Not allowing for null alleles in the model slightly reduced consistency across runs in inferring K and in individual assignment, an effect described by [66]. Individuals were assigned to subpopulations based on their probability of population membership, using thresholds for assignment of 50% (to assign all individuals in a subpopulation) and 90% (allowing samples with a probability of membership <90% in all populations to remain unassigned).

In structure, we performed 5 runs at each value of the fixed parameter K (the number of subpopulations, from K = 1 to K = 6). Each run consisted of 500,000 replicates of the MCMC after a burn-in of 100,000 replicates. We used the admixture model and allowed allele frequencies to be correlated. All other parameters were set to default values [67]. This configuration is thought to provide the best resolution in the case of potentially subtle population structure [31], [65]. We observed the common phenomenon that once the true K is reached, likelihoods for larger Ks plateau and the variance among runs increase [67]. Thus, we also used a ΔK measure that has been proposed to alleviate this problem and provide a more robust estimate of K [68]. To assign individuals to subpopulations, we performed a final run (1,000,000 burn-in and 1,000,000 replicates) at the inferred K. Values of q, the proportion of an individual's sampled genome characteristic of each subpopulation, were used to assign individuals to subpopulations at both a 50% and 90% threshold for assignment.

We assessed genetic diversity within each subpopulation (as determined by the 50% threshold geneland assignment) and the divergence among subpopulations with several standard measures; the number of alleles per locus (A), expected heterozygosity (He), and genetic variation among individuals within subpopulations (FIS) for each locus and the average multilocus value for each population. To compute FST, we used the method of [69] as implemented in spagedi (version 1.2; [70], but we used random permutations (n = 20,000) to assess significance rather than bootstrapping across loci. We also calculated a scaled measure of genetic differentiation, G′ST [71] that accounts for the high level of heterozygosity inherent in microsatellite markers.

To assess the spatial pattern of genetic variation in desert tortoises, we performed a spatial autocorrelation analysis with genalex version 6.2 [72]. Analyses were based on matrices of pairwise inter-individual genetic distances and pairwise topographic distances. We used multilocus assessments to improve replication for a pair of individuals and provide a more precise analysis of genetic pattern [73], [74], [75], [76]. An important consideration for spatial autocorrelation analysis is the choice of distance classes because this choice can influence the outcome and interpretation [77], [78]. In our analyses, we adopted distance classes that allowed relatively even sample size per distance class. We also tested other distance class options and found they revealed similar patterns. For all analyses, the pattern had levelled out within 10 km, thus we present the results up to 10 km only. We tested for significant deviation of spatial autocorrelation patterns from the random distribution of genotypes by plotting the 95% confidence intervals for the null hypothesis (estimated by 1000 random permutations of individuals genotypes among geographic locations).

We used Mantel tests to determine correlational significance between matrices of inter-individual genetic distance (Rousset's a; [79]) and several environmental factors: topographic distance, elevation, roads (including Manix Trail, which is unpaved and has low traffic volume, and Fort Irwin Road, which is paved and has high traffic volume), power lines, slope, and soil type. Pairwise distance matrices were calculated for topographic distance, elevation, and slope. For roads, power lines, and soil type, matrices were binary indicators of whether individuals were found on the same side (0) or opposite sides (1) of the putative barrier, or were sampled on the same (0) or different (1) soil types. Since only one of the two subpopulations (‘West’) contained tortoises on both sides of the power lines, we included only tortoises assigned to that subpopulation in comparisons involving the power lines. For all significant correlations, we used partial Mantel tests to determine whether correlations remained significant when the influence of other variables were removed. There has been debate regarding the use of partial Mantel tests, particularly with the use of randomizations to assess significance, which can be problematic when both variables are spatially non-random [80], [81], [82], [83]. However, this approach remains the most commonly used method for examining the association among distance matrices in population genetics (e.g., [15], [84]). In addition, recent work has shown that individual-based analyses such as Mantel tests are vastly more responsive than population-based approaches for detecting barriers to gene flow [44]. All Mantel and partial Mantel tests were performed using zt software [85]. Significance was assessed using 100,000 randomizations. We controlled for multiple tests using a false discovery rate procedure [86], [87], [88], a powerful alternative to the Bonferroni correction that seeks to minimize both type I and type II errors, with the allowed proportion of false positives set at 0.05.