Range-Wide and Regional Patterns of Population Structure and Genetic Diversity in the Gopher Tortoise

Gopher tortoises (Gopherus polyphemus) are restricted to fire-dependent, xeric longleaf pine (Pinus palustris) savannahs of the southeastern United States. In these ecosystems, wildlife biologists consider them to be keystone species as their burrows provide shelter for over 300 different species (Witz et al. 1991). Over the past century, habitat destruction and degradation (Kelly and Bechtold 1989) have led to an 80% reduction in the total number of gopher tortoise (Auffenberg and Franz 1982) with the greatest reduction found in populations west of the Tombigbee and Mobile rivers (USFWS 1987). Recent studies in this part of the range have shown population declines, poor hatching success (Noel et al. 2012), and low recruitment (Epperson and Heise 2003), even in protected areas. The dramatic decline led to the federal listing of threatened (USFWS 1987) under the US Endangered Species Act (ESA 1973, as amended) for populations west of the Tombigbee and Mobile rivers (hereafter referred to as the western or listed portion of the range). Recently, the U.S. Fish and Wildlife Service (USFWS 2011) reported that the area east of the Tombigbee and Mobile rivers also warranted listing as threatened, but because this listing was “precluded by higher priority actions” gopher tortoises were added to the candidate list. These eastern populations of gopher tortoise remain candidates for listing given that they face numerous threats, but the magnitude of these threats is currently perceived to be moderate to low (USFWS 2015).

Characterizing population genetic structure across the range of a threatened species can inform management decisions by identifying genetically distinct groups, thereby helping prevent the loss of genetic diversity and potential adaptations to the local environment (Allendorf and Luikart 2009). The criteria of genetic and ecological distinctiveness have been widely used by researchers to define units for conservation through designations such as evolutionary significant units (e.g., Ryder 1986; Waples 1991) and management units (Moritz 1994). The value of these approaches has also been codified by the federal government in an amendment to the Endangered Species Act of 1973 that allows for the recognition of distinct population segments (USFWS and NOAA 1996). Previous genetic studies of population structure in gopher tortoises have used both mitochondrial DNA (Osentoski and Lamb 1995; Clostio et al. 2012; Ennen et al. 2012) and microsatellite loci (Schwartz and Karl 2005; Clostio et al. 2012) to define genetically distinct groups and identify barriers to gene flow that would be important in a conservation and management context. In most of these studies, sampling tended to have gaps in coverage, with the focus either being on populations east (Osentoski and Lamb 1995; Schwartz and Karl 2005) or west of the Tombigbee and Mobile rivers (Clostio et al. 2012). The two most robust studies to date in terms of sampling across the range found two genetically distinct groups based on mitochondrial haplotypes (Clostio et al. 2012; Ennen et al. 2012) or five genetic groups when considering microsatellite data (Clostio et al. 2012). The authors recognized, however, that further sampling was needed. Restricted geographic sampling has also limited our capacity to detect population structure at the smaller spatial scales, with the only studies to date being the work of Schwartz and Karl (2005) in southern Georgia and Florida and Clostio et al. (2012) within the listed portion of the range.

Not only is identifying genetically distinct populations important for conservation, but understanding the patterns of genetic diversity across the landscape can provide insight into the demographic history and levels of gene flow among populations. Ennen et al. (2010) found that several populations of gopher tortoises in the listed portion of the range had significantly lower levels of genetic diversity than populations east of the Tombigbee and Mobile rivers, and Clostio et al. (2012) detected evidence of recent and historical genetic bottlenecks in the western (i.e., listed) portion of the range. While reduced genetic variation could be the product of population bottlenecks, Ennen et al. (2010) pointed out that historically this area may have possessed lower levels of genetic diversity. Presumably both effective population sizes and levels of gene flow are higher at the center of a species range, acting as the genetic drivers behind lower levels of genetic diversity and higher levels of genetic differentiation in the periphery (Eckert et al. 2008; Dai and Fu 2011).

While a number of studies have looked at the population genetics of gopher tortoises, a thorough assessment of population structure is lacking at both the range-wide level as well as finer spatial scales. Herein, we attempt to comprehensively address both spatial scales by doubling the data available for the number of sampling sites, individuals, and microsatellite loci compared to previously published work. The first objective of this study was to identify population structure across the range using a variety of analytical approaches and then use our more intensive geographic sampling to characterize potential population structure within each of the genetically defined groups. Secondly, we wanted to compare patterns of genetic diversity, gene flow, and demographic history across the range. The ultimate goal is to use the insights obtained from this work to guide management decisions for the species.

