High inbreeding and low connectivity among Ambystoma texanum populations in fragmented Ohio forests

Abstract Habitat loss and fragmentation negatively impact the size and diversity of many natural populations. Woodland amphibians require connected aquatic and terrestrial habitats to complete their life cycle, and often rely on metapopulation structure for long‐term persistence. Wetland loss and deforestation fragment amphibian populations, which may result in population isolation and its negative effects. The aim of this research was to analyze the population genetic structure of small‐mouthed salamanders (Ambystoma texanum) in western Ohio, where agriculture is now the dominant land use. Salamander tail tissue was collected from eight breeding pools. Three pools occur in the same forest; the other five are in forest patches at distances ranging from 250 m to 20 km from one another. Eight microsatellite loci were amplified by PCR and genotyped for allele size. Observed heterozygosities were lower than expected in all sampled populations; the two most isolated sites (Ha1, Ha2) had the highest inbreeding coefficients. Ha2 also had the lowest mean number of alleles and was found to be genetically differentiated from populations to which our data analysis indicates it was historically connected by gene flow. The most distant site (Ha1) had the highest number of private alleles and showed genetic differentiation from other populations both historically and currently. Geographic distance between pools was strongly correlated with the number of private alleles in a population. The results suggest that population isolation results in decreased genetic diversity and that a breakdown of metapopulation structure due to landscape change may contribute to differentiation between once‐connected populations.

genotyped for allele size. Observed heterozygosities were lower than expected in all sampled populations; the two most isolated sites (Ha1, Ha2) had the highest inbreeding coefficients. Ha2 also had the lowest mean number of alleles and was found to be genetically differentiated from populations to which our data analysis indicates it was historically connected by gene flow. The most distant site (Ha1) had the highest number of private alleles and showed genetic differentiation from other populations both historically and currently. Geographic distance between pools was strongly correlated with the number of private alleles in a population. The results suggest that population isolation results in decreased genetic diversity and that a breakdown of metapopulation structure due to landscape change may contribute to differentiation between once-connected populations.

K E Y W O R D S
Ambystoma texanum, habitat fragmentation, landscape genetics, microsatellites, mole salamanders, population structure

