Population structure, gene flow, and sex‐biased dispersal in the reticulated flatwoods salamander (Ambystoma bishopi): Implications for translocations

Abstract Understanding patterns of gene flow and population structure is vital for managing threatened and endangered species. The reticulated flatwoods salamander (Ambystoma bishopi) is an endangered species with a fragmented range; therefore, assessing connectivity and genetic population structure can inform future conservation. Samples collected from breeding sites (n = 5) were used to calculate structure and gene flow using three marker types: single nucleotide polymorphisms isolated from potential immune genes (SNPs), nuclear data from the major histocompatibility complex (MHC), and the mitochondrial control region. At a broad geographical scale, nuclear data (SNP and MHC) supported gene flow and little structure (F ST = 0.00–0.09) while mitochondrial structure was high (ΦST = 0.15–0.36) and gene flow was low. Mitochondrial markers also exhibited isolation by distance (IBD) between sites (p = 0.01) and within one site (p = 0.04) while nuclear markers did not show IBD between or within sites (p = 0.17 and p = 0.66). Due to the discordant results between nuclear and mitochondrial markers, our results suggest male‐biased dispersal. Overall, salamander populations showed little genetic differentiation and structure with some gene flow, at least historically, among sampling sites. Given historic gene flow and a lack of population structure, carefully considered reintroductions could begin to expand the limited range of this salamander to ensure its long‐term resilience.

. However, modern anthropogenic landscape changes have made it impossible for some species to naturally recolonize extirpated areas of their historic range and translocations may be necessary to reintroduce species with low vagility.
Together with population structure, estimates of contemporary and historic migration rates are important for estimating connectivity among breeding sites. Historically connected sites with little genetic structure are straightforward candidates for translocations and reintroductions. For species with population structure, moving animals to places where the species has been extirpated or where populations have been severely reduced becomes a more complex challenge. The International Union for Conservation of Nature (IUCN) defines translocations as the human-mediated movement of living organisms from one area, with release in another. It defines reintroductions as the intentional movement and release of an organism inside its indigenous range from which it has disappeared (IUCN, 2013). Several studies have tested the effectiveness of translocations and reintroductions in order to conserve endangered amphibian species. Eastern hellbenders (Cryptobranchus alleganiensis) have survived several years post-release in sites where they were translocated to supplement existing populations (Kraus et al., 2017;McCallen et al., 2018) and spotted salamanders (Ambystoma maculatum) have been reintroduced by moving egg masses into artificial ponds to accelerate restoration without compromising larval survival (Sacerdote, 2009). Overall, promoting dispersal has been a tool in amphibian conservation and may be used to bolster and expand populations of imperiled species.