Blood samples were collected between 2000 and 2012 by the authors or obtained from colleagues for 933 individuals collected from 47 sampling sites across the range (Table 1; Figure 1) following the methods outlined in Ennen et al. (2012). This included 14 sampling sites reported on by Schwartz and Karl (2005). Given the history of gopher tortoise translocation, especially in Florida (Mushinsky et al. 2006), we made every effort to avoid collecting in areas known to have received relocated individuals. However, in two cases involving translocations we either used individuals known to be prior residents (site 33) or we used the donor site in the analysis (site 30). We stored each blood sample in a 2.0-mL vial with approximately 0.5 mL of SED tissue preservation buffer (Seutin et al. 1991). We extracted genomic DNA using a DNeasy Tissue Kit (Qiagen) and genotyped each individual for 20 of the loci (Table A1, Archived Material in Dryad, doi:10.5061/dryad.nk064) described by Kreiser et al. (2013). We performed polymerase chain reactions on a Veriti 96-well Thermal Cycler (Applied Biosystems) in 12.5-μL reactions consisting of 50 mM KCl, 10 mM Tris-HCl (pH 8.3), 0.01% gelatin, 2.0 mM MgCl₂, 200 μM dNTPs, 0.1875 units of Taq polymerase (New England Biolabs), 0.3 μM of the M13 tailed forward primer (Boutin-Ganache et al. 2001), 0.3 μM of the reverse primer, 0.1 μM of the M13 labeled primer (LI-COR), 20–100 ng of template DNA, and water to the final volume. Polymerase chain reaction cycling conditions consisted of an initial denaturing step of 94°C for 2 min. We followed this with 35 cycles of 30 sec at 94°C, 1 min at an annealing temperature of 50–60°C depending on the locus, and ended with 1 min at 72°C. We completed the cycle with a final elongation step of 10 min at 72°C. We visualized microsatellite alleles using a LI-COR 4300 DNA Analyzer along with a 50–350-bp size standard (LI-COR) and scored them using Gene Image IR v. 3.55 (LI-COR).