| INTRODUCTION
Habitat destruction and fragmentation are the leading cause of species declines and extinctions worldwide (Fahrig, 2003). Habitat patches may contain small populations that are isolated from conspecifics in patches that are too far away or unreachable due to migration barriers.
Isolated populations are at great risk of extinction from either stochastic or deterministic causes, as their genetic diversity and biological fitness are reduced over time (Allentoft & O'Brien, 2010).
The eastern North American landscape has been significantly altered in the last two hundred years, mainly through deforestation. In western Ohio today, small forest patches are common in an agricultural landscape (Gordon, 1969). These patches often contain vernal pools that serve as breeding habitat for woodland amphibians. Loss of both terrestrial and aquatic habitats for amphibians has occurred, as wetland drainage accompanied forest clearing for agriculture (Dahl, 1990).
Many amphibian species in this region are in decline (Lannoo, 1998), and spatial isolation due to fragmentation may be a serious threat to their survival as woodland amphibians prefer not to traverse nonforested landscapes (Pittman & Semlitsch, 2013;Rothermel, 2004;Rothermel & Semlitsch, 2002).
Amphibian extinctions are occurring at an alarming rate around the globe (Blaustein, Wake, & Sousa, 1994;McCallum, 2007), and the principal driver of extinction is habitat change (Brooks et al., 2002). Many amphibian species have metapopulation structure that requires habitat connectivity between geographically close breeding populations (Marsh & Trenham, 2001). Low vagility, and sensitivity to migration barriers produced by habitat fragmentation and other anthropogenic development, limits migration between breeding populations (Semlitsch & Bodie, 2003). Amphibian species conservation must utilize molecular tools to understand the genetic impact of habitat change (McCartney-Melstad & Shaffer, 2015).
Identifying the geographic scale at which gene flow is interrupted by habitat fragmentation is crucial for prescribing conservation and management plans for amphibians. Scale can vary widely depending on the species' sensitivity, migratory ability, and the effective population sizes of isolated populations (Jacobs & Houlahan, 2011;Richardson, 2012;Scott, Komoroski, Croshaw, & Dixon, 2013). A consensus of landscape genetics studies suggests that a species-specific approach is necessary for assessing the effect of habitat change on population genetics (Storfer, Murphy, Spear, Holderegger, & Waits, 2010). Factors including geographic distance between breeding pools, intervening land use, and occurrence of migration barriers are important in influencing the connectivity among populations of amphibians at both a landscape and fine scale, and may impact the genetic diversity of each population.
Mole salamanders (Ambystoma) inhabit terrestrial areas that contain ephemeral breeding pools, to which adults exhibit high fidelity (Shoop, 1968;Shoop & Doty, 1972). Forest fragmentation may disrupt migration between breeding populations due to their low vagility. A study of the spotted salamander (A. maculatum) in northeastern Ohio found that gene flow may occur in a highly fragmented landscape, possibly along riparian zones (Purrenhage, Niewiarowski, & Moore, 2009), while others have found genetic differentiation among breeding populations in the same forest patch (Bartoszek, 2009) and across a small migration barrier such as a railroad track (Bartoszek & Greenwald, 2009).
The aim of this study was to investigate the impact of forest fragmentation on the genetic diversity and population structure of the small-mouthed salamander (Ambystoma texanum) at both the landscape and fine spatial scales in western Ohio. This species' range extends throughout most of eastern North America, where it resides in forests that contain vernal pools (Petranka, 1998). Small-mouthed salamanders may disperse up to 125 m from their natal pool (Williams, 1973). This species is one of the smallest Ambystomids (Petranka, 1998), and its dispersal ability may be more limited than larger species (Semlitsch & Bodie, 2003).
Given the highly fragmented nature of forest patches in western Ohio, the main concern is that a lack of connectivity between breeding populations in separate forest patches eliminates metapopulation structure, results in reduced genetic diversity, and may lead to local extinctions. In the mid-1800s, continuous forest in Hardin County, Ohio, was cleared for farming. Since settlement, the county is heavily agricultural (USDA 2009), yet contains multiple forest patches with vernal pools (Figure 1). The spatial arrangement of the forest patches in the study area provided an opportunity to approach this question from both a landscape and fine scale. Prior studies have shown the importance of both of these scales for Ambystoma salamanders (Greenwald, Gibbs, & Waite, 2009;Porej, Micacchion, & Hetherington, 2004), and here, data are presented from populations within the same forest patch and from populations at increasing distances away from this focal forest patch.
It was hypothesized that: (1) Genetic diversity would be lower (fewer alleles and higher inbreeding coefficients) in isolated populations, as compared to those in forest patches with more than one breeding pool. (2) Populations in separate forest patches would show genetic differentiation from one another, while those in the same patch would not. (3) Migration is not occurring between populations in separate forest patches that were historically connected by forested habitat and presumably gene flow.

| Study sites and field sampling
Salamanders were collected from eight breeding pools in five forest patches in rural Hardin County, Ohio ( Figure 1). Three pools are within a 283-ha forest (Lawrence Woods State Nature Preserve), and five additional pools are in separate forest patches at varying distances.
Landscape measures determined using Google Earth (Google Inc. 2009) were forest patch area, pool area, distance between forest patches, and distance between pools (Table S1). For each sampled pool, the average distance to all other study pools and distance to nearest sample pool were recorded.
Adult A. texanum were captured in screen funnel traps placed in breeding pools in March 2009 and 2010. Tail tissue from up to 20 individuals was preserved in ethanol, and salamanders were released back into the pools. For sites without 20 adult samples, larvae were used so that we had a total of 20 individuals for each pool.
Ambystoma texanum larvae were collected with dip nets during the summer of 2008.

| DNA extraction and marker amplification
Genomic DNA was extracted and purified from tail tissue using standard proteinase-K digestion and Wizard SV Genomic DNA Purification System (Promega). Ambystoma texanum larvae are difficult to distinguish from unisexual Ambystoma larvae (Brandon, 1961), which are common in some sample pools (Pfingsten & Downs, 1989). Larval species identification was confirmed using a molecular technique previously described (Rhoads, Krane, & Williams, 2010).
Eight polymorphic microsatellite loci were used as genetic markers in this study (Table S2); seven of the loci were characterized for A. texanum (Williams & Dewoody, 2004), and one locus was originally identified in A. jeffersonianum (Julian, King, & Savage, 2003) and found to cross-amplify in A. texanum (J. P. Bogart, personal communication).
The microsatellite loci were amplified by polymerase chain reaction, and fragment analysis of the products was carried out on an ABI 3100 and compared to a ROX 500 standard. PCR product size was determined using Genescan v.3.1 and Genotyper v.2.5 (Applied Biosystems).

| Null alleles and relatedness
Each microsatellite locus was tested for conformance to Hardy-Weinberg equilibrium and for linkage disequilibrium using Genepop v.4 (Rousset, 2008) (Table S2) and were detected to have significant heterozygote deficiency in Genepop. Genotypes were adjusted in Microchecker according to Chakraborty, De Andrade, Daiger, and Budowle (1992).
To rebut the possibility that larval samples from the same pool were siblings, Queller and Goodnight's (1989) method of calculating pairwise relatedness (r qg ) between individuals within populations was performed in Genalex v.6.3 (Peakall & Smouse, 2006). Bootstrap resampling (1,000 bootstraps) provided 95% confidence intervals for r qg for each population. A permutation test (999 permutations) for significant differences among mean population relatedness was performed, and 95% confidence intervals were derived for the expected range of r qg under the null hypothesis of random mating. This method determines whether populations with larval samples showed higher within-population relatedness than those with only adult samples, which would be expected if siblings had inadvertently been collected.

| Within-population diversity and population structure
Genalex was used to determine population genetic measures for each population including mean number of alleles, number of private The program structure 2.3.1 (Pritchard, Stephens, & Donnelly, 2000) was used to assign individuals to inferred genetic clusters based on multilocus genotypes. Ten independent runs of K = 1-8 were performed with 10,000 Markov Chain Monte Carlo (MCMC) repetitions and a 10,000 burn-in period assuming no admixture and independent allele frequencies and using sampling locations as priors. K was identified as the mean maximal value of ln P(D) and confirmed with calculation of ΔK (the second-order rate of change between successive K values) (Evanno, Regnault, & Goudet, 2005). Once the most likely number of genetic clusters was identified, a final run of K was conducted with 10,000 MCMC repetitions and a 10,000 burn-in period with the same parameter set.
In order to determine the amount of variation explained by the genetic clusters identified in structure, a hierarchal AMOVA (analysis of molecular variance) was performed (999 permutations) in Genalex with populations assigned to regions based on the genetic clusters (region 1-Ha1, region 2-Ha2, region 3-all other sites). The F ST analog Φ ST , which suppresses within-individual variation, was used in order to detect subtle regional population structure.

| Population size and migration rates
The Garza and Williamson (2001)  Maximum likelihood analysis was performed in MiGrate (Beerli, 2008) to estimate both migration rates and effective population sizes.
The first run consisted of 10 short chains (50,000 sampled, 100 recorded) and one long chain (500,000 sampled, 5,000 recorded). Fourchain adaptive heating at default temperatures (1, 1.5, 3, and 10,000) was used to increase efficiency. Two additional runs consisted of 10 short chains (10,000 sampled, 500 recorded) and three long chains (100,000 sampled, 5,000 recorded), combining the long chains and no adaptive heating.
Output from MiGrate gives ϴ values that can be converted to effective population size (N e ) from the equation ϴ = 4N e μ. N e was calculated assuming a microsatellite mutation rate (μ) of 5 × 10 −4 (Garza & Williamson, 2001).
In addition to comparing bidirectional migration rates (M) for all population pairs, rates were averaged for groups of populations in two ways. First, pools that are in the same forest patch (Ha10 and Ha11, Lawrence Woods sites Ha3L1, 4, 5) were combined.
All of the populations that comprise region 3 from the AMOVA (Ha3L1, 4, 5 and Ha7, Ha10, Ha11) were then combined in order to assess migration rates among the genetic clusters identified in structure.

| Landscape genetics
A bivariate correlation matrix was created in pasw v.18 (2009) that included F IS values, number of private alleles (N P ), effective population sizes (N e ), forest patch area, pool area, distance to nearest sample pool, and average distance to other sample pools. Several linear univariate and multiple regressions with F IS , N P , or N e as the dependent variable and the landscape traits as independent variables were performed. Distance to nearest sample pool and average distance to other sample pools were significantly correlated (p < .001), so they were not included together in a multiple regression model. In Genalex, a Mantel test (999 permutations) was used to test for a relationship between geographic and population genetic distances.

| RESULTS
One hundred and sixty small-mouthed salamanders were genotyped from eight breeding populations at eight highly polymorphic microsatellite loci (Table S2). A total of 143 alleles were detected, with the number of alleles per locus ranging from 12 (AjeD422) to 26 (Atex65 and Atex87). These results expand the size ranges for some Atex loci from those previously published (Williams & Dewoody, 2004) and are the first published for A. texanum from AjeD422.

| Null alleles and relatedness
Microchecker found no evidence of large allele dropout in our data set, but did detect the possibility of null alleles for six loci due to homozygous excess (Table S2). Genepop corroborated in finding heterozygote deficiency for the same loci. The Microchecker adjusted genotypes with added dummy alleles increased the observed heterozygosities so that they were very close to the original expected heterozygosities.
Briefly, this method changes some homozygotes to heterozygotes by adding a dummy allele to some of the individuals in populations showing homozygous excess for a given locus. As explained in the Section 4, the data analysis was conducted with the original data set.
Increased relatedness (r qg ) was not detected in populations with larval samples (Fig. S1). Three populations (Ha2, Ha3L1, Ha11) had significantly higher r qg than expected under panmixia (p < .05), but two of those were comprised of mainly adult samples. Hardin 10, which was a mixture of adult and larval samples, had the only negative mean r qg (−0.011).

| Within-population diversity
A summary of genetic results by population is presented in Table 1.
The mean number of alleles per locus in each population ranged from 9.0 (Ha2) to 11.1 (Ha10). The overall mean observed heterozygosity

| Population structure
The greatest number of private alleles (12) was found in Ha1; all other populations had 4 (Ha7 and Ha10) or fewer (Table 1). Pairwise F ST values indicate that Ha1 has highly significant genetic differentiation relative to all other sites (p = .001), and Ha2 has significant genetic differentiation from all sites (p ≤ .009) except Ha10 (p = .148), which is the closest geographically (Table 2). There is very little to no genetic differentiation among Ha7, Ha10, and Ha11, and between those sites and the Lawrence Woods sites (Ha3L1, 4, 5). There is, however, significant genetic differentiation among the Lawrence Woods sites. The hierarchal AMOVA was performed with the populations assigned to three regions (Ha1, Ha2, all other sites). The AMOVA detected significant structure at the regional and within-population scale, but not among the populations in region 3 (Table 3).

| Population size
The effective population size (N e ) estimates are shown in Table 4

| Migration rates
Bidirectional migration rates generated by MiGrate for all sampled populations are presented in Table 5 Comparing migration rates at the forest patch scale shows that all rates to and from Ha1 are ≤0.91 (Table 6). Rates to and from Ha2 are all ≥0.97, and rates to and from Ha7 are all ≥0.92 (both excluding Ha1).
Ha 10,11 have migration rates ≥1.03 going to Ha2, Ha7, and the Ha3 sites; this patch has a rate coming from Ha2 of 1.28 and from Ha3 sites of 0.97. The Lawrence Woods sites (Ha3) have migration rates that range from 0.97 (to Ha 10,11) to 1.28 (from Ha2). T A B L E 2 Pairwise population F ST values a (below diagonal) and geographic distances in km (above diagonal) Migration rates were also compared among the regions from the AMOVA (region1-Ha1, region 2-Ha2, region 3-all other sites) ( Table 6). The migration rates between regions 1 and 2 (Ha1 to Ha2: 0.56, Ha2 to Ha1: 0.73) and between regions 1 and 3 (Ha 1 to region 3: 0.70, region 3 to Ha1: 0.77) were lower than between regions 2 and 3 (region 2 to region 3: 1.09, region 3 to region 2: 1.20).

| Landscape genetics
The univariate regression analysis found no significant relationship between any of the landscape variables and F IS values, or between number of private alleles and either forest patch area or pool area (  Figure 3).
The regression analysis using N e as the dependent variable and landscape traits as independent variables found that pool area has a significant positive correlation with N e (r 2 = .544, p = .037), but forest patch area does not (r 2 = .013, p = .788) ( distance to other sample pools and distance to nearest sample pool were significantly correlated with N e (r 2 = .673, p = .013 and r 2 = .682, p = .012, respectively). Pool area was combined with each distance measure into two separate multiple regression models, and both models showed improved fit (r 2 = .805, p = .017 and r 2 = .804, p = .017).

| DISCUSSION
The There are two possibilities that could explain the low observed heterozygosities in this data set: either there is a high frequency of null alleles in the loci or inbreeding has resulted in greater homozygosity than expected given the number of alleles present in the populations.
It is expected that as a population is reduced in size, heterozygosity decreases more rapidly than the loss of alleles due to genetic drift (Hartl & Clark, 1997). Reproduction among close relatives results in an increase in homozygosity and can be detected by calculating the inbreeding coefficient for the population. The low frequency of potential homozygous null genotypes in our data set (no PCR product -see Table S2) suggests that null alleles are not present at a frequency high enough to be driving the low observed heterozygosities and high inbreeding coefficients. Furthermore, departure from Hardy-Weinberg equilibrium in six of the eight genetic markers suggests that the pattern of nonrandom mating occurs across loci, and not specifically in one or two loci, in which case null alleles would be a more plausible explanation.

It was hypothesized that breeding populations in forest patches
with only one pool would show signs of decreased genetic diversity such as fewer alleles and lower observed heterozygosity. This was expected to be especially prominent in patches that appear isolated, Woods are all about 900 m from each other, which is a substantial distance for A. texanum to disperse, even within forest habitat. However, the presence of other unsampled pools within Lawrence Woods could serve as stepping stones that facilitate gene flow between our sites (Kimura & Weiss, 1964). This level of inbreeding in the Lawrence Woods sites may provide a baseline for population structure that is naturally present even in Ambystoma breeding populations that are in relatively close proximity to each other in forested habitat.
The mean observed heterozygosity was 0.543, which was very similar to the value found by Zamudio and Wieczorek (2007) of 0.534 for populations of A. maculatum in New York. Their study also detected higher than expected relatedness within ponds, although the range of r qg (0.02 to 0.17) was higher than the r qg values in this study (−0.011 to 0.040).
The measure of relatedness is determined by comparing each population to the entire group of populations and therefore cannot be directly compared between studies. If a group of populations contains one or more populations with higher relatedness than the others, then they will have higher r qg values and the range of r qg values will be larger. The populations in this study may have relatedness measures that are more similar and therefore have r qg values that are lower and within a smaller range.
T A B L E 6 Migration rates (M) a from MiGrate between forest patches and between regions from AMOVA In the year prior to sampling, there had been development within this forest patch between the two pools (personal observation), but it is unlikely that the disruption of gene flow between them would be detected this quickly. There are a few private alleles in the Lawrence Woods pools, which supports the F ST values in indicating population structure within this forest patch.
The identification of genetic clusters by structure largely corroborates the F ST data and supports the idea of past gene flow among the Lawrence Woods sites and Ha7, 10, and 11. All of these sites are included in a single cluster that is separate from Ha1 and Ha2. The AMOVA found significant structure at the regional level (regions were based on the genetic clusters identified in structure), while no significant variation was explained at the population level for region 3. So while there may be some population structure among the Lawrence Woods sites, it is not great enough for structure to be able to assign stone model prior to fragmentation between the forest patches that are left today (Kimura & Weiss, 1964 MiGrate also provided estimates of historical effective population sizes that may be most useful for comparisons between populations and are not necessarily accurate current population sizes. The range of effective population sizes calculated from the ϴ values generated by MiGrate is reasonable (Pfingsten & Downs, 1989 Population size is important in determining whether an isolated population's genetic diversity is depressed by inbreeding and genetic drift. According to Garza and Williamson (2001) historical N e estimate from MiGrate, which strongly suggests that this population was once larger in size and more genetically diverse than it is today.
Landscape factors that may limit population size for Ambystoma species include the size of the breeding pool and/or size of the forest patch. The regression analysis found that neither forest patch area nor pool area was correlated with inbreeding coefficients; however, effective population size did have a significant positive relationship with pool area-a result supported by other amphibian studies (Greenwald, Gibbs et al., 2009;Wang, Johnson, Johnson, & Shaffer, 2011). Overall, this suggests that while population size may be limited by breeding pool size, higher inbreeding coefficients are not directly related to the physical area of available breeding or terrestrial habitat. F IS may be more influenced by the isolation of populations. Landscape change, such as habitat fragmentation that disrupts connectivity among populations, reduces genetic diversity in amphibians (Greenwald, Purrenhage, & Savage, 2009;Greenwald, Gibbs et al., 2009;Richardson, 2012).
In the regression analysis, effective population size was also significantly correlated with the two landscape distance measures, although the relationship was positive and not expected. This positive correlation was driven by Ha1 having the highest N e and also being the most distant population from the others sampled. There was also a positive relationship between the landscape distance measures and inbreeding coefficients, but they were not significantly correlated.
Lack of significance was probably driven by Ha7, which had one of the lowest F IS values and is geographically distant from the other sites. The distance measures did, however, have a highly significant positive correlation with the number of private alleles in a population, which suggests that more isolated populations have less gene flow than among those closer together. This was also supported by the results of the Mantel test for isolation-by-distance that found that genetic distance between populations was positively related to their geographic distance (although it was not significant). Zamudio and Wieczorek (2007) found that among A. maculatum populations sampled in New York, those populations ≥4.8 km from each other were genetically independent. This may suggest that metapopulations for Ambystoma species operate within this scale, as populations within this distance were more similar. The importance of this scale for amphibians is supported by a study on Hyla arborea (Arens et al., 2006), but a study on A. maculatum and Rana sylvatica in a landscape with abundant forest and wetland habitat found high connectivity up to 17 km (Coster, Babbitt, Cooper, & Kovach, 2015). The results of this study show that Ha2 is today genetically differentiated from Ha11 and the Lawrence Woods populations, which are all within 5 km.
However, the historic migration rate estimates suggest that Ha 2 and Ha7 had gene flow at the same level that each had with the Lawrence Woods populations, even though Ha2 and Ha7 are almost 8 km apart. The estimates of historical gene flow in and around Lawrence

| CONCLUSIONS
Woods suggest that most of the populations sampled were once part of a metapopulation in a continuously forested habitat that had numerous pools with subtle population structure. Migration between populations in the metapopulation may have occurred at different rates between some of the populations, as data from this study suggest. Varying gene flow rates could have been influenced by geographic distance between populations and landscape traits, such as topography and soil moisture that impact vegetative composition.
It cannot be assumed that breeding populations in the same forest patch today necessarily have high rates of gene flow between them.
The Lawrence Woods sites illustrate this because there was greater genetic differentiation among those populations than between any of the Lawrence Woods sites and Ha7, Ha10, or Ha11. Conservation efforts should protect networks of breeding pools in close proximity, while considering natural dispersal routes and potential migration barriers, factors were not investigated in this study.
Landscape change, including deforestation and potentially loss of vernal pools, has disrupted the metapopulation structure detected to be present historically among our populations (excluding Ha1). Ha3L5 today is significantly differentiated from Ha2; but in the past, these two populations were connected by a level of gene flow higher than between Ha3L5 and the other Lawrence Woods populations we sampled, demonstrating that gene flow may have been highest between Lawrence woods pools and those that today are outside the forest.
Today, the landscape between these two sites is agricultural and likely inhospitable to migrating amphibians. Even among the sampled sites without significant pairwise F ST values (Ha10, 11, 7, and between those and Lawrence Woods), it is unlikely that gene flow between forest patches is occurring today because of both geographic distance and intervening land use, as supported by private alleles in Ha7 and Ha10. These populations may become genetically differentiated in the future, as mutations continue to introduce new alleles.

ACKNOWLEDGMENTS
Undergraduate students Andy Lewis, Sarah Timko, Mark Lauber, and Anna McCrate provided research assistance for this study.

AUTHOR CONTRIBUTIONS
E.R. and P.K.W. initiated and designed the project, and performed field sampling. E.R. performed the laboratory and data analyses, and wrote the first draft of the manuscript. C.K. advised on the laboratory and data analyses. All authors contributed to the final manuscript.

DATA ACCESSIBILITY
Microsatellite genotypes will be submitted to Dryad for data archiving.