Post‐Pleistocene differentiation in a Central Interior Highlands endemic salamander

Abstract Aim For many endemic species with limited dispersal capacities, the relationship between landscape changes and species distributions is still unclear. We characterized the population structure of the endemic ringed salamander (Ambystoma annulatum) across its distribution in the Central Interior Highlands (CIH) of North America, an area of high species endemism, to infer the ecological and evolutionary history of the species. Methods We sampled 498 individuals across the species distribution and characterized the population genetic structure using nuclear microsatellite and mitochondrial DNA (mtDNA) markers. Results Ambystoma annulatum exist in two strongly supported nuclear genetic clusters across the CIH that correspond to a northern cluster that includes the Missouri Ozark populations and a southern cluster that includes the Arkansas and Oklahoma Ozarks and the Ouachita Mountains. Our demographic models estimated that these populations diverged approximately 2,700 years ago. Pairwise estimates of genetic differentiation at microsatellite and mtDNA markers indicated limited contemporary gene flow and suggest that genetic differentiation was primarily influenced by changes in the post‐Pleistocene landscape of the CIH. Main Conclusions Both the geologic history and post‐European settlement history of the CIH have influenced the population genetic structure of A. annulatum. The low mtDNA diversity suggests a retraction into and expansion out of refugial areas in the south‐central Ozarks, during temperature fluctuations of the Pleistocene and Holocene epochs. Similarly, the estimated divergence time for the two nuclear clusters corresponds to changes in the post‐Pleistocene landscape. More recently, decreased A. annulatum gene flow may be a result of increased habitat fragmentation and alteration post‐European settlement.