At each sampling site we tested loci for Hardy–Weinburg equilibrium and linkage disequilibrium using GENEPOP on the web (Raymond and Rousset 1995) with the alpha level of these tests adjusted by a sequential Bonferroni correction (Rice 1989). We also tested loci for the presence of null alleles using MICROCHECKER (Van Oosterhout et al. 2004). We calculated basic summary statistics (allelic richness: A, observed heterozygosity: H_o, and expected heterozygosity: H_e) with FSTAT 2.9.3 (Goudet 2001). We excluded sampling sites 4 and 10 (Theodore Mars Wildlife Management Area and Lynn's House) from the calculation of allelic richness due to low sample sizes. We also tested for significant differences in observed heterozygosity (H_o), expected heterozygosity (H_e), and allelic richness (A) among the genetic groups defined by our Bayesian analysis of population structure (see below). We tested all data for normality using a Shapiro–Wilk test, and conducted the appropriate parametric (analysis of variance with Tukey HSD) or nonparametric (Kruskal–Wallis with a nonparametric comparison for each pair using a Wilcoxon method) using JMP 7.0.1 (SAS Institute Inc. 2007) with an alpha level of 0.05.

We first characterized population structure across the range of gopher tortoises by employing an analysis of molecular variance (AMOVA; Excoffier et al. 1992) as implemented in ARLEQUIN v. 3.5 (Excoffier and Lischer 2010) with significance assessed by 1,000 random permutations. We grouped sampling sites based on putative barriers to dispersal (i.e., rivers) and ecological gradients (i.e., physiographic provinces; Table 2). The first model (“Tombigbee–Mobile”) assigned sampling sites into two groups as delineated by the Tombigbee and Mobile rivers, which are used to define the listed portion of the range (USFWS 1987). The second model (“Apalachicola–Chattahoochee”) also tested two groups separated by the phylogeographic break at the Apalachicola River previously identified by an analysis of gopher tortoise mitochondrial data (Ennen et al. 2012). We also assigned sampling sites into three groups based on divisions of the Coastal Plain physiographic province (model 3, “Physiographic Provinces”) or the groups were prescribed by the Apalachicola and Chattahoochee rivers and the Suwannee and Okefenokee–St. Mary's rivers (model 4, “Apalachicola-Suwannee”). The fifth model (“Mobile–Apalachicola–Suwannee”) grouped sampling sites on basis of all of the major river systems described in models 1, 2, and 4, while the sixth model (“Mobile–Apalachicola–Physiographic”) divided sampling sites on the basis of both rivers and physiographic regions. The last model (“Mobile–Apalachicola–Physiographic–Ridge”) was built on model 5 through the inclusion of two groups within peninsular Florida defined as east and west of the central ridge system.

The program STRUCTURE uses a Bayesian approach to assign individuals to some number of groups that are in Hardy–Weinberg and linkage equilibrium. We used STRUCTURE v. 2.3.4 (Pritchard et al. 2000) to test values of K from 1 to 10 using a model of admixed ancestry and assuming correlated allele frequencies between groups without using population information as a prior. For each value of K, we ran 20 replicates with 500,000 Markov chain Monte Carlo iterations and a burn-in of 200,000. We selected the best value of K by examining the likelihood values for each K and by using the ΔK method (Evanno et al. 2005) as calculated by the program Structure Harvester v. 6.92 (Earl and von Holdt 2012). We then averaged the 20 runs for the best value of K with CLUMPP v. 1.1.2 (Jakobsson and Rosenberg 2007) and visualized the average ancestry coefficient (q) for each individual with DISTRUCT v. 1.1 (Rosenberg 2004). We selected the best value of K as the smallest value that seemed to represent biologically meaningful genetic groups. Hereafter, we also use the term region to refer to the areas delineated by the sampling sites comprising these genetic groups. Since intraregional population structure was still possible, we performed another set of analyses on individuals within each of the genetic groups (e.g., hierarchical analysis; Vähä et al. 2007). We used the same parameters for the STRUCTURE analyses with 20 iterations of each value of K from 1 to 6.

We used FSTAT 2.9.3 to calculate the unbiased estimator of F_ST (θ; Weir and Cockerham 1984) among sampling sites. We also calculated the mean pairwise F_ST values for sampling sites within and among the genetic groups as defined by the STRUCTURE analysis. Lastly, we pooled individuals within each of the genetic groups and calculated pairwise F_ST values.

We employed two methods for estimating gene flow among the regions defined by the Bayesian analyses. The first was the coalescent-based approach of Migrate-n that examines historic levels of gene flow on the order of approximately 4N_e (effective population size) generations (Beerli 2009). The second was the disequilibrium-based method employed by BAYESASS that measures gene flow over the past several generations (Wilson and Rannala 2003). Both methods have been widely used in studies with a conservation focus, some of which contrast the two measures to make inferences about the impact of recent changes in the landscape on the organism (Epps and Keyghobadi 2015). However, recent simulation work (Samarasin et al. 2016) suggests that both the magnitude and timing in changes of connectivity between populations can considerably influence the accuracy of the migration rates obtained through these analyses. Given these potential biases, we interpret our results with caution.

We estimated contemporary migration rates (m_ij), defined as the proportion of individuals in population i that arrived from population j, with BAYESASS v. 3.0.4 (Wilson and Rannala 2003). Preliminary runs determined that the default conditions for the mixing parameters resulted in acceptance rates between 0.2 and 0.6 for the allele frequency, migration rate, and inbreeding coefficient as per the authors' recommendations. We performed three independent runs using different initial seed values for a total of 10,000,000 iterations with a burn-in of 1,000,000. After removing the burn-in, we collected samples every 100 generations, and used these to infer the posterior probability distributions of migration rates among the genetic groups. To ensure our analyses were run long enough for convergence, we examined the trace files using Tracer v. 1.5 (Rambaut and Drummond 2007) and we compared parameter estimates across the three independent runs. We considered migration rates significantly different at alpha of 0.05 if the 95% credible sets did not overlap between the genetic groups.

We estimated historical migration rates among the genetic groups using Migrate-n (version 3.6.4; Beerli and Felsenstein 2001; Beerli 2006) where movement was constrained to occur between geographically adjacent areas. We used the Brownian mutation model, a starting genealogy derived from a Unweighted Pair Group Method with Arithmetic Mean tree and initial θ and M values derived from the F_ST calculation. Additionally, the analysis used uniform priors for θ with minimum, maximum, and delta values of 0.01, 100.0, and 9.99, respectively. We ran four independent chains with static heating using temperature settings of 1.0, 1.5, 3.0 and 1,000,000.0. The Markov chain Monte Carlo search included a total of 5,000,000 steps, recorded every 100 generations, with the first 10,000 discarded as burn-in. We assessed convergence by determining if there was an effective sample size > 1,000 and a unimodal posterior distribution for each parameter.

We checked individual sampling sites and Bayesian-defined regions for evidence of recent genetic bottlenecks using the program BOTTLENECK (Piry et al. 1999) and the M ratio test (Garza and Williamson 2001). Peery et al. (2012) cautioned that the mutation model parameters used in these tests can influence the ability to detect genetic bottlenecks. Hence, for each analysis we used a range of parameters including the ones suggested by the authors of the method as well as values identified as appropriate by Peery et al. (2012) from their simulations. For the BOTTLENECK analysis, we first removed all sampling sites with fewer than 10 individuals (Table 1). We ran the two-phased model of mutation test for 1,000 iterations. Variance values included 10, 12, and 30% while we set the proportion of stepwise mutations to 70, 78, 80, or 90%. We determined significance of heterozygosity excess using the two-tailed Wilcoxon sign rank test and a sequential Bonferroni correction (Rice 1989). For the M ratio test we first used some of the values suggested by the authors with proportions of one-step mutations of 70 and 90%, an average size of non–one-step mutations (Δg) of 3.5, and a value of 10 for θ. We then ran the analysis with the proportion of one-step mutations set to 78%, an average size of non–one-step mutations (Δg) of 3.1, and θ values of 5 and 10. We compared the calculated M ratios to a critical value for M for each population that was obtained from the 95% threshold of 10,000 simulations of an equilibrium population.

We genotyped 933 individuals from 47 sampling sites located across the range (Table 1; Table A1, Archived Material in Dryad, doi:10.5061/dryad.nk064; Table A2, Archived Material in Dryad, doi:10.5061/dryad.nk064). None of the 20 loci consistently deviated from Hardy–Weinburg equilibrium or demonstrated linkage disequilibrium after a sequential Bonferroni correction nor did we detect null alleles for any locus. Average observed and expected heterozygosity values across collections ranged from 0.388 to 0.743 and from 0.415 to 0.703, respectively, while average allelic richness ranged from 2.23 to 4.33 (Table 1). We made multiple attempts to obtain complete genotypes for each individual if necessary, and ultimately 855 individuals (92%) possessed 10% or less missing data.

Each AMOVA model explained a significant portion of the molecular variance among groups of sampling sites. However, the Mobile–Apalachicola–Physiographic model (model 6; five groups) and the Mobile–Apalachicola–Physiographic–Ridges model (model 7; six groups) had the highest among and lowest within-group variances (18.0/7.0% and 18.1/6.5%, respectively) compared to the other models (Table 2). Admittedly the amount of among-group variance partitioned by the models did not differ greatly. For example, the among-group variance for the models with two groups of sampling sites (Tombigbee–Mobile and Apalachicola–Chattahoochee) was similar (15.7 and 14.6%, respectively) to the variance explained by models with three groups of sampling sites (Physiographic Provinces = 14.5%; Apalachicola–Suwannee = 15.9%).

The ΔK analysis of the STRUCTURE results recognized the deepest level of population subdivision at K = 2, corresponding to the split in the range as defined by the Tombigbee and Mobile rivers. However, the likelihood values began to plateau around K = 5. We have designated regions (Western, Central, West Georgia, East Georgia, and Florida) across the range reflecting the distribution of these genetic groups (Table 1; Figure 3). Twenty-six of the 47 sampling sites possessed an average ancestry coefficient (q) score of > 0.9 (Table 3), although there were multiple individuals with q scores > 0.9 at 43 of them. For collections located at the boundary between regional groups it was not uncommon to find individuals with mixed ancestry (Figure 2) such as between the Western and Central (13–16), Central and West Georgia (20 and 21), West Georgia and East Georgia (28), and East Georgia and Florida regions (35). We detected at least three cases of sampling sites with an individual with a q score > 0.9 in a different genetic group, putatively the result of being transported there by humans (i.e., waifs). This included one individual from a Western sampling site with ancestry in the Florida–East Georgia group, one in an East Georgia sampling site with ancestry in the Florida group and one in a West Georgia sampling site with ancestry in the East Georgia group. At two West Georgia sampling sites we found four individuals with a high degree of ancestry in the Central group (q > 0.9). However, since both of these locations (20 and 21) are on the border between the Central and West Georgia regions they may not necessarily be the result of human-mediated movement.

Varying degrees of regional population structure was evident within each of the five regions identified by the range-wide analysis (Figure 4). Within the Western region, the analysis detected two genetic groups. The first genetic group was represented by sampling sites (1–8, 11, and 12) where all individuals possessed q scores > 0.9. The remaining sampling sites located east of the Pascagoula River or just to the west of the river (9, 10, 13, and 14) had ancestry in the second genetic group, although individual q scores were much lower. The Central region likewise possessed two genetic groups that fell out along a west–east split with sampling sites 15 and 16 in one group and sites 17–19 in the other, and most individuals demonstrated a high degree of ancestry (q > 0.9) in their respective group. We also detected two genetic groups within the West Georgia region represented by southwestern sampling sites (20 and 21; 23–25) in one while the eastern sampling sites (26 and 27) and extreme northwestern sampling site (22) were in the other. At the eastern periphery of the range, we detected four genetic groups representing collections in the western (28), northern (33), eastern (31 and 32), and southern (34) parts of the East Georgia region. The remaining two sampling sites showed admixture between either the western–eastern groups (29) or eastern–southern groups (30). Within the Florida region, we also detected four genetic groups. The first group was represented by the northern locale (35), the second group contained three sampling sites along the eastern coast of Florida (38, 40, 45), the third group was comprised of collections in the central and western coast of Florida (36, 39, 41–43), and the fourth group was represented by the southwestern island (47). The remaining sampling sites (37, 44, and 46) mostly reflected varying degrees of admixture among geographically adjacent groups.

Pairwise F_ST values among sites ranged from 0.014 to 0.309 (Table A3). When we pooled these F_ST values according to the five regions identified by the STRUCTURE analysis, mean pairwise F_ST values among regions ranged from 0.098 to 0.223 (Table 4). For comparisons among sampling sites within a region, the Western region had the lowest mean F_ST (0.042) while the East Georgia region had the highest value (0.086). When sampling sites were pooled at the regional level, the lowest mean F_ST value was between the Central and West Georgia regions (0.085) with the highest between the Western and East Georgia regions (0.288).

The contemporary migration rates estimated by the three independent runs of BAYESASS were concordant so we only report one of them. In general, migration rates between regions were low with the highest proportion (> 0.9) of individuals derived from the source region (Table 5). We found the highest migration rates among adjacent regions, but these rates were not always symmetric. For example, migration was greater from the Central to Western regions and the East Georgia to West Georgia regions than it was in the opposite directions. In both cases the 95% credible sets overlapped, although there was only a marginal degree of overlap in the Western–Central comparison. Migration rates between the Western Georgia and Florida regions were small in both directions with both containing zero within the 95% credible set. We appear to have achieved convergence in the Migrate-n analysis as parameters possessed effective sample size scores over 1,000 and were represented by a unimodal distribution. Modal mutation-scaled migration rates were variable across the range and tended to be asymmetric between regions (Figure 5). Migration rates between Florida and East Georgia were the highest. Between other regions the modal values were lower, although it should be noted that the 95% highest probability density in these cases included zero.

In general, sampling sites from the extreme western and eastern parts of the range as well as the southern Florida collection (site 47, Cayo Costa) had lower levels of genetic diversity. We averaged genetic diversity values across sampling sites grouped into the five genetically defined regions (see above). Both H_o and H_e differed among regions using a Kruskal–Wallis test (χ² = 31.97 and 22.32, df = 4 for both, P = 0.0002 and < 0.0001, respectively; see Table 1 for post-hoc comparisons). Also, allelic richness (A) differed among regions using an analysis of variance (F_4,44 = 27.20, P < 0.001, see Table 1 for post-hoc comparisons). After applying Bonferroni corrections, our two bottleneck analyses (i.e., BOTTLENECK and M ratio) found no evidence of bottlenecks at the regional level (K = 5). One sampling site (33) possessed significant bottleneck test results for all analyses while two others (1 and 4) were also significant for runs that did not include a 90% proportion of stepwise mutations (Table 1). Three different sampling sites (24, 29, and 35; Table 1) showed evidence of a bottleneck in the M ratio test for runs that included the following parameters: 1) 90% proportion of one-step mutations, Δg of 3.5, θ = 10 or 2) 78% proportion of one-step mutations, Δg of 3.1, θ = 5.

