INTRODUCTION

Rapid environmental change poses a significant concern for the persistence of wildlife populations because of the loss of genetic diversity and adaptive potential (Pauls et al. 2013, Miraldo et al. 2016, Leigh et al. 2019). Levels of genetic variation differ among populations and species and can be influenced by many factors, including historical events, geographic features, and life history and ecological traits (Avise 2000, Leffler et al. 2012, Ellegren and Galtier 2016). These factors may act alone or in combination to generate recognizable patterns of genetic variation both among and within populations. Assessing genetic diversity patterns within co-distributed species sampled across varied landscapes can reveal similarities or differences in these genetic patterns and elucidate the capacity of each species to adapt to changing environments.

Small populations are of concern in conservation biology because they have reduced potential for adaptation (Willi et al. 2006). Under neutral theory, genetic diversity is expected to be positively correlated with effective population size (Ne; Kimura 1979). Populations can experience occasional crashes in population size (i.e. bottlenecks), in which Ne becomes smaller and genetic diversity decreases over a relatively short timescale (Gillespie 2004), or species may persist with naturally small populations for long periods of time (e.g. Stacey and Taper 1992). Small populations can deviate from idealized random mating and may experience a loss of genetic diversity via inbreeding, and also greater risk of genetic drift (Ellstrand andElam 1993, Furlan et al. 2012). However, genetic variation in small populations can be maintained or even increased through migration (i.e. gene flow).

Understanding spatial population genetic structure can lead to conservation decisions that promote maintenance of genetic variation within populations (Chambers 1995, Hohenlohe et al. 2021). Evolutionary processes such as migration and ecological factors such as landscape variation can strongly influence the spatial patterns of genomic variation. A common method to examine factors influencing genetic variation is to test for isolation by distance (IBD; Wright 1943), in which genetic differences between populations increase as geographic distance increases. Whether a species exhibits this pattern depends on its ability to disperse across the landscape, with greater connectivity and gene flow expected to decrease the strength or presence of IBD. Species with lower dispersal ability or species distributed across barriers to gene flow may also exhibit population structure or clustering, which can be inferred using methods to estimate admixture coefficients between individuals (e.g. the sparse non-negative matrix factorization (sNMF) algorithm applied to population genetics; Frichot and François 2015) or through principal components analysis (PCA) of genomic data (e.g. Liu et al. 2020). Identifying spatial patterns within multiple co-distributed species is an important first step to determine the influence of landscape features on population connectivity and the potential for gene flow to maintain genetic variation.

Highly variable landscapes, such as coastal habitats, provide an opportunity to investigate how evolutionary processes shape genetic variation within species and how different species might respond to environmental change. Freshwater organisms in coastal habitats are threatened by urbanization, rising sea levels, and saltwater intrusion (Mitsch andHernandez 2013, Herbert et al. 2015). Amphibians are especially sensitive to increases in salinity because of their highly permeable skin and eggs, along with a limited ability to tolerate hyperosmotic internal conditions (Shoemaker and Nagy 1977, Bentley and Yorio 1979, Lillywhite 2006). However, many more amphibian species than previously recognized can persist within brackish or saline environments and have mechanisms for salt tolerance (Hopkins and Brodie 2015). Investigating the genetic variation among populations within these landscapes is important to document the potential interplay between local adaptation and on-going gene flow.