| INTRODUC TI ON
Species distributions are governed by factors that operate across ecological and evolutionary time frames, such as changes in temperature and precipitation patterns, shifts in species interactions, and demographic stochasticity (Brown, Stevens, & Kaufman, 1996;Hewitt, 1996;Roy, Valentine, Jablonski, & Kidwell, 1996). Integrating information about the structure and dynamics of species' contemporary geographic ranges with the geologic history of a landscape allows for a better understanding of the ecology and evolutionary history of a species. Processes acting over evolutionary time frames (e.g., shifts in climate associated with glaciation events, speciation events associated with changing landscape dynamics during geological formations, shifts in paleodrainages) have lasting effects on the structure of species ranges (Hewitt, 1996). For instance, contraction of species distributions into refugial areas during the glaciation cycles of the Pleistocene (2.6-0.012 Ma) led to high rates of genetic differentiation within species (Hewitt, 2000;Lee-Yaw, Irwin, & Green, 2008;Puckett, Etter, Johnson, & Eggert, 2015;Zamudio & Savage, 2003) and in some cases allopatric speciation (Avise & Walker, 1998;Knowles, 2000;Maddison & McMahon, 2000;Smith & Farrell, 2005).
The Central Interior Highlands (CIH) ecoregion of North America, which includes the Ozark Highlands and Ouachita Mountains of Missouri, Arkansas, and Oklahoma, has experienced shifts in species assemblages over ecological and evolutionary time frames, resulting in high levels of endemism (Soltis, Morris, McLachlan, Manos, & Soltis, 2006). This ecoregion was once part of a larger highland region that included the Appalachian Mountains to the east (Mayden, 1985). Geologic evidence suggests that separation of the highlands into Eastern and Central regions occurred during the Cretaceous to mid-Miocene (65-12 Ma; Cushing, Boswell, & Hosman, 1964;Galloway, Whiteaker, & Ganey-Curry, 2011). Although the Ozarks and Ouachitas date to the Ouachita orogeny of the Pennsylvanian epoch (318-299 Ma), they were formed by distinct geological processes and have been disjunct since their formation (Thornbury, 1965). The Ozarks are a northeast to southwest oriented limestone uplift with large amounts of karst habitat whereas the Ouachitas are an east-west-oriented range of folded mountains (Thornbury, 1965).
During the Pleistocene, the CIH were further bisected by the development of the Arkansas River system between the Ozarks and Ouachitas (Mayden, 1985;Wiley & Mayden, 1985). As a result of its complex geologic history, the CIH have served as a refugium for species during glacial cycles as well as the center of origin for many other species, notably fishes (Mayden, 1985;Near & Keck, 2005). During the Pleistocene glacial cycles, temperate habitats shifted south to the Gulf Coast and the CIH was covered with boreal flora (Soltis et al., 2006). At present, this ecoregion is composed of temperate broadleaf forest and mixed broadleaf forest (Flader, 2004) that contains more than 200 endemic species, of which over 160 are Ozark endemics and 48 are Ouachita endemics (Buhay & Crandall, 2005;Crandall, 1998;Mayden, 1988;Ouachita Ecoregional Assessment Team, 2003;Ozarks Ecoregional Assessment Team, 2003).
Phylogeographic studies of CIH species have found little concordance in the observed patterns of intraspecific genetic structure (Hardy, Grady, & Routman, 2002). For some taxa, the Arkansas River Valley is a barrier to movement, creating distinct Ozark Highland and Ouachita Mountain groups (Grady, Cashner, & Rogers, 1990;Herman & Bouzat, 2016;Puckett et al., 2014). In plethodontid salamanders (Herman & Bouzat, 2016;Martin, Shepard, Steffen, Phillips, & Bonett, 2016) and freshwater fishes (Grady et al., 1990;Hardy et al., 2002), genetic differentiation was found to be lower between the Eastern Highlands and the CIH than observed between the Ozarks and Ouachitas; thus, the two CIH regions were likely colonized independently and have experienced low rates of gene flow. However, other taxa, especially aquatic taxa, cluster more strongly based upon watersheds within the CIH than geologic formation, creating phylogeographic breaks between the northern Ozark Highlands of Missouri, where rivers drain into the Missouri and Mississippi Rivers, and southern populations in the Arkansas and Oklahoma Ozarks, where rivers drain into the Arkansas River (Hardy et al., 2002;Hutchison & Templeton, 1999;Routman, Wu, & Templeton, 1994).
Although aquatic and terrestrial species exhibit similar patterns of gene flow across broad scales that largely correspond to major drainages, there is less concordance at finer scales (i.e., subdrainage basins within a watershed). For aquatic organisms, gene flow is higher within subdrainage basins than among subdrainage basins as movement is constrained within aquatic ecosystems (Crowhurst et al., 2011;Hardy et al., 2002). Movements of terrestrial species, however, are less constrained, and greater genetic differentiation has been observed along geological formations (Herman & Bouzat, 2016;Martin et al., 2016) and state boundaries (Hutchison & Templeton, 1999), even within reintroduced species (Puckett et al., 2014). The observation of biogeographic breaks at political boundaries not designated based on the geology of the region (such as major rivers) suggests that differences in historic land use may have affected upland terrestrial habitats. After European settlement, the CIH experienced extensive habitat alteration (Cutter & Guyette, 1994;Nigh, 2004) that negatively influenced many endemic species through decreased dispersal among isolated glade habitats (Hutchison & Templeton, 1999) and increased siltification and eutrophication of streams leading to decreased reproductive success (Wheeler, Prosen, Mathis, & Wilkinson, 2003).
Previous genetic studies indicate that A. annulatum have been influenced by processes across evolutionary and ecological time frames. Using mitochondrial DNA (mtDNA) restriction fragment length polymorphisms (RFLPs) on samples collected across the Missouri portion of the species range, Phillips, Suau, and Templeton (2000) found three haplotypes with a south to north gradient of decreasing genetic diversity. The most common haplotype was found across all sampling sites, whereas the other two haplotypes were only found in southern Missouri, toward the center of the species distribution. Their results suggested that the fluctuating margins of temperate forest and prairie habitats during the Hypsithermal intervals of the Holocene (12 ka -present) resulted in repeated expansions and contractions of the northern periphery of the range and thus lower genetic diversity in the northern A. annulatum populations.
More recent studies using microsatellite loci observed two genetic clusters over an approximately 7,000 ha landscape (Burkhart et al., 2017;Peterman et al., 2015), which indicates the potential for dispersal limitation over more ecological time frames.
Herein, we employ nuclear microsatellite genotypes and mitochondrial DNA sequences to estimate the genetic diversity and differentiation across ecological and evolutionary time scales for A. annulatum. We increase the breadth of geographic sampling to encompass the entire distribution of A. annulatum which allowed us to characterize the patterns of genetic population structure and diversity of A. annulatum across its distribution, infer phylogeographic breaks throughout the range, and infer the species demographic history. We hypothesized that (a) overall genetic diversity across the range will be low as a result of the range contraction during the Pleistocene; (b) the southern Ozark Highlands or Ouachita Mountains will have higher genetic diversity as would be expected given a northern recolonization of the CIH following the Pleistocene F I G U R E 1 Known occurrence locations and sampling sites for Ambystoma annulatum within the International Union for the Conservation of Nature (IUCN) distribution overlaid on a simplified 2011 Land Use, Land Cover layer. "Forest" habitats include broadleaf, mixed, and coniferous forests, and other shrub/scrub land use classes. "Herbaceous" habitats include agricultural fields, prairie, and other herbaceous habitats. Sampling sites are indicated by the large circles glacial cycles; and (c) gene flow will be limited as a result of land use changes post-European settlement (e.g., timber harvest, fire suppression, reforestation) given that A. annulatum is a forest-dependent species and is sensitive to landscape perturbations.