Similar to previous work (Clostio et al. 2012), our range-wide Bayesian analyses clustered individuals into five genetically distinct regions (i.e., Western, Central, West Georgia, East Georgia, and Florida) throughout the distribution of the gopher tortoise. Likewise, the results of the AMOVA analysis seemed to suggest that grouping sampling sites into five or six regions best partitioned the genetic variation among regions. Within these regions, most populations were comprised of individuals with > 90% ancestry in a particular genetic group, although admixture did tend to occur within collections at the boundaries of adjacent genetic groups. Major river systems and physiographic provinces formed the boundaries between the genetically defined regions. The two major riverine boundaries were the Tombigbee–Mobile and Apalachicola–Chattahoochee, which separated three of the regions (i.e., Western, Central, and West Georgia).

The Tombigbee and Mobile rivers defined the Western group and this river system is also used (USFWS 1987) to define the listed portion of the range. The Western group was distinct with only four sampling sites where the average q score was < 90%, two of which (13 and 14) were located along this boundary. The second barrier, as previously reported (Clostio et al. 2012; Ennen et al. 2012), was formed by the Apalachicola and Chattahoochee rivers. Sampling sites along the border between the Central and West Georgia regions (sites 20 and 21) demonstrated a higher degree of admixture. These two sampling sites were likely influenced by the construction of the reservoir (see Ennen et al. 2011) and may possess individuals from both sides of the river. We had limited samples along the Suwannee River; therefore, we could not fully assess the importance of this river as a barrier to gene flow. However, one of the groups identified in Schwartz and Karl's (2005) AMOVA was comprised of collections on both sides of the Suwannee River suggesting that this river does not represent a long-term barrier to gene flow like the other rivers (i.e., Tombigbee–Mobile and Apalachicola–Chattahoochee river systems).