The Florida Panhandle is one such coastal region of high interest to study genetic variation within and among populations. Located within the North American Coastal Plain, which has one of the highest degrees of amphibian biodiversity in the United States ( Jenkins et al. 2015), this region contains four barrier islands with some species found both on the islands and mainland. The two larger barrier islands, St. Vincent Island and St. George Island, are estimated to have formed approximately 4000 years ago via sediment deposition from the Apalachicola River during sea-level fluctuations (Forrest 2007, López andRink 2007). At present, St. Vincent is approximately 50 km 2 in size, triangular in shape, and 0.5 km from the mainland at its nearest point (Davis and Mokray 2000). In contrast, St. George Island is a long, thin barrier island approximately 73 km 2 in area, 45 km long, < 2 km wide, and < 6 km from the mainland at its nearest point. Whereas St. George Island is largely developed, with a bridge connecting to the mainland, the majority of St. Vincent Island is a National Wildlife Refuge that is only accessible by boat. Non-endemic island taxa tend to have significantly lower genetic diversity than their mainland counterparts, and this pattern has been demonstrated in other barrier islands in Florida (Frankham 1997, Kalkvik et al. 2018). Within the Florida Panhandle, differences in venom protein expression have been found between mainland and island populations of snakes (Margres et al. 2016(Margres et al. , 2017)), but studies of genetic variation in other taxa in this region are lacking. Of the 36 known native amphibian species (20 frogs and 16 salamanders) from the counties adjacent to the barrier islands in this region (Krysko et al. 2011), St. Vincent Island has 11 frogs and one salamander, and St. George Island has eight of those frogs and the same one salamander (the two-toed amphiuma, Amphiuma means; Krysko et al. 2011, iNaturalist 2024).

In this study, we focus on four frog species from three families that have established populations on these barrier islands: the southern toad-Anaxyrus terrestris (Bonnaterre, 1789) (Family Bufonidae), the southern leopard frog-Rana sphenocephala (Cope, 1886) (Family Ranidae), and the green treefrog-Hyla cinerea (Schneider, 1799) and squirrel treefrog-Hyla squirella (Daudin, 1800) (Family Hylidae). Studying four species with differing ecologies and evolutionary histories allows us to compare patterns of genetic variation within coastal and barrier island habitats and determine whether landscape effects are consistent across species. These four species are abundant and widespread throughout the Southeastern Coastal Plain but differ in their genetic divergence patterns on a broad geographic scale (Barrow et al. 2017(Barrow et al. , 2018)). The two species with more shallow genetic divergences, A. terrestris and H. squirella, were predicted to maintain relatively stable ranges since the Last Glacial Maximum (LGM). The other two species, R. sphenocephala and H. cinerea, exhibited deep (> 2 Mya) divergences and were predicted to occur in isolated regions during the LGM (Barrow et al. 2017). These contrasting responses to past climate change may relate to ecological differences; for example, A. terrestris and H. squirella are habitat generalists that may be more resilient to environmental change, whereas R. sphenocephala and H. cinerea tend to rely more on permanent bodies of water and may have more limited dispersal (Conant andCollins 1998, Lannoo 2005). On the other hand, both R. sphenocephala and H. cinerea are known to tolerate brackish wetlands in coastal habitats (Hopkins andBrodie 2015, Albecker andMcCoy 2019), and all four species are commonly found on barrier islands, suggesting they are similarly resilient to such dynamic landscapes. On a fine spatial scale, dispersal ability and patterns of genetic structure in these species are unknown and remain to be investigated.

Here, we sequenced population genomic data from the barrier islands and adjacent mainland of northern Florida for our four focal species. Our primary goal was to investigate the influence of a dynamic coastal landscape on genetic variation in a comparative framework. We addressed three main objectives. First, we compared genomic diversity estimates between mainland and island population within species. Second, we determined spatial patterns of variation within species including the strength of IBD and evidence of genetic clustering. Third, we estimated effective population sizes, inbreeding values, and migration rates between mainland and island populations of each species to better understand the evolutionary mechanisms underlying the observed genetic patterns. We predicted that island populations would have lower genomic diversity compared to mainland populations within species. In addition, we expected to find differences in patterns of spatial genetic variation among species associated with their differing ecologies and dispersal abilities. Within each species, we expected to find some degree of structure between island and mainland Genetic variation of frogs on barrier islands • 3 populations. Alternatively, given the relatively recent formation of the barrier islands and potential on-going migration, there may not be an apparent effect of these landscape features on genetic patterns within species.

MATERIAL AND METHODS

Data collection and assembly