The reticulated flatwoods salamander (RFS, Ambystoma bishopi)
is a federally endangered species with a highly fragmented range.
Historically, RFS occurred throughout the southeastern United States in fire-maintained longleaf pine (Pinus palustris) ecosystems (Palis, 1996;Petranka, 2010). Over the last several decades, this species has declined for many reasons including climate change, drought, continued loss of flatwoods habitat, and a change in fire regime from summer fires when breeding habitat is dry to prescribed fires in winter when habitat is wet (Bishop & Haas, 2005;Gorman et al., 2013). Winter fires consume less vegetation and leave undesirable plants in breeding habitats, which reduces recruitment and ultimately causes localized extinctions and increased fragmentation.
The RFS now only occurs in a fraction of its former range with just six known breeding sites (Farmer et al., 2016;. The RFS breeds in ephemeral, fishless ponds dominated by grasses and forbs (Palis, 1996(Palis, , 1997. These ponds dry during the summer months but fill with rain in the winter (Chandler et al., 2016).
During rain events in October and November, adults migrate to dry pond basins and lay their eggs in anticipation of the seasonal inundation (Brooks et al., 2019b;Palis, 1996Palis, , 1997. Once the basins fill and the eggs are inundated, the larvae hatch, and eventually metamorphose when the ponds begin to dry in the spring (Palis, 1996).
This breeding strategy is highly dependent on predictable seasonal rainfall, and if the ponds do not fill or do not stay full all winter, the entire larval cohort can be lost (Chandler et al., 2016;Palis, 1996).
Breeding sites consist of several separate ponds that are occupied each year and are connected by suitable habitat. Individual salamanders appear to disperse among ponds within a breeding site but do not currently disperse between breeding sites, as estimated from occupancy-based metapopulation models (Brooks, Smith, Frimpong, et al., 2019) and the lack of suitable connective habitat between breeding sites. Estimates of dispersal in relation to geographical distance have led the United States Fish and Wildlife Service (USFWS) to define any occupied pond (and a 460-m radius around it) as a breeding population, and any grouping of ponds within 3.2 km (2 miles) of each other as a breeding site (USFWS, 2015). Populations of RFS are typically managed by USFWS at the breeding site level.
Genetic population structure, effective number of breeders, and gene flow have been previously estimated for RFS with microsatellite data, suggesting that there is some structure among ponds within breeding sites (F ST , 0.004-0.112), limited migration, and a low mean number of breeders per pond (N b : 12.5-30.1, Wendt et al., 2021). Low N E has also been observed in another endangered amphibian, the California tiger salamander (Ambystoma californiense), where N E estimates ranged from 11 to 64 per pond with an average of 30 (Wang et al., 2011). To effectively estimate population structure at small scales, a larger number of genetic markers may be needed since subtle patterns can be obscured when few makers are used (Rittmeyer & Austin, 2015). For example, several hundred single nucleotide polymorphisms (SNPs) were better able to detect fine-scale structure in tiger salamanders (Ambystoma tigrinum) than 12 microsatellites (McCartney-Melstad et al., 2018;Titus et al., 2014). Furthermore, mitochondrial (mtDNA) data can be included in analyses of population structure and gene flow to examine sexbiased dispersal. Mitochondrial DNA is inherited via the matrilineal line, and by comparing estimates obtained for mtDNA and nuclear data, it may be possible to elucidate patterns of migration for males and females.
Understanding patterns of gene flow and population structure at the landscape and local scale is necessary to determine the suitability of translocations for the RFS. The aims of this research were to estimate genetic structure, effective population size, genetic variation, and historic and contemporary gene flow among five extant breeding sites of RFS. We generated potential immune gene SNP data, which were complemented by sequence data generated for two major histocompatibility (MHC) exons and one mitochondrial marker in Williams et al. (2020). Our goal was to leverage these three marker types, at the local and landscape scale, to make estimates between and within breeding sites, and to ultimately explore a more complete image of RFS population structure, genetic diversity, gene flow, and sex-biased dispersal. Immune genes were targeted because pathogens are a major threat to amphibians globally (McCallum, 2007), and variation at these genes is associated with broad resistance to pathogens (Sommer, 2005). These genes may show whether RFS populations have responded to past infections in a way that would create population structure or indicate (through low variation) populations that may be especially vulnerable to disease. Assessing connectivity and genetic differentiation among the remaining breeding sites can inform potential reintroductions and translocations of RFS in order to begin re-establishing this species across its former range and thus reduce its risk of extinction .

| Sample collection
Tissue samples were collected from five RFS breeding sites found on public lands (Figure 1)

| SNP sequencing
DNA was extracted from tissues using a DNEasy Blood and Tissue Kit (Qiagen). Hundreds of RFS immune genes were sequenced in a target enrichment experiment using 19,339 RNA-biotinylated baits designed to target RFS immune genes. To create custom RNA baits, the axolotl (Ambystoma mexicanum) transcriptome (Bryant et al., 2017) was filtered by the gene ontology term "immune response" using the GO2TR 1.0.8 pipeline (Gene Ontology to Target Region, Elbers & Taylor, 2015). This pipeline retained all exons related to immune response (Ortutay & Vihinen, 2006). Using these exons, Arbor Raw sequence reads were demultiplexed using the Illumina BaseSpace platform permitting zero mismatches in the barcodes.
Adapter quality trimming was performed on demultiplexed pairedend reads using all adapters sequences supplied with BBDuk 38.16 (https://sourc eforge.net/proje cts/bbmap/) and the following options: ktrim = right, which trims kmers from 5′ to 3′ direction of a read, k = 23 which sets the kmer length to 23 bp, mink = 11 which sets the minimum kmer matching length to 11, hdist = 1 that uses a hammering distance of one to increase the number of kmers stored, tpe which trims both reads to the same length, tbo which trims adapters based on pair overlap, qtrim = rl for quality trimming right and left, and trimq = 15 which sets the quality trim to 15 using the Phred algorithm. Next, quality and adapter trimmed paired-end reads were mapped to the axolotl 3.0.0 reference genome (https:// axolo tl-omics.org/dl/AmexG_v3.0.0.fa.gz) using BBMap 38.16 (https://sourc eforge.net/proje cts/bbmap/) and the vslow and usejni options. SAM files were converted to BAM files with SAMtools 1.9 (Li et al., 2009). Picard 2.18.10 (http://broad insti tute.github.io/ picard) was used to clean, sort, add read groups, and mark duplicates. The program CallVariants 38.16 (https://sourc eforge.net/ proje cts/bbmap/) was used to call SNPs by ignoring duplicates and keeping SNPs with quality scores greater than or equal to 27. Finally, this dataset was reduced to only di-allelic SNPs with genotype missingness of 5% or less were retained using VCFtools 0.1.15 (Danecek et al., 2011).
Retained SNPs were tested for linkage disequilibrium and Hardy-Weinberg equilibrium (HWE) using VCFtools, geno-chisq, and hardy tests. The p.adjust function in R-4.0.2 (R Core Team, 2018) was then used to correct for multiple tests using the false discovery rate (FDR, Benjamini & Hochberg, 1995). Gene duplication in the RFS, a pattern observed in other amphibians (McCartney-Melstad et al., 2018;Waples, 2015), can lead to heterozygote excess causing This possible gene duplication prompted filtering of the SNP data using stringent conditions. Any SNP that had a FDR less than 0.05 in one or more ponds was discarded from the dataset. Next, BayeScan 2.1 (Foll & Gaggiotti, 2008) was used to predict if SNPs were putatively under selection. As no SNPs were identified as being under selection, SNPs were used to estimate population structure and gene flow under neutral expectations. Finally, blastn v. 2.2.31+ (Altschul et al., 1990) was used to identify bait locations in the axolotl genome, and custom gene annotations were used to predict the genes each SNP might occur in (Sergej Nowoshilow, pers. comm.).
Allelic richness (AR), observed (H O ), expected heterozygosity (H E ), and pairwise F ST (Weir & Cockerham, 1984) were estimated at the pond and breeding site level for SNP data using hierfstat v. 0.04-22 (Goudet, 2005). Heirfstat was also used to conduct principal component analysis (PCA) for SNP data. GenePop v4.6 (Raymond & Rousset, 1995;Rousset, 2008) was used to calculate the inbreeding coefficient for each breeding site (F IS , Weir & Cockerham, 1984) along with an isolation by distance analysis (IBD). IBD was conducted using the decimal degree location for the center of each pond, 10,000 Mantel test permutations, and a minimum distance of 0.001 decimal degrees between samples. Effective population size was calculated for each breeding site using the linkage disequilibrium method in program NeEstimator v2.1 (Do et al., 2014.).
Recent patterns of gene flow were estimated with BayesAss 3.0 (Wilson & Rannala, 2003). BayesAss estimates migration rate over the last three generations by calculating the probability of migrant ancestry for all individuals, assigning them as either a non-migrant, 1 st generation, or migrants that are 2 nd generation or greater. For all analysis, 3 × 10 6 iterations were run starting with a burn-in of 2 × 10 6 followed by 1 × 10 6 iterations with a sample taken every 100 steps for a posterior dataset of 1 × 10 4 samples. These analyses were repeated 10 times using different seeds to ensure convergence of models. STRUCTURE v2.3.4 (Pritchard et al., 2000) and Admixture 1.3.1 (Alexander et al., 2009) were run to assess population admixture using ancestry models to infer the number of populations (K) from the data. In program Admixture, the lowest cross-validation value was used to select K between K = 1 to K = 7. For STRUCTURE, the model that maximizes marginal likelihood was used and results were confirmed by comparing distruct plots (Rosenberg, 2003).

| MHC and mtDNA sequencing
To include other breeding sites and multiple marker types for population structure analyses, we used sequence data generated in Williams et al. (2020) for major histocompatibility complex (MHC) class Iα exon 3 and class IIβ exon 2 and the mitochondrial control region (Williams et al., 2020). An analysis of molecular variance (AMOVA) was conducted for all sampling locations with MHC and mtDNA sequence data using Arlequin 3.5.2.2 (Excoffier et al., 2005).
Samples from Mayhaw and Garcon were removed from migration, STRUCTURE, and IBD analyses because of small sample sizes.
For MHC data, pairwise F ST (Weir & Cockerham, 1984) was calculated in hierfstat v.0.04-22, where MHC sequences were coded as alleles (i.e., each distinct sequence at a locus was given a number, e.g., allele 01). MHC data were only used at the site level because sample sizes at individual ponds were low (≤12) and the SNP data provided a more robust analysis for nuclear data at the smaller spatial scale. For mtDNA, Φ ST was calculated in Arlequin 3.5 using mitochondrial haplotype frequencies. Unlike the MHC data, this analysis was conducted at both the pond and breeding site level.
Principal component analysis, effective population size, and program STRUCTURE v2.3.4 (Pritchard et al., 2000) were run for MHC markers as described for SNP data, but these analyses were not implemented for mitochondrial DNA because they require diploid data.  (Hime et al., 2021;Williams et al., 2013). Consequently, we consider the SNPs generated to be anonymous.

| Genetic variation and effective population size
In the SNP dataset, allelic richness was 1.38 at East Eglin and 1.44 at West Eglin (   (Table 4).

| Population structure
AMOVA results with MHC data indicated little population structure: Only 9.53% of the total variation was among breeding sites (Table 6). Greater population structure was detected with mitochondrial data: 25.49% of the total variation was among breeding sites (

| Gene flow
For program Migrate-N, the mode migration estimate for SNP, MHC, and mitochondrial data was reported, because for SNP and MHC markers, the mean, median, and mode were almost identical (Table S3), but for mitochondrial estimates, the mode better represented the peak of the distribution on the curve (Figures S4-S6).
For mitochondrial data, migration rates were not multiplied by four  Figure S5). For mitochondrial data, models gave wide estimates with long tails. Estimates of migration with mitochondrial data were highest from Escribano to East and West Eglin at 759.0 and 921.0, respectively, while migration out of West Eglin was the lowest at 63.0 to East Eglin and 77.0 to Escribano (Table S3 and Figure S6). Although the lowest confidence intervals for many mitochondrial migration rates were at or near zero, mean, median, and mode estimates were always more than 60 individuals. MHC class Iα estimates were uninformative, probably because range-wide genetic variation was very low.
Migration rates as estimated with SNP data in BayesAss showed an asymmetrical pattern: 32.0% of individuals in the West were migrants from the East while 0.6% of individuals in the East were migrants from the West (Table S4) (Tables 4 and 5).
At a broad geographical scale, the combined data from ponds within each breeding site (i.e., East Eglin, West Eglin, and Escribano) continued to support gene flow and a lack of genetic structure as estimated with SNP and MHC data. However, structure was high between all breeding sites as estimated with mitochondrial DNA (Table 4). The disparity in genetic structure and dispersal estimates for nuclear and mitochondrial DNA is indicative of male-biased dispersal, a pattern observed in other amphibians including the redbacked salamander (Plethodon cinereus) and the alpine salamander (Salamandra atra, Helfer et al., 2012, Liebgold et al., 2011. Females may be capable of moving long distances, but as fall breeders, they rely on cues like previous experience to select suitable egg laying habitat. As a result, they may have a stronger philopatric connection to their natal pond in contrast to males that orient based on receptive females. Thus, females do not appear to disperse as widely given the pressure to select suitable egg laying habitat (Burkhart et al., 2017;Moore & Whiteman, 2016;Peterman et al., 2015).
Although there was some mitochondrial genetic structure over short distances, analyses indicated that Φ ST increased considerably when ponds were separated by more than 1.5 km (Tables 3 and   4), a value similar to the dispersal distances calculated using spatially explicit occupancy models (Brooks, Smith, Frimpong, et al., 2019). This was apparent on Eglin where IBD analysis of mitochondrial data showed no relationship between distance and structure within breeding sites, but did show a significant correlation between East and West Eglin (Tables 3 and 4). In contrast, MHC and SNP data showed no isolation by distance within or between Eglin sites indicating that males may move, or historically moved, between these ponds at a higher rate than females. Notably, although Escribano is treated as one breeding site by the USFWS, mitochondrial and IBD data showed that highly philopatric females may effectively occupy two distinct breeding sites at this location ( TA B L E 6 AMOVA as estimated with MHC and mtDNA data F I G U R E 2 Distruct plots for SNP (a) and MHC class Iα and IIβ (b) data. For K = 2 and 3, individual membership does not align with sampling location indicating a lack of population structure and supporting K = 1. Distruct plots of SNP data broken down by pond found in Figure S3. MHC, major histocompatibility complex; SNP, single nucleotide polymorphisms between sampled ponds have a clumped distribution (currently occupied ponds are either 0.1-1.5 km apart, 2.5-3.5 km apart, or >10.0 km apart) and do not fall on a continuum. Mitochondrial data suggested that females typically do not disperse more than 1.5 km but that may be an underestimate because there were no sampled ponds in the 1.5-2.5 km range.
Other studies have found greater population structure in fall breeding ambystomatids, such as the RFS, than in spring breeders on the same landscape. For example, ringed salamanders (Ambystoma annulatum) and marbled salamanders (Ambystoma opacum), both fall breeders, showed more population structure on Fort Leonard Wood, MO than spotted salamanders, which breed in the spring (Burkhart et al., 2017, Peterman et al., 2015. Fall breeders may show increased philopatry to breeding ponds because they mate when ponds are dry and must deposit eggs in the dry basin. Accordingly, selecting suitable egg laying habitat is more difficult because females must anticipate inundation as compared to the spring where water is present when breeding begins. Consequently, for females, there is a strong selective pressure to choose breeding sites based on cues such as previous experience or pond-associated vegetation (Brooks, Smith, Frimpong, et al., 2019). This pressure to select breeding sites may contribute to greater philopatry in fall breeding salamanders thus decreasing gene flow and increasing genetic structure (Burkhart et al., 2016, Peterman et al., 2015.
We attempted to isolate SNPs from immune genes, and although most were off-target and effectively anonymous, it is possible that selection may play a role in population structure despite our BayeScan results indicating SNP neutrality. However, had there been selection on specific SNPs with different disease challenges at the various sampling sites, we would expect to find population structure rather than its absence. It would seem unlikely that 90 SNPs would all be under selection in the same way at multiple locations so as to cause a consistent lack of population structure across sampling sites.
Future work with a larger number of SNP loci located throughout the genome may help to clarify population structure in RFS.
Measures of SNP genetic diversity were low and mostly uniform  (Wang et al., 2011). Although estimates obtained with microsatellites cannot be directly compared with SNP and MHC markers, it appears that diversity is low across the RFS's sampled range, but N E may be comparable to other endangered ambystomatids.
F IS values were below zero for SNP data ( Alternatively, a more severe bottleneck at one site may have removed genetic variants found at other sites, creating an uneven pattern of diversity that may be interpreted as asymmetrical migration (Wilson & Rannala, 2003 In conclusion, at a broad geographical scale, RFS populations showed little population structure with nuclear SNP data (F ST = 0.00-0.09) but higher structure with mitochondrial DNA (Φ ST = 0.15-0.36). At this same scale, mitochondrial DNA exhibited isolation by distance, whereas nuclear markers did not. These discordant results suggest that males drive dispersal whereas females are more philopatric. Using the distance between sampled sites, and sharp increases in Φ ST , it is likely that females disperse less than 1.5 km from their natal pond (Tables 3-5). Finally, program Migrate-N and BayesAss supported some migration, at least historically, between East Eglin, West Eglin, and Escribano breeding sites.
Future amphibian conservation will require both molecular and traditional conservation work to inform appropriate management strategies, especially in fragmented or shrinking habitats with increasingly less connected populations and fewer active breeding sites. Combining molecular methods with traditional conservation techniques can help identify appropriate source populations for reintroductions, identify sex differences in gene flow, and preserve potentially unique diversity and adaptive potential. Finally, time is a resource that is in short supply: the increasing frequency of strong hurricanes (e.g., Hurricane Michael in 2018 and Hurricane Sally in 2020) as well as the risk posed by global amphibian diseases and inbreeding in small populations underlines the need to begin reestablishing populations in order to decrease the risk of extinction.

ACK N OWLED G M ENTS
We are thankful to Dr.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest for this article.

DATA AVA I L A B I L I T Y S TAT E M E N T
MHC and mtDNA data for this study are available at the NCBI