The remaining delineations among groups appeared to be associated with transition areas between the East Gulf Coastal Plain, Sea Island, and Floridian physiographic provinces (Figure 2). Similar to the regional boundaries defined by rivers, individuals in populations at the boundaries of the physiographic provinces demonstrated admixture with the adjacent genetic groups. At the University of Northern Florida collection (site 35), located along the transition between the Sea Island and Floridian provinces, all but three individuals showed mixed ancestry (q < 0.9) between the East Georgia and Florida groups. Likewise, along the East Gulf Coastal Plain and Sea Island transition, at sampling site 28 in the East Georgia region, 20 of 28 individuals demonstrated varying degrees of admixture with the Western Georgia group. We are cautious about interpreting the degree to which these provincial boundaries serve as barriers to gene flow for two reasons. First, this association could be due to our limited sampling along these transition areas. Second, these physiographic provinces do not possess well-defined boundaries (Thornbury 1967), thus explaining some of the admixture at locations along their transitions.

We recovered some degree of additional genetic structuring (i.e., intraregional groups) within each of the five regional genetic groups: K = 2 in the Western, Central, and Western Georgia regions, but K = 4 in the East Georgia and Florida regions. The weakest level of intraregional population structure was in the west where, like Clostio et al. (2012), we recovered a genetic break that roughly coincided with the Pascagoula and Chickasawhay rivers. Similarly, the Pascagoula River has been suggested as a contact zone between distinct groups of the ground skink (Scincella lateralis) from opposite sides of the river (Jackson and Austin 2010). However, the intraregional group east of the river was much more weakly defined in our study, with high levels of admixture in individuals from these sampling sites (q usually between 0.40 and 0.60) compared to q scores > 0.90 reported by Clostio et al. (2012). Perhaps this discrepancy reflects differences in the number of loci genotyped and the geographic extent of sampling. We do believe that the intraregional differentiation detected within the Western genetic group is real. In the range-wide analysis, all individuals within the Western group had high assignment values (q > 0.9), and therefore the intraregional genetic group we detected east of the Pascagoula and Chickasawhay rivers does not appear to be the result of the presence of individuals with admixture from the Central group.