We sampled individuals from two barrier islands (St. George and St. Vincent Islands) and the adjacent mainland in the Gulf Coast region of northern Florida (Table 1; Fig. 1). Frogs were collected by hand according to approved permits from January 2011-July 2014 (Supporting Information, Table S1). Individuals were either salvaged as roadkill, euthanized, and dissected, or toe clipped and released following approved animal care protocols from the Florida State University Animal Care and Use Committee (Protocol Numbers 1017 and 1313). Voucher specimens are archived at the Florida Museum of Natural History and remaining tissues are archived in the Museum of Southwestern Biology Division of Genomic Resources (Supporting Information, Table S1). We sequenced a total of 36 individuals (nine per species), including N = 4 per species from the barrier islands and N = 5 per species from the mainland, spanning a maximum distance of 118 km between localities.

Genomic DNA was extracted from each sample using an E.Z.N.A Tissue DNA Kit (Omega Bio-tek, Inc.) following the manufacturer's protocol. We prepared DNA libraries using a double-digest restriction site-associated DNA (ddRAD) sequencing protocol (Peterson et al. 2012). We digested 250 ng of input DNA with the restriction enzymes SbfI-HF and MspI, while ligating adapters with unique barcode sequences for each sample within an index group. We then pooled samples within index groups and added a second index primer for each group using PCR with 2× Phusion High Fidelity Master Mix (New England Biolabs) and 18 cycles of 98 °C for 5 s, 60 °C for 25 s, and 72 °C for 10 s, followed by a final extension at 72 °C for 5 min. Pooled reactions were cleaned with Sera-Mag Speedbeads (Rohland and Reich 2012), size selected to 300 to 450 bp using a BluePippin system (Sage Science, Inc.) to extract the selected DNA fragments, and quantified them with a Bioanalyzer system (Agilent Technologies). We sequenced libraries on an Illumina HiSeq 4000 with 150 base pair paired-end sequencing at the Ohio State University Comprehensive Cancer Center Genomics Shared Resource.

Sequenced reads were cleaned and assembled using the bioinformatics pipeline Stacks v.2.61 (Catchen et al. 2013). First, we used the 'process_radtags' program to demultiplex reads by individual barcode, removing reads with uncalled bases, trimming reads to a length of 138 bp, and removing low-quality reads (a raw phred score below 10). Second, we used the 'denovo_map. pl' program to assemble the R1 reads and call variants for each species to generate four within-species datasets. We required a minimum of 10 reads to make a stack (-m) and allowed three mismatches between stacks within (-M) and between (-n) individuals. Sequence data were processed using resources from the University of New Mexico Center for Advanced Research Computing.

Genomic diversity within and across species

We used the 'populations' program in Stacks v.2.61 (Catchen et al. 2013) to obtain genomic diversity metrics within individuals and populations. Nuclear diversity from one or few individuals can be used to represent species-level or population-level diversity when thousands of loci are sampled from the genome (e.g. Chen et al. 2017, Grundler et al. 2019). We computed population genetic statistics including nucleotide diversity (π) and expected heterozygosity at the individual-level including both variant and fixed positions. Subsequent analyses were conducted in R (R Core Team 2023).

Within each species, we compared genomic diversity estimates between individuals sampled from the barrier islands (hereafter, 'Island') and individuals sampled from the mainland (hereafter, 'Mainland'). We visualized the differences among species and the variation within species using box plots. Welch Two Sample t-tests assuming unequal variances between groups were performed within each species to test whether there was a significant difference in π between Island and Mainland populations.

Genetic structure and patterns of variation