| Study area and sample collection
We surveyed breeding ponds where A. annulatum are known or presumed to occur across the Ozark Highlands and Ouachita Mountains in Missouri, Arkansas, and Oklahoma ( Figure 1). These regions are dominated by temperate broadleaf and mixed forests (Ouachita Ecoregional Assessment Team, 2003;Ozarks Ecoregional Assessment Team, 2003). To balance representing the maximal A. annulatum genetic diversity across the range and the number of sampling sites, we sampled up to 20 individuals per pond from up to two ponds per county that were separated by >1.5 km. In total, we in Ambystoma spp. salamanders shows that effective population size and genetic diversity does not significantly differ among life stages (Peterman, Brocato, Semlitsch, & Eggert, 2016) and genetic inference is largely concordant among breeding seasons across the same landscape (Burkhart et al., 2017;Peterman et al., 2015). All

| Molecular analysis
Microsatellite DNA-We extracted DNA from tissue samples using InstaGene (Bio-Rad) following the protocols outlined by Peterman, Connette, Spatola, Eggert, and Semlitsch (2012). We genotyped individuals from 17 of the 19 sampling sites at 14 microsatellite loci (Peterman et al., 2012(Peterman et al., , 2013. We did not include the Lincoln and Le Flore sites because they had insufficient sample sizes for microsatellite analyses (n = 6 and n = 2, respectively). Forward primers were fluorescently labeled and amplified in two multiplexed polymerase chain reactions (PCR; Peterman et al., 2013) following the protocol outlined in Peterman et al. (2012). We sized the amplified products on an ABI 3730xl DNA Analyzer (Applied Biosystems) using Liz 600 size standard at the University of Missouri DNA Core Facility and scored genotypes using GeneMarkerv.1.97 (Softgenetics). Every multiplex plate included a negative control to detect contamination and a positive control to standardize scoring among plates. We genotyped 10% of samples twice to ensure quality of genotype calls.
We tested for the presence of full siblings using Colony v2.0.5.9 (Jones & Wang, 2010). Within Colony, we set male and female mating to polygamous without inbreeding and ran a long run with full likelihood, high precision, and no sibship prior. We tested for significant linkage disequilibrium (LD) and deviations from Hardy-Weinberg Equilibrium (HWE) with Genepop on the Web (Raymond & Rousset, 1995;Rousset, 2008). All tests were conducted with 1,000 dememorization steps and 100 batches with 1,000 iterations per batch, and we assessed significance using a Bonferroni correction for the number of comparisons (Rice, 1989). We tested for the presence of null alleles using the "PopGenReport" package (Adamack & Gruber, 2014) for R v3.3.3 (R Core Team, 2018). We used hp-RaRe (Kalinowski, 2005) to calculate rarefied allelic richness (A R ) and GenoDive v2.0 (Meirmans & Van Tienderen, 2004) to calculate observed and expected heterozygosity (H E and H O ), pairwise fixation index values (F ST ; using 9,999 permutations), and inbreeding coefficients (F IS ).
We tested for isolation by distance (IBD) using the "vegan" package (Oksanen et al., 2016) for R, and we tested for significant differences in our genetic diversity estimates between genetic clusters and core versus periphery populations using Kruskal-Wallis tests in R. Finally, we tested for the presence of a recent genetic bottleneck using the program BottleneCk v1.2.02 (Cornuet & Luikart, 1996) using the Wilcoxon sign rank test for heterozygote excess under 10,000 iterations of a two-phase mutation model (TPM) with 95% single-step mutations and 5% multistep mutations and a variance among multiple steps of approximately 12 (Piry, Luikart, & Cornuet, 1999).
To assess population genetic clustering, we used the Bayesian clustering algorithms implemented in the program StRuCtuRe v2.3 (Pritchard, Stephens, & Donnelly, 2000). In StRuCtuRe, we used 2.5 × 10 5 burn-in steps followed by 5.0 × 10 5 MCMC iterations for ten replicates of each K = 1-20 under an admixture model with correlated allele frequencies and no location prior. We estimated the most likely number of genetic clusters inferred by StRuCtuRe using the ΔK algorithm of Evanno, Regnaut, and Goudet (2005) as implemented in StRuCtuRe haRveSteR v0.6.94 (Earl & vonHoldt, 2012), averaged the model outputs for the ten technical replicates using Clumpp v1.2 (Jakobsson & Rosenberg, 2007), and visualized the results using DiStRuCt v1.1 (Rosenberg, 2004). We then tested for hierarchical genetic substructure within each putative cluster in a separate analysis for ten replicates of K = 1-15 following the same protocols as the initial analyses. We conducted a hierarchical analysis of molecular variance (AMOVA) in aRlequin v.3.5 (Excoffier & Lischer, 2010) with data grouped into STRUCTURE clusters.
Mitochondrial DNA-We randomly selected a subset of 77 individuals for mtDNA sequencing analysis from 11 populations distributed relatively evenly across the range of A. annulatum to maximize the amount of genetic diversity represented (Figure 1). We amplified a 900 basepair (bp) fragment consisting of the 3′ end of the cytochrome b gene, tRNA Thr , tRNA Phe , the intergenic space (IGS), and tRNA Pro to the 5′ end of the D-loop (control region) using the THR and DL1 primers (Shaffer & McKnight, 1996). This fragment has been found to be informative in previous phylogeographic studies of Ambystoma spp. (Church et al., 2003;Pauly, Bennett, Palis, & Shaffer, 2012;Pauly, Piskurek, & Shaffer, 2007;Shaffer & McKnight, 1996;Zamudio & Savage, 2003). We conducted the PCR in 25 µl volumes using 1X AmpliTaq Gold PCR Buffer (Applied Biosystems), 200 uM dNTPs, 2 mM MgCl 2 , 0.8 mM BSA, 0.4 uM of each PCR primer, 1 U AmpliTaq Gold (Life Technologies), and ~15 ng of DNA using thermocycler conditions set by Shaffer and McKnight (1996) with a 10-min initial denaturation step at 95°C and an annealing temperature of 55°C. We bi-directionally sequenced PCR products at the University of Missouri DNA Core Facility on an ABI 3730 DNA Analyzer (Applied Biosystems).
We edited and aligned mitochondrial sequences in GeneiouS v8.1.9 (Kearse et al., 2012). To verify the concordance of sequence calling, we used the ClustalW method (Thompson, Higgins, & Gibson, 1994) to align raw sequence data for each individual. After verifying individual sequences, we aligned consensus sequences for all individuals using a ClustalW alignment. To guarantee equivalent comparison of haplotypes across individuals, we trimmed our alignments to a 776 bp segment of our 900 bp region which included 233 bp of the intergenic spacer, the tRNA Pro , and 470 bp of the Dloop. We quantified the number of distinct haplotypes using Fabox v1.41 (Villesen, 2007;accessed 11 December 2017). We visualized haplotype networks in popaRt v1.7 (Leigh & Bryant, 2015) using the median-joining network method setting epsilon = 0 (Bandelt, Forster, & Rohl, 1999). Finally, we calculated haplotype diversity (h), nucleotide diversity (π), and Φ ST in aRlequin.

| Demographic modeling
We inferred A. annulatum demographic history in the CIH using approximate Bayesian computation (ABC) implemented in DiyabC v2.1.0 (Cornuet et al., 2014) using our microsatellite data. We tested a single model to estimate the divergence time between the northern and southern CIH clusters. In this model, the northern CIH and the southern CIH merge t 1 generations before present. Across all simulations, we allowed effective population size (N e ) of each cluster to vary as amphibian population sizes are highly variable in both time and space (Werner, Skelly, Relyea, & Yurewicz, 2007). We simulated 10 7 data sets drawing from uniform prior distributions of N e and divergence times. We calculated A, H O , mean allele size variation, pairwise F ST (Weir & Cockerham, 1984), and genetic distance between populations (δμ) 2 (Goldstein, Linares, Cavalli-Sforza, & Feldman, 1995) for each simulation to compare with the observed data set in the model selection.
Our simulations assumed that microsatellites followed a generalized stepwise mutation model (Estoup, Jarne, & Cornuet, 2002) with the mean mutation rate drawn from a uniform prior distribution between 10 -5 and 10 -3 . As suggested in the DIYABC manual, we allowed microsatellites to vary up to 40 alleles and we used default settings for all other microsatellite mutation parameters. We estimated the posterior distributions of N e and divergence time using a local linear regression of the closest 1% of simulated data sets to the observed data following logistic transformation of parameter values (Cornuet et al., 2014(Cornuet et al., , 2008. Finally, we estimated bias in the parameter estimates using DIYABC's built-in function and simulated 500 pseudo-observed data sets from the 90% highest posterior density distributions for each parameter. DIYABC computes the true value and compares the estimated values to the variance in simulated data sets and reports root of the relative mean integrated square error (RRMISE), average relative bias, and factor 2 score of the mode. The factor 2 score is the proportion of simulated data sets for which the estimate is at least half and at most twice the true value and values closer to 1.0 are desirable.

| Microsatellite DNA
We found no individuals within the same sampling location that had a probability of being related as full siblings >95%. Similarly, we observed no sampling site or loci that deviated significantly from expectations under HWE and no loci were significantly linked. Further, we did not detect the presence of null alleles within our data set.
Thus, all individuals and loci were retained for downstream analyses.
We observed two well-supported genetic clusters (ΔK = 154.49) in our StRuCtuRe analyses (Figure 2). The northern CIH cluster includes every Missouri sampling location and the southern CIH cluster includes all Arkansas and Oklahoma sampling sites (Figure 2). In our tests for hierarchical structure, we observed 11 well-supported clusters within the northern cluster (ΔK = 81.22) and two well-supported clusters within the southern cluster (ΔK = 991.54; Figure 2).
Within the northern cluster, average H O was 0.603 ± 0.076 (SD) and H E was 0.625 ± 0.082 (SD), whereas in the southern cluster diversity average H O was 0.623 ± 0.075 (SD) and average H E was 0.686 ± 0.069 (SD) ( Table 1). Within the northern cluster, the average allelic richness was 5.01 ± 1.44 (SD) and average rarefied allelic richness was 2.16 ± 1.18 (SD), whereas the average allelic richness TA B L E 1 Genetic diversity for Ambystoma annulatum populations using 14 microsatellite loci including sample size (N), mean number of alleles per locus (A), mean rarefied allelic richness (A R ), observed (H O ) and expected heterozygosity (H E ), and inbreeding coefficient (F IS ). Cluster membership is based upon STRUCTURE assignments shown in Figure 2 Population  Table 1) and (b) from each of the two supported genetic clusters identified in our full analysis. Spatial arrangement of genetic structure for A. annulatum across its distribution (c). Ellipses color corresponds to the admixture assignment scores from the full analysis (a) and pie colors correspond to admixture assignment scores from the two independent analyses for hierarchical substructure of individuals within each of the two presumed genetic clusters (b) was 5.64 ± 1.03 (SD) and average rarefied allelic richness was 1.68 ± 0.07 (SD) in the southern cluster (Table 1). Finally, within the northern cluster the average F ST was 0.191 ± 0.074 (SD) and F IS was 0.034 ± 0.044 (SD) whereas average F ST was 0.190 ± 0.065 (SD) and F IS was 0.091 ± 0.057 (SD) in the southern cluster (Table 1). For all genetic diversity metrics, neither group had significantly greater diversity (all p > .05). The AMOVA results indicated significant subdivision among groups (and among populations within groups (Table S1).
Across the distribution of A. annulatum, all but one pairwise comparison of sampling locations revealed significant genetic differentiation (Table 2); the St. Louis #1 and St. Charles sampling locations in northeastern Missouri (24.6 km apart) were not significantly different (Table 2). Tests for IBD indicate that distance was a significant predictor of genetic differentiation across the entire range (Mantel's r = .363, p = .001; Figure S1); however, distance was not a significant predictor of genetic structure within each cluster (northern clus-  (Table S3).

| Mitochondrial DNA
We observed nine mtDNA haplotypes based on seven polymorphic basepairs across the 776 bp mtDNA sequence in the 77 individuals sampled (Accession No.: MN242388-MN242396;  Table 3 and Figure 3). No clear patterns emerged in our median-joining haplotype network (Figure 3) as most haplotypes were 1-3 bp diverged from haplotype A, the most common haplotype (Table 3).
Across all populations, haplotype diversity (h) ranged from 0.000 to 0.694 and nucleotide diversity (π) ranged from 0.000 to 0.002 (Table 3). Populations toward the edge of the range were fixed for a haplotype whereas populations closer to the range center had haplotype diversity between 0.250 and 0.694 (Table 3). Pairwise Φ ST values based on mtDNA range from −0.040 to 1.00 (Table S4) and the average pairwise Φ ST was 0.542 ± 0.370 (SD). Tests for IBD indicate that distance was a significant predictor of genetic differentiation based on mtDNA haplotypes (Mantel's r = .362; p = .032; Figure S1); however, distance was not a significant predictor of genetic differentiation in either cluster (northern CIH: Mantel's r = −.315; p = .950; southern CIH: Mantel's r = .390; p = .500, respectively).

| Demographic modeling
We estimated that the divergence time between the northern and the southern CIH clusters was approximately 1,381 generations (90% HPD = 771-5,966 generations; Table 4). Given an average generation time of 2 years , we estimated that these

| D ISCUSS I ON
The patterns of A. annulatum genetic diversity and differentiation relate to both geologic and anthropogenic processes. Overall, we observed two strongly supported nuclear genetic clusters that correspond to a northern Ozark cluster and a southern Ozark and Ouachita Mountain cluster. Their estimated divergence time was 2,762 years before present (90% HPD = 1,542-11,932 years ago), corresponding to the early to middle Holocene (Pielou, 2008). Effective population size estimates were higher in the southern CIH cluster (N e = 8,201) than the northern CIH cluster (N e = 4,148). We acknowledge that the effective population sizes may be overestimated since we observed substructure within A. annulatum; simulation studies have shown that estimates of effective population size are higher in populations with genetic substructure than in randomly mating populations of similar size (Petit & Excoffier, 2009). For A. annulatum specifically, breeding census population sizes ranged from 336 to 1,971 adults at a single pond during a single breeding season . Similarly, estimates of effective population size ranged from 2,596 to 34,310 for two congeneric species at a single sampling location (Nunziata, Lance, Scott, Lemmon, & Weisrock, 2017). Taken altogether, the higher effective population size and haplotypic diversity observed in the southern Ozarks suggests that the south-central CIH was the center of diversity for A. annulatum and potential location of an ancestral population, during the Pleistocene and climatic fluctuations during the Hypsithermal intervals of the Holocene.
Further, we observed significant pairwise F ST values between our northern and southern CIH clusters (Table 2 and Table S4) and low migration rates among populations (Table S2) we found concordance in the approximate location of the A. annulatum phylogeographic break with those of aquatic and terrestrial fauna that exist between the northern and southern Ozarks (Hardy et al., 2002;Hutchison & Templeton, 1999;Phillips, Fenolio, Emel, & Bonett, 2017;Ray, Wood, & Simons, 2006;Routman et al., 1994).
Within aquatic taxa, most studies found significant genetic differentiation between populations located in watersheds draining southward into the Arkansas River and watersheds draining northward in the Missouri or eastward into the Mississippi (Figure 2; Hardy et al., 2002;Ray et al., 2006). In some instances, such as hellbender salamanders (Cryptobranchus alleganiensis spp.), genetic differentiation between watersheds is pronounced enough to merit subspecific status (Crowhurst et al., 2011;Routman et al., 1994).
Further, cave-adapted species, such as grotto salamanders (Eurycea F I G U R E 3 Spatial arrangement of mtDNA haplotypes by sampling location for Ambystoma annulatum. For each color, the proportion of the pie corresponds to the frequency of the haplotype in that population as outlined in Table 3 Note: Bias estimates for the root of the relative mean integrated square error (RRMISE), average relative bias, and the factor 2 score of the mode of the posterior distribution are also included. spelaea), show similar phylogeographic breaks between the northern and southern Ozarks that were attributed to the palaeodrainages, contemporary drainages, and subplateaus in the Ozarks (Phillips et al., 2017). Similarly, some terrestrial taxa exhibit similar patterns of genetic diversity despite these organisms not being constrained to riparian networks (Hutchison & Templeton, 1999).
We observed nine unique mtDNA haplotypes across the range of A. annulatum (four private to the northern CIH, four private to the southern CIH, and one shared), one of which was found within most Missouri sampling sites (Table 3 and Table S2), in comparison with the three previously observed haplotypes found in Missouri using RFLP data (Phillips et al., 2000).  (Table 3 and Figure 3). However, Haplotype C was only 1 bp different from the most common haplotype (Haplotype A) and was intermediate to haplotypes that were found within populations in the northern CIH cluster (Figure 3).
Taken altogether, we suggest that A. annulatum populations experienced a bottleneck during Pleistocene glacial cycles and Hypsithermal intervals of the Holocene, with populations contracting into two refugial areas in the south-central Ozarks and differentiated as they re-colonized the peripheries of their contemporary range during postglacial dispersal. We observed a trend toward more haplotypes and higher genetic diversity in the center of the species distribution than toward the periphery ( Post-European settlement, the CIH experienced extensive timber harvest (Cunningham, 2007;Flader, 2004), land conversion into agriculture (Jacobson, 2004;Jacobson & Primm, 1997), channelization and degradation of riparian habitats (Jacobson, 2004;Jacobson & Primm, 1997;Nigh,2004;Swift, 1984) and changes in the natural wildfire regime (Cutter & Guyette, 1994;Jacobson, 2004), and more recent reforestation (Flader, 2004). Despite the general records of anthropogenic disturbance in this area, there is a paucity of detail regarding the spatial and temporal scale of disturbance which precludes the ability for detailed landscape genetic analyses. However, these landscape perturbations likely decreased the amount of suitable aquatic breeding and terrestrial nonbreeding habitat available for A. annulatum. Thus, with decreased habitat suitability and low levels of A. annulatum gene flow, genetic drift in small isolated populations may have led to the observed patterns of genetic differentiation and substructure as suggested for C. c. collaris (Hutchison & Templeton, 1999). Although the high rates of anthropogenic disturbance in the CIH have only occurred over the last 200 years, for a species with an average generation time of 2 years, such as A. annulatum , up to 100 generations may have been sufficient for the observed patterns of genetic differentiation to emerge post-European settlement.
Worldwide, amphibian populations are experiencing declines as a result of increased rates of anthropogenic disturbance (i.e., habitat alteration, degradation, and loss), global climate change, and increased prevalence of disease (Wake & Vredenburg, 2008). Despite their designation as a species of conservation concern, A. annulatum populations are presumed to be relatively stable by state wildlife agencies and thus our data provide a record of the genetic diversity that is essential for evaluating the potential genetic consequences of declining populations. Although our data did not support the presence of a recent genetic bottleneck, the observed genetic substructure and differentiation present across the distribution of A. annulatum indicates that these populations have experienced population declines from presettlement population sizes. Further, our inability to attain sufficient sample sizes from historic localities ( Figure 1) suggests that more detailed study of A. annulatum occupancy, abundance, and niche requirements across its distribution would provide important data for the management of this species.

ACK N OWLED G M ENTS
We

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest associated with this work.

AUTH O R CO NTR I B UTI O N S
JJB, EEP, LSE, and RDS conceived the ideas for this project; JJB, CJK, and CNS conducted the research and data analyses; LSE and RDS directed the project; and JJB led the writing of the manuscript with assistance from EEP, CJK, CNS, and LSE.

DATA AVA I L A B I L I T Y S TAT E M E N T
All mtDNA sequences are available on GenBank (Accession No.: MN242388-MN242396). Locality data are only provided to the county level as exact coordinates are sensitive due to the conservation status of this species; however, further information is available from the authors upon request.