In the Central and Western Georgia genetic groups, our results support a K = 2 for both regions. In the Central region, the two intraregional groups appear to be identified by a well-defined split between populations east (sites 17–19) and west (sites 15 and 16) of the Escambia River. However, with only five collections represented within this region we are cautious about interpreting the importance of the Escambia River, and sampling more locations in closer proximity to the eastern side of the river is necessary to clarify its impact on gene flow. Within the Western Georgia region, sampling sites fell into southwestern (20, 21, 23–25) and northeastern (22, 26, 27) intraregional groups, but the break between the two did not appear to correspond to any particular riverine or physiographic boundary. While the division of individuals into the two intraregional groups is supported by relatively high ancestry coefficients, there is still admixture evident at most locales, suggesting contemporary gene flow between these intraregional genetic groups.

With a K = 4, both the East Georgia and peninsular Florida genetic groups showed greater intraregional genetic differentiation than the previous three regions. In the East Georgia region, not all sampling sites cleanly fell into one of these intraregional four groups, but varying degrees of admixture demonstrated that gene flow does connect the regions. Notably, the Aiken Gopher Tortoise Preserve sampling site (33) formed its own intraregional genetic group. This collection is at the northern extent of the range for gopher tortoises and is isolated by > 120 km from other native populations within South Carolina. The resident population only comprised approximately 12 adults (Clark et al. 2001; T. Tuberville, unpublished data), but despite the small population size it has the highest level of heterozygosity of any sampling site. Furthermore, in the range-wide STRUCTURE analysis these individuals show admixture of both Florida and East Georgia groups. Although Holbrook (1842) noted that the species was historically common in the region and the population was presumed to be native (Clark et al. 2001), historical releases cannot be ruled out. In peninsular Florida, sampling sites with strong ancestry in any one of the four intraregional groups tended to be found at coastal locales either along the Atlantic (35, 38, 40, and 45) or Gulf of Mexico (41–43, 47). The remaining sampling sites reflected varying degrees of admixture among adjacent groups. These four intraregional genetic groups roughly correspond to those described in Schwartz and Karl (2005). Our two northeastern collections (35 and 37) are in a similar geographic area as those that made up a north Florida intraregional group in their study. Their central–west coast group is very similar to ours, except ours includes a more northern sampling site (36) in Alachua County. Schwartz and Karl (2005) grouped Cayo Costa (47) and Highland Hammock (44) and our analysis agreed that Highland Hammock does show some ancestry from the Cayo Costa genetic intraregional group. While Schwartz and Karl (2005) only had one collection representing the southern Atlantic coast (Jonathan Dickenson, site 45) our additional sampling indicated that this intraregional group extends farther north along the coast to include sampling sites 38 and 40. This seems to correspond to the Atlantic Coast Ridge, which has influenced population genetic structure in other upland vertebrates in Florida (Clark et al. 1999; Branch et al. 2003).