To examine spatial patterns of genetic variation within species, we first tested for signatures of IBD. For each species, we calculated Euclidean distances between all sampled individuals utilizing the 'distm' function in the geosphere package (Hijmans 2022). We then used the dartR package (Gruber et al. 2018, Mijangos et al. 2022) to calculate Euclidean genetic distances between individuals based on allele frequencies using the 'gl.dist. ind' function. We tested for IBD using the 'mantel.randtest' function of the ade4 package (Thioulouse et al. 2018) by performing We next ran genetic clustering analyses to examine intraspecific population structure and to determine whether Mainland and Island populations for each species were different. We created population structure plots with the R package LEA (Frichot and François 2015), which is designed for population genomics, landscape genomics, and genotype-environment association tests. Using the 'snmf ' function, we estimated ancestry coefficients from the genotypic matrix (indicative of gene flow) and evaluated the number of ancestral populations based on genetic similarity. We conducted runs setting the number of ancestral populations from K = 1 to K = 5, each with 10 replicates, and estimated the cross-entropy criterion to determine the best value of K. For each species, we visualized the inferred ancestry proportions when K = 2 and K = 3.

To further investigate patterns of genetic variation within species, we used the 'gl.pcoa' function in the R package dartR to identify genetic clusters using Pearson Principal Component Analyses (Gruber et al. 2018, Mijangos et al. 2022). Pearson Principal Component Analysis (PCA) is a dimensionality reduction method that takes an input matrix of data and reduces variables into a few principal components (PCs) that maintain most of the variation in the dataset. We used a matrix of single nucleotide polymorphisms (SNPs) as the input data for the analysis and ran each species separately. We then used the function 'gl.pcoa.plot' to visualize the results on two-dimensional plots, examining variation along the first three PC axes.

Effective population size, inbreeding, and migration rate estimates

We estimated the effective population size (Ne) for each species using the Linkage Disequilibrium (LD) method implemented in NeEstimator v.2.1 (Do et al. 2014). We assumed random mating and removed singleton alleles. NeEstimator was run with the 'gl.LDNe' function in the dartR R package (Gruber et al. 2018, Mijangos et al. 2022). Confidence intervals (CI) were calculated with the jack-knife method ( Jones et al. 2016). Within each species, we estimated Ne for all nine samples as a single population, and for datasets with Island (N = 4) and Mainland (N = 4 or 5) individuals considered as separate populations. We tried datasets with and without an individual from the mainland of Franklin County, Florida, sampled from the coast near St. George Island for all species.

We estimated contemporary migration rates between and within the Island and Mainland populations for the four species and the mean posterior estimates of inbreeding coefficients (F) of Island and Mainland populations using BayesAss v.3.0.5 (BA3; Wilson and Rannala 2003). We used PGDSpider v.2.1.1.5 (Lischer and Excoffier 2012) to convert the .vcf file obtained in STACKS from each species to Immanc format (.inp) to run 1 and Supporting Information, Table S1. The map was generated using QGIS v.3.10.4.

BA3. We ran BA3 using the program executable file BA3SNP for use with SNPs, for each species using 40 million MCMC iterations. We discarded the first 4 million iterations as burn-in and sampled every 100 iterations. To check for convergence of each MCMC run, we visualized the BA3 trace results in Tracer v.1.7.2 (Rambaut et al. 2018).

RESULTS

Data summary

Across the 36 individuals, we sequenced an average number of 5.317 (range: 3.036-8.329) million reads. Using STACKS, we assembled an average of 3.827 (1.841-5.615) million reads, with an average coverage of 171 (99-283) reads per locus (Supporting Information, Table S1). Within species, we assembled an average of 0.809 (0.778-0.822) million sites for Anaxyrus terrestris, 1.226 (1.169-1.248) million sites for Rana sphenocephala, 1.304 (1.284-1.318) million sites for Hyla cinerea, and 1.532 (1.336-1.565) million sites for Hyla squirella, from which we estimated genomic diversity metrics. The final datasets used for subsequent population genetic analyses included 5343 SNP loci for A. terrestris, 8266 SNPs for R. sphenocephala, 8202 SNPs for H. cinerea, and 8137 SNPs for H. squirella. The total amount of missing data, calculated as the percentage of missing alleles in the SNP datasets, was 5.525% for A. terrestris, 5.295% for R. sphenocephala, 4.596% for H. cinerea, and 4.237% for H. squirella.

Genomic diversity within and among species

We found that all four species followed the expected pattern of higher genetic diversity for the Mainland population compared to the Island population (Fig. 2). The π values were significantly higher in Mainland individuals compared to Island individuals for A. terrestris (Welch Two Sample t-test with unequal variances: t = -2.879, d.f. = 6.38, P = .026); R. sphenocephala (t = -3.149, d.f. = 4.94, P = .026); H. cinerea (t = -4.061, d.f. = 4.86, P = .010); and H. squirella (t = -2.868, d.f. = 4.84, P = .036).

The average overall π was highest for R. sphenocephala (0.00391), followed by H. cinerea (0.00364), and A. terrestris (0.00330), while H. squirella had a lower average π of 0.00254 (Fig. 2). Average expected heterozygosity followed the same pattern (0.00195 for R. sphenocephala, 0.00182 for H. cinerea, 0.00165 for A. terrestris, and 0.00127 for H. squirella). We observed more variation in π values among Mainland individuals compared to Island individuals for H. cinerea, H. squirella, and R. sphenocephala, whereas it was approximately equal for A. terrestris. The results for A. terrestris also demonstrated more variation in π among Island individuals compared to the other three species (Fig. 2).

Genetic structure and patterns of variation

We found significant signatures of IBD in all species except for A. terrestris (Mantel test statistic r = 0.144, simulated P-value = 0.230). We observed the strongest correlations between genetic and geographic distance in H. squirella (r = 0.713, P-value = 0.001) and R. sphenocephala (r = 0.701, P-value = 0.001), followed by H. cinerea (r = 0.647, P-value = 0.001) (Fig. 3; Supporting Information, Fig. S1).

The best fit number of ancestral populations for each species estimated by sNMF was K = 1 (Supporting Information, Fig. S2). When we visualized the patterns for K = 2, however, Mainland and Island populations for A. terrestris and R. sphenocephala separated into different clusters with some admixture (Fig. 4A, C). The Mainland and Island populations within H. cinerea and H. squirella also separated when K = 2, but a single Mainland individual of both species from coastal Franklin County (Locality 11; Fig. 1) grouped with the Island population (Fig. 4E, G).

The PCA largely demonstrated similar patterns of separation between Mainland and Island populations within species (Fig. 4). The first two PCs explained 31.5% of the variance for A. terrestris, 32.1% for R. sphenocephala, 29% for H. cinerea, and 32.9% for H. squirella. Mainland and Island populations formed Genetic variation of frogs on barrier islands • 7

separate clusters for all four species along the first and second PCs, with A. terrestris and R. sphenocephala having the most distinct clusters. Hyla cinerea and H. squirella also clustered into Mainland and Island groups except for the single individual from each species in coastal Franklin County (Locality 11). This sample clustered more closely with individuals from St. George Island than other Mainland individuals. Within the Island populations, St. George Island and St. Vincent Island individuals clustered separately for all four species. We also plotted the third PC against the first PC and found similar clustering patterns compared to the plots of the first and second PC axes for A. terrestris, R. sphenocephala, and H. squirella (Supporting Information, Fig. S3). For H. cinerea, the coastal Mainland individual no longer clustered with the Island population when plotted along the third PC.

Effective population size, inbreeding, and migration rate estimates

We found contrasting Ne estimates across species, with the largest population size estimated for H. cinerea (Ne = 399. All MCMC analyses in BA3 converged when we visually inspected the trace plots for each parameter and observed high Effective Sample Size (ESS) values. All four species had higher inbreeding values (F values) in the Island compared to Mainland populations (Table 2). We found the highest inbreeding value in A. terrestris (F ISLAND = 0.465) and the lowest inbreeding value in H. squirella (F MAINLAND = 0.186; Table 2). The BA3 analysis identified evidence of recent individual migrants within and between Island and Mainland populations (Table 2). However, for all four species, migration rates between the Mainland and Island populations were considerably smaller (m = 5-6%) compared to migration rates within either the Mainland or Island (m = 94-95%; Table 2). There was no noticeable difference in migration rates depending on whether the direction was to or from the Island population.

DISCUSSION

We generated new genome-wide sequence data for four co-distributed frog species and investigated genetic variation in a dynamic coastal and barrier island landscape of northern Florida. All four species exhibited higher average genomic diversity in the mainland compared to the barrier islands, suggesting a consistent effect of small island area on diversity within multiple species. We also found differences in average genomic diversity among species, with the lowest estimates for Hyla squirella, which corresponded with the lowest effective population size (Ne) estimated for that species. We found significant IBD patterns for all species except Anaxyrus terrestris, and we found contrasting clustering patterns across species. Although K = 1 was the best sNMF model for each species, indicating limited genetic structure, mainland and island individuals were somewhat differentiated in the K = 2 sNMF model and the PCA. All four species also showed evidence of recent migration between island and mainland populations, but with substantially higher migration rates within island and within mainland populations. We discuss these new results in the context of previous work, address potential sampling limitations, and highlight potential avenues for future investigation in this and other coastal-barrier island systems.

All four species in our study showed the predicted pattern of lower genomic diversity within Island populations when compared to their Mainland counterparts. Though this biogeographic pattern has been well documented across other taxa (Cardoso et al. 2009, Hufford, Mazer and Hodges 2014, Lourenço et al. 2018) and island types (Sonsthagen et al. 2012, Altamirano-Ponce et al. 2023), fewer studies have documented this variation across amphibians over such newly formed coastal islands (e.g. Velo-Antón et al. 2012, Duryea et al. 2015). Potential factors contributing to low genetic diversity in Island populations could be related to small effective population sizes (Ne), founder effects, and genetic drift that can act strongly on small and newly arrived populations (Wright 1931, Mayr 1942, Frankham 1997). However, because of the limited number of individuals sampled from each species and population, we were not able to estimate Ne within populations, which would provide insight into potential colonization histories that can cause lower nucleotide diversity of barrier island frogs. The datasets estimating Island and Mainland populations separately resulted in Infinite values which appear when sampling error is larger than expected (Do et al. 2014). Small sample sizes of fewer than five individuals were also unreliable in estimating Ne in other amphibian populations (Trumbo et al. 2023). In our full datasets that estimated Ne for each species, the upper confidence interval was Infinite in three of the four species. Simulations suggest that larger census population sizes result in higher chances of negative/infinite Ne estimations, with at least 1% of the census population needing to be sampled to obtain precise but potentially biased estimates (Marandel et al. 2019). While census population sizes of our frog species in this area are unknown, it is likely that our sample sizes of nine individuals are less than 1% of the census population (Pham et al. 2007). Violation of the assumption of nonoverlapping generations can also lead to bias in Ne estimation (Waples et al. 2014). Increasing sample sizes may allow for more accurate estimates in future studies.

Average genomic diversity differed among the four species in our study, with the highest values for R. sphenocephala and the lowest for H. squirella. Levels of genomic diversity differ among species across the Tree of Life, and this topic has received renewed attention in recent years (Leffler et al. 2012, Ellegren andGaltier 2016). In comparisons of genomic diversity across broad taxonomic groups such as animals (Romiguier et al. 2014) or plants (Chen et al. 2017), key species traits related to reproductive strategy and lifespan have been associated with differing levels of diversity-species with more parental investment, fewer offspring, and longer lifespans were genetically less diverse. Studies conducted within taxonomic groups, where life history variation is often narrower, have found mixed results, ranging from few or no predictors of genomic diversity (e.g. Genetic variation of frogs on barrier islands • 9

Australian lizards- Singhal et al. 2017), or identifying species traits such as body size as negatively correlated with diversity (e.g. butterflies- Mackintosh et al. 2019;mammals-Brüniche-Olsen et al. 2018). Comparisons of genomic diversity within anurans are lacking, but recent studies comparing within-species diversity for > 100 species based on a single mitochondrial gene have found little evidence for a link between diversity and reproductive or ecological traits (Barrow et al. 2021, Amador et al. 2024). Neutral genetic diversity levels are expected to be positively correlated with effective population size, and our results from this initial comparison of four species partially align with this expectation in that H. squirella had both the lowest average genomic diversity and the lowest Ne. Despite having the smallest Ne, H. squirella populations had the lowest inbreeding coefficients of the four species. Lower genomic diversity and increased risk of inbreeding are unlikely to be a concern for any of these species, which are all ranked as 'Least Concern' and are typically abundant and widespread (IUCN 2023). Future studies incorporating more anuran species with genome-scale data will enable more robust tests of whether ecological or life history traits correspond with differences in genomic diversity among species.

Geography can also influence genomic variation through different processes, for example isolation via geographic distance or barriers to dispersal. Our results showed strong signatures of IBD across our barrier island system for three of the four studied species, except for A. terrestris. These results could be explained by the better dispersal abilities of A. terrestris promoting gene flow between populations. Dispersal rates of A. terrestris have previously been estimated over a much broader geographical and temporal scale (Barrow et al. 2017). In that study, dispersal parameters were estimated for the four species in our study based on mitochondrial (mtDNA) phylogenies sampled across the Southeastern Coastal Plain. According to that study, R. sphenocephala had the lowest dispersal rates, while the two treefrogs had intermediate dispersal rates, with H. cinerea having a higher dispersal rate than H. squirella, and A. terrestris having the highest dispersal rate. In a follow-up study, IBD was also tested for the four species on a broad geographic scale across the Southeastern Coastal Plain, and again, these results corresponded to our findings. Anaxyrus terrestris exhibited the weakest genetic structure, which may be attributed to its high dispersal ability and physiological tolerance allowing it to migrate farther from aquatic habitats (Thorson andSvihla 1943, Lemckert 2004). Future work with larger sample sizes and additional populations would be useful to test alternative models such as isolation by resistance (IBR) to predict barriers such as the ocean interfering with dispersal between suitable habitats (Cushman et al. 2006), and isolation by environment (IBE) to correlate environmental and genetic differences spatially (Wang and Bradburd 2014). Testing for all three models (IBD, IBE, and IBR) would provide additional insights into the distribution of habitats within the landscape and how they can limit dispersal and create population structure.

Regardless of the expected recent colonization of these barrier islands, given both their young age and proximity to the coast, we found some evidence of genetic differentiation between Mainland and Island populations. Our sNMF analyses suggested a panmictic population (K = 1) was the best model for all species (Supporting Information, Fig. S2). Panmixia has been proposed previously as the best migration model in this barrier island system for venomous snakes (Margres et al. 2017). This model is expected when differentiation and colonization is recent, the distances between mainland and island habitats are short, or the organisms exhibit some degree of saltwater tolerance (e.g. Zavodna et al. 2005). However, the sNMF K = 2 model and PCA results showed similar clustering patterns that indicated subtle genetic differentiation between the Mainland and Island populations for all four species (Fig. 4). These two methods are functionally different, with sNMF being modelbased and inferring ancestry coefficients based on cross-entropy, while PCA is a dimensionality-reduction method that does not account for ancestry. The PCAs differed primarily from the sNMF results by further separating the Island individuals into two clusters corresponding to each island for all four species, whereas sNMF K = 3 models only showed this pattern for H. cinerea (Supporting Information, Figs S3-S4).

The genetic patterns we observed suggest that the body of saltwater between mainland Florida and the St. George and St. Vincent Islands is acting to some extent as a barrier to dispersal in these species (e.g. Lourenço et al. 2018). However, in the case of H. cinerea and H. squirella, one individual from the coastal Mainland in Franklin County clustered with the Island populations possibly due to the proximity between these populations (see Fig. 1). This pattern was identified both in the PCA and in the sNMF analyses with small admixture proportions (Fig. 4). Previous research on other island systems has shown that saltwater is not an absolute barrier to amphibian dispersal (Duryea et al. 2015), indicating that the continued dispersal for these four frog species is plausible.

We expected that ongoing migration could be high between populations for all species because of the proximity between the islands and mainland Florida, river flooding and rafting, or even human-mediated movement. Interestingly, migration rates between Mainland and Island populations were relatively low for all four frog species (Table 2), which contrasts with the panmictic model suggested by the sNMF analyses. This result might be explained by our low sample sizes that are not capturing accurate migration rates or it could indicate that the low number of migrants is sufficient to connect these populations. Additional individuals and different tools (e.g. Gronau et al. 2011, Pickrell andPritchard 2012) could be used to measure past and current gene flow between these populations in future studies, improving our understanding of the effects of gene flow on frog populations in a barrier island system. For example, it could be important to know whether populations are locally adapted to islands or coastal environments, in which case gene flow from inland populations could be maladaptive. These populations could provide an opportune system for examining local adaptation vs. gene flow (or maladaptive gene flow), where gene flow can reduce local adaptation and prevent population differentiation (e.g. Fitzpatrick and Reid 2019).

Our study examined thousands of loci from a small number of individuals sampled per species, thus our results should be interpreted with this in mind. Despite our sampling effort and approach, we can make comparisons and identify similarities with previous results. The population structure for the four species of interest has previously been studied across a broad geographic range in the Southeastern United States using mtDNA, microsatellites, and SNPs from target capture sequencing (Newman and Rissler 2011, Hether and Hoffman 2012, Barrow et al. 2017, 2018). These studies also found patterns of IBD and population structure in the species of interest, albeit on a much larger geographical scale than ours. Here, we demonstrate that subtle population structure can be detected within these species on a fine spatial scale and provide the first comparison of anuran genomic variation in this barrier island system.

One area of interest for further study is whether selection is acting on coastal and barrier island populations and contributing to differentiation of adaptive traits. For example, coastal populations of H. cinerea that occur in brackish environments exhibited differences in gene expression that suggest local adaptation to salt tolerance (Albecker et al. 2021). Our use of ddRAD enabled cost-effective assessments of genome-wide variation for multiple species, but these approaches only sample a small portion of the genome and are likely to miss selective sweeps (Tiffin and Ross-Ibarra 2014). The frog species we studied currently have no reference genomes, but future whole genome sequencing of individuals from the Island and Mainland populations could be used to identify adaptive loci. Additionally, the use of adaptive loci rather than genome-wide genetic variation for conservation genetic approaches has been debated. Teixeira and Huber (2021) argue that neutral genetic diversity is not a good predictor of extinction risk and that understanding functional genetic diversity is more meaningful. In contrast, Kardos et al. (2021) argue that genome-wide genetic variation plays a crucial role in conserving biodiversity and that focusing on conserving a few functional loci and ignoring genome-wide variation could be detrimental. Whole genome data may be useful for evaluating the roles of both functional loci and genome-wide neutral genetic variation in the conservation of small populations such as those found on islands.

CONCLUSION

Our results show that barrier island populations have lower genetic diversity than their mainland counterparts and this pattern is consistent across species. Island and coastal populations are likely to experience dramatic habitat alterations in the coming years associated with urban development, climate change, and sea level rise (Mitsch andHernandez 2013, Herbert et al. 2015). An improved understanding of the capacity for different species to persist in changing landscapes is essential for efforts to maintain biodiversity. Our study provides an initial picture of the genomic variation within four species that are currently considered widespread and abundant throughout the Southeastern United States (IUCN 2023). In addition to sequencing more individuals per population in the future, resampling these populations over time would provide insight into how genomic diversity and connectivity may change in this dynamic system. Although these species appear resilient to challenging landscapes, such as barrier islands with relatively little freshwater, increased habitat fragmentation and decreased genetic diversity and overall fitness still threaten future population survival. Continued monitoring of common species is critical to help keep them common and provide insight into the evolutionary processes shaping natural populations.