We found significantly lower genetic diversity at the periphery of the range of the gopher tortoise, a pattern that has also been reported in other Gopherus species (Fujii and Forstner 2010; Hagerty and Tracy 2010). Specifically, like Ennen et al. (2010) and Clostio et al. (2012), we found reduced genetic diversity in the Western region, but we also found significantly lower levels in the East Georgia region as well. One possible explanation for the reduction in genetic diversity is that populations in the Western and East Georgia regions have experienced extensive population bottlenecks due to anthropogenic disturbances. Our bottleneck tests showed no evidence of an overall signature of genetic bottlenecks for sampling sites within these two regions. However, we interpret this result with caution as simulation studies indicate that bottleneck tests can be sensitive to a number of factors, including low sample sizes, and empirically a number of studies have failed to detect signatures of genetic bottlenecks even in cases where there were documented population declines (Peery et al. 2012).

Contemporary levels of gene flow were asymmetric from central to peripheral regions with higher levels of gene flow from the Central to Western regions. However, gene flow was greater between the East Georgia and West Georgia regions, although the 95% credible set did overlap. Historic rates of gene flow tended to show different patterns, although their interpretation is complicated by the broadly overlapping 95% highest probability densities. The higher level of genetic diversity in the Central region might be a consequence of historic gene flow from adjacent groups (Western and West Georgia). This would agree with Ennen et al.'s (2012) suggestion that two mitochondrial clades diverged in isolation, but later the two groups expanded across their current range, producing a large contact zone in Alabama, Western Georgia, and the panhandle of Florida. Two other observations of note are that gene flow (historic and contemporary) has been low between the Florida and Western Georgia groups (i.e., through the Florida panhandle), while gene flow at both time scales has been the highest between the East Georgia and Florida regions.

Our findings have important conservation and management implications. Analysis of our more extensive geographic sampling supported Clostio et al. (2012) in the recognition of five genetically distinct groups across the distribution of the gopher tortoise. The regions delineated by these five groups would seem to warrant recognition as management units in terms of conservation planning. In addition, we found additional population structure within each of these five regions that point to the potential significance of rivers (e.g., the Pascagoula and Escambia) and physiographic features (e.g., the Atlantic Coast Ridge in Florida) as barriers to gene flow. However, additional work in evaluating this finer-scale structure is needed. Wildlife managers should consider the presence of regional and intraregional differentiation among populations as part of management plans that involve translocations. Wildlife managers have relocated gopher tortoises, especially those in Florida, more extensively than any other amphibian or reptile (Dodd and Seigel 1991). Researchers have devoted considerable effort to evaluating the methodology and success of head-starting efforts (rearing hatchlings until they reach a larger size) and translocations (Heise and Epperson 2005; Tuberville et al. 2005, 2008, 2015) in augmenting or reestablishing populations, and wildlife managers have recognized both as potential conservation strategies for gopher tortoises (USFWS 2012).

Please note: The Journal of Fish and Wildlife Management is not responsible for the content or functionality of any archived material. Queries should be directed to the corresponding author for the article.

To cite this archived material, please cite both the journal article (formatting found in the Abstract section of this article) and the following recommended format for the archived material.

Gaillard D, Ennen JR, Kreiser BR, Qualls CP, Sweat SC, Birkhead R, Tuberville TD, Aresco M, McCoy ED, Mushinsky HR, Hentges TW. 2017. Data from: Range-wide and regional patterns of population structure and genetic diversity in the gopher tortoise. Journal of Fish and Wildlife Management, 8(2):xxx–xxx. Archived in Dryad Digital Repository: doi:10.5061/dryad.nk064

Table A1. Genotype data for 20 microsatellite loci for each of the 933 gopher tortoise (Gopherus polyphemus) individuals collected between 2000 and 2012 from 47 sampling sites located across the southeastern United States. Found at DOI: http://dx.doi.org/10.5061/dryad.nk064 (221 KB XLSX).

Table A2. Allele frequencies by sampling site for each of the 20 microsatellite loci genotyped in this study for the 933 gopher tortoise (Gopherus polyphemus) individuals collected between 2000 and 2012 from 47 sampling sites located across the southeastern United States.

Found at DOI: http://dx.doi.org/10.5061/dryad.nk064 (112 KB XLSX).

Table A3. Pairwise genetic and geographic distances for each of the 47 gopher tortoise (Gopherus polyphemus) sampling sites across the southeastern United States collected between 2000 and 2012 with F_ST values listed below the diagonal and distance in kilometers above the diagonal. Found at DOI: http://dx.doi.org/10.5061/dryad.nk064 (69 KB XLSX).

Reference A1. Kelly JF, Bechtold WA. 1989. The longleaf pine resource. Pages 11–22 in Proceedings of the symposium on the management of longleaf pine. USDA Forest Service, Southern Forest Experiment Station, Gen. Tech. Report SO-75.

Found at DOI: http://dx.doi.org/10.5061/dryad.nk064 (439 KB PDF).

Reference A2. [USFWS] U.S. Fish and Wildlife Service. 2012. Range-wide conservation strategy for the gopher tortoise. Jackson, Mississippi.

Found at DOI: http://dx.doi.org/10.5061/dryad.nk064(827 KB PDF).

We would like to thank the USFWS (S. Ginger and M. Hinderliter), the Gopher Tortoise Council, Mississippi Department of Wildlife, Fisheries, and Parks (R. Jones), and U.S. Army Corps of Engineers- Engineer Research and Development Center-Construction Engineering Research Laboratory (T. Smith) for providing funding and support for this project. Partial support for sample collection and manuscript preparation by TDT was also provided by the Department of Energy (DE-FC09-07R22506) to the University of Georgia Research Foundation. We thank N. Sharp, D. Hartley, L. McCoy, K. Amatuli, J. Goynyer, C. Walters, K. Shelton, and R. Bolt for support with tissue collection. We thank J. Schaefer for assistance with Migrate-n runs on the cluster maintained by the School of Computing at the University of Southern Mississippi. The Associate Editor and two anonymous reviewers provided useful comments that improved the earlier version of this manuscript.

Any use of trade, product, website, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. Government.