Strong population genetic structure and cryptic diversity in the Florida bonneted bat (Eumops floridanus)

Knowledge of the genetic structure and cryptic diversity is essential for the conservation of endangered species. We conducted a genetic survey of the federally endangered Florida bonneted bat (Eumops floridanus) sampled from its USA range in southern Florida. Florida bonneted bats are primarily found in four regions separated by approximately 100 to 250 km, including three western natural areas: Babcock Webb WMA (BW), Polk County (PC), and Collier County (CC) and one urban population on the east coast, Miami-Dade County (MD). We used 22 microsatellite loci and cytochrome b sequences to assess the extent of connectivity and levels of genetic diversity. Populations were highly differentiated at microsatellite loci (overall FST = 0.178) and model-based and ordination analyses showed that MD was the most distinct among pairwise comparisons. Regional populations were small (Ne < 100) with no evidence of inbreeding. Contemporary migration and historic gene flow suggested that regional populations have not frequently exchanged migrants, and thus the divergence among western regions was likely a result of genetic drift. Significantly, mitochondrial DNA revealed that haplotypes from MD were similar or shared with those recognized as Eumops ferox from Cuba and Jamaica, and divergent (1.5%) from the remainder of bonneted bats in Florida. Our data support the management of each of the four populations as distinct population segments, and that BW, PC and CC combined are on an independent evolutionary trajectory from bats in MD. Bonneted bats in Florida appear to harbor cryptic diversity that will require a reassessment of their taxonomy.


Introduction
A recognized early step in formulating management plans for endangered species is to quantify genetic diversity (Coates et al. 2018). This information allows managers to prioritize populations, such as those with high diversity, unique polymorphisms, or that are demographically uncoupled (Taylor et al. 2020). Additionally, measuring genetic connectivity among populations can have direct conservation and management implications (Lowe and Allendorf 2010) including modeling how species will respond to environmental change. For example, sustained immigration into small populations can increase demographic stability and by extension population persistence (Gonzalez et al. 1998;Hufbauer et al. 2015). In contrast, small, disconnected populations may lose genetic variability via drift, lessening the impact of selection and reducing the ability to adapt (Abbott 2011;Hanski and Mononen 2011). Overall, isolated populations are expected 1 3 to display greater levels of genetic differentiation and population genetic structure.
Animals capable of flight are subject to fewer barriers to dispersal, and by extension typically show less genetic structure, than species that only walk or swim (Medina et al. 2018). Among bat species (order Chiroptera), the only mammalian group with active flight, there is substantial variation in levels of philopatry, with implications for population genetic structure, despite their ability of sustained flight. Dispersal in bats is linked to differences in social structure, mating system, and individual behavioral responses to geographic barriers (reviewed in Moussy et al. 2013). Difficulty in obtaining direct estimates of inter-population dispersal in bats has meant that population genetic studies have been the predominant means of estimating their spatial dynamics (Moussy et al. 2013). In non-migratory species, genetic structure has been explained as an outcome of the extent of gene flow resulting from either contemporary processes such as natal dispersal, mating dispersal or recent colonization, or from historical isolation (i.e., phylogeographic processes; Avise 2000). Distinguishing historical processes (i.e., predating anthropogenic impacts) from more recent processes (i.e., contemporary colonization and founder events) can be difficult because the latter may be strongly influenced by the former (Keyghobadi 2007). However, fast-evolving nuclear markers like microsatellites are expected to reflect the finest-scale contemporary processes whereas mitochondrial sequence data provides a better picture of deeper historical processes (Sunnucks 2000).
The Florida bonneted bat (Eumops floridanus, family Molossidae) is a federally endangered bat that is endemic to southern Florida, USA. Despite having one of the most geographically-restricted ranges of any bat in the United States, and recent efforts to examine occupancy across its range (Bailey et al. 2017a), its distribution within Florida remains poorly understood. Initially known only from Miami-Dade County in southeastern Florida (Barbour 1936), Florida bonneted bats were first observed on the west coast of Florida in 1979 and are currently known to roost in Charlotte, Lee, Collier, and Polk counties (Belwood 1981;Angell and Thompson 2015;Braun de Torrez et al. 2016. The Florida bonneted bats in Miami-Dade County roost primarily in buildings (Gore et al. 2015;Webb et al. 2021), and artificial roost boxes. Florida bonneted bats from the other counties typically roost in groups of < 25 in tree cavities or artificial roost boxes in natural areas, but there has been limited evaluation of urban areas (Braun de Torrez et al. 2016). The largest known population is currently in Babcock-Webb Wildlife Management Area (Charlotte County) where > 300 unique individuals have been PIT tagged since 2014 (Braun de Torrez et al. 2020); however, there is little information about the abundance of Florida bonneted bats outside of Babcock-Webb WMA.
There are fundamental aspects of Florida bonneted bat ecology and evolution that must be determined in order to effectively manage and conserve the species. Namely, there is a lack of information on their dispersal ability, the extent of connectivity between populations, and standing levels of genetic diversity. Bonneted bats can cover large distances when foraging (Ober et al. 2017a;Webb 2018), but the extent to which neighboring populations exchange migrants and associated levels of gene flow remains unknown. Because of their fragmented distribution and presumed small population sizes, Florida bonneted bat are thought to be at an elevated risk of extinction due to inbreeding depression, genetic drift, and stochastic events such as major storms (USFWS 2013). Further complicating conservation and management efforts, the evolutionary relationship between Florida bonneted bat closely related species remains unresolved. Until recently, the Florida bonneted bat was considered a subspecies of E. glaucinus, a widespread species complex ranging across Central and South America and the Caribbean (Eger 1977;Best et al. 1997). The Florida bonneted bat was elevated to species status based on morphological differences (primarily overall size) from Eumops in the Caribbean and Mexico (Timm and Genoways 2004;Wilson and Reeder 2005); however, genetic evidence (mtDNA and AFLPs) has not supported species distinction (i.e., monophyly) of Florida bonneted bats from closely related Eumops in the Caribbean and Mexico (McDonough et al. 2008;Bartlett et al. 2013).
Here we estimate population structure of the Florida bonneted bat in southern Florida. Using a combination of microsatellite markers and mitochondrial sequence data we identify the relative roles of contemporary processes and historical isolation in shaping genetic patterns of the Florida bonneted bat. Our primary objectives were to: (1) describe the genetic structure of regional collections; (2) quantify levels of genetic diversity within and among collections; (3) measure levels of genetic differentiation between sites; (4) test for evidence of migration and estimate gene flow; (5) generate estimates of effective population size and model their stability through time; (6) and test alternative historical population models using an Approximate Bayesian Computational (ABC) approach.

Florida bonneted bat sampling and DNA extraction
We obtained Florida bonneted bat samples from four core regions of the species' range in Southern Florida (see Supplemental File for details). These four regions were: 1) the vicinity of Avon Park Air Force Range, Polk County (PC), 2) Fred C. Babcock/Cecil M. Webb Wildlife Management Area, Charlotte County (BW), 3) Collier County (CC), and 4) areas in Miami-Dade County (MD) (Fig. 1). The first three regions are in central and southwestern Florida (hereafter referred to as "western populations") and represent largely natural habitat types including pine flatwoods, freshwater prairies, cypress swamp, and hardwood hammock. Sampling in MD occurred in the developed environment of the Miami metropolitan area.
DNA was obtained from live bats from oral swabs (Catch-All™ Sample Collection Swab, Epicentre cat. QEC89100), wing punches (4 mm biopsy punch of wing membrane, Worthington Wilmer and Barratt 1996), or both. After optimizing PCR conditions, we did parallel DNA isolations on both swab and wing punch samples collected simultaneously from 20 bonneted bats to test for genotyping consistency across sampling methods. Detailed sampling methods are provided in the Supplemental File.

Primer development and marker assessment
We included one published locus (TabrD15; Russell et al. 2005) and developed microsatellite loci de novo from PacBio sequence reads following the methodology outlined in Saarinen and Austin (2010). Primer information, including details on development, Genbank Accession numbers, and PCR conditions are provided in the Supplemental File.
We performed a suite of tests to ensure the suitability of our microsatellite panel to infer population genetic parameters. First, frequencies of null alleles were estimated within each regional sample and globally using the method of Brookfield (1996) implemented in the R (vers. 3.6.1, R Core Team 2019) program PopGenReport vers. 3.0.4 (Adamack and Gruber 2014). Second, pairwise tests for genotypic linkage disequilibrium (LD) were done using Fisher's exact test implemented in Genepop R v. 1.13 (Rousset 2008). Departures from the null hypothesis of no LD were evaluated using sequential Bonferroni correction at α = 0.05. Third, the power of the 22 loci to detect genetic heterogeneity among regional samples was tested using POWSIM vers. 4.0 (Ryman and Palm 2006). We simulated empirical sample sizes at effective population sizes of 500, 1000, 3000, and 5000, with F ST calculated at generation 0, 20, 35, 50, 75 and 100. Fisher's exact tests were used to test for homogeneity with tests at generation 0 expected to give false significance (α) at approximately 0.05. The proportion of significant tests over 1000 replicates at each value of F ST was used to indicate power of the markers to detect differentiation due to drift. Lastly, we tested for departures from Hardy Weinberg Equilibrium (HWE) within regional samples and globally using the hw.test function in the R program pegas vers. 0.13 (Paradis 2010). Significance was assessed using exact tests (Guo and Thompson 1992) with 1000 permutations.

Assessing genetic structure and population genetic diversity
We used multiple analytical approaches to infer genetic structure among regional collections of samples (Fig. 1). Specifically, we used both a model-based approach implemented in STRU CTU RE vers 2.3.4 (Pritchard et al. 2000) and two non-model multivariate approaches; Principal Component Analysis (PCA) and Discriminant Analysis of Principal Components (DAPC). In STRU CTU RE, we applied the linkage model and evaluated K ranging from 1 to 10 with 20 replicate runs per K, 100,000 burn in iterations and 200,000 sampling iterations. Model (K) selection included 1) Fig. 1 South Florida indicating the regional locations of known Florida bonneted bat breeding locations. In some locations only general coordinates were known and some sample coordinates were jittered to avoid overlap. Inset contains areas being considered for critical habitat designation by the USFWS the maximal estimate of the posterior probability of the data for a given K, Pr(X | K), hereafter L(K); 2) The adhoc ∆K method (Evanno et al. 2005); and 3) The supervised measures (MedMedK, MedMeanK, MaxMedK, and MaxMeanK) proposed by Puechmaille (2016). For the latter, we evaluated inclusion thresholds of 0.5, 0.6, and 0.7 (Puechmaille 2016). Finally, based on the results from L(K), ∆K, and the supervised measures, we inspected the results from a range of K values (Novembre 2016) to evaluate their potential implications for the interpretation of latent genetic structure.
Principal component analysis was conducted using the function dudi.pca in the R package ade4 vers. 12.7-16 (Dray and Dufour 2007) applying weighted vectors and mean centering. We transformed absolute variance of principal component axes to percentage of total variance represented in the genotype data. DAPC was conducted in the R package adegenet vers. 2.1.3 (Jombart 2008; Jombart and Ahmed 2011), following Miller et al. (2020) (Supplemental File). DAPC reduces intra-group principal component variation to maximally explain between group variation (Jombart et al. 2010). We first performed crossvalidation for discriminant analysis using the xvalDapc function. We tested up to a maximum of 200 PCA components, retaining three axes for each DA (one less than the number of a priori regional groups). DA clusters were plotted in adegenet with an overlayed minimal spanning network, implemented in the R package ade4, and we examined posterior assignment values for individual bats to assess structure.
Once the most likely number of genetic clusters (above) were identified, we then described genetic variation within and among clusters. Observed (H o ) and expected heterozygosity (H e ), unbiased heterozygosity (uH e ), inbreeding coefficient (F IS ), number of alleles per locus (A) and the 95% confidence interval (CI) based on 999 bootstrap intervals were estimated using the R program diveRsity vers. 1.9.90 (Keenan et al. 2013). Allelic richness (A r ) was calculated using rarefaction to the smallest regional sample size (PC, n = 15). We evaluated the proportion of total genetic variation that was contained within and among regional populations using analysis of molecular variance (AMOVA) implemented in GENODIVE vers 3.04 (Meirmans and Van Tienderen 2004). We performed a standard AMOVA to assess the proportion of variance among the four sampling regions and a second AMOVA in which western populations (BW, CC, PC) were grouped together to evaluate the variance between MD and western locations (assuming hierarchical genetic structure, see "Results" section).
Levels of genetic differentiation between clusters were estimated via F ST (Nei 1987) with 95% CI calculated using 9,999 bootstraps replicates. We calculated pairwise F-statistics and D Jost using the function fastDivPart in PopGen-Report (Adamack and Gruber 2014). Estimates of fixation (and gene flow) may be biased when subpopulations have high levels of unique allelic variation (i.e., display negative dependence; Jost 2008). To test for negative dependence, we used the corPlot function in diveRsity to plot F ST (Θ)/D Jost calculated for individual loci against the number of alleles per locus. Contrasting slopes between estimators (e.g., negative F ST and positive D Jost ) was taken as evidence of bias stemming from elevated mutation rates (Keenan et al. 2013), while similar slopes (e.g., all positive) suggest that drift, not mutation, was primarily shaping measures of fixation (F ST ) and differentiation (D Jost ).
We also tested for isolation by distance among regional populations to assess whether diagnosed structure could reflect a gradient of differentiation. We regressed (F ST / (1 − F ST )) on log geographic distance (centroid of sampling location within each region) using 1000 bootstrap permutations to estimate the 95% CI of the regression slope in Genepop. We used a Mantel test to examine for significant correlation of genetic (Euclidian) and geographic distance matrices, and we ran a stratified Mantel test permuting individuals within BW + CC + PC to allow for the potential influence of hierarchical structure as inferred from PCA, DAPC and DIYABC (see below) results. Both Mantel tests utilized 9,999 permutations and were conducted in GENODIVE.

Estimates of migration and gene flow
To test for recent migration among populations we used the GENSBACK option in STRU CTU RE to evaluate the probability of recent migration or mixed ancestry. This analysis used USEPOPINFO = 1 to identify which individuals were unlikely to be residents of their sampled population or were a descendant of a recent migrant. STRU CTU RE was run for K = 4 based on outputs from previous STRU CTU RE and multivariate analyses (see "Results" section). We set GENBACK = 2 to assign migrants (0) and recent migrant ancestry (parental = 1, grandparent = 2). We evaluated different migration rate priors (MIGPRIOR) to assess sensitivity to this variable as an initial condition (Pritchard and Wen 2004), using a range of values (0.01, 0.05, and 0.1) and ran 200,000 and 500,000 iterations for burnin and run length, respectively.
We further quantified recent migration using Bayes-Ass ver. 3.0 (Wilson and Rannala 2003). BayesAss estimates the posterior mean migration rates based on posterior individual ancestry. We used a total of 1.0 × 10 7 iterations, discarding the first 1.0 × 10 6 as burnin and a sampling frequency of 100 to avoid autocorrelation between samples. Final parameters were dF = 0.3, dM = 0.1, dA = 0.3. Migration rate estimates were considered non-zero where the 95% credibility set (mean ± 1.96 × stdev) did not overlap zero.
While the previous approaches can identify contemporary (0-2 generations) migration, we also estimated the directional relative migration among populations and the symmetry of gene flow using the divMigrate function in diveRsity (Alcala et al. 2014;Sundqvist et al. 2016). We estimated effective number of migrants (Nm) to quantify gene flow and assessed the asymmetry of gene flow among pairs of populations by looking for overlap of 95% CI around estimates of migration generated using 1000 bootstrap values (Sundqvist et al. 2016). Overlapping confidence intervals were interpreted as evidence for symmetrical gene flow.

Demographic analyses
To inform conservation and planning efforts we generated multiple N e estimates: (1) contemporary (i.e., previous 1-2 generations) N e based on patterns of LD; and (2) historical (i.e., over the past 500 generations) changes in N e based on microsatellite repeat motif frequencies and differences between alleles. Contemporary N e was estimated using the LD method (Waples and Do 2008) implemented in NeEstimator vers. 2 (Do et al. 2014). We calculated N e under both random mating and monogamy models as mating system can have a large impact of inferences of N e (Weir and Hill 1980). We used VarEff vers. 1.0 (Nikolic and Chevalet 2014) to test for changes in N e over the past 500 generations, encompassing the period of major anthropogenic change in south Florida. For the final analyses we applied a generic mutation rate prior of 0.001 (Ellegren 2004) and used a two-phase mutation model ('T 0.15'). Additional population-specific parameters are in Supplemental File (Table S2). Following 10,000 burnin iterations we ran 10,000 steps with 10 steps per batch and 100 steps between sampled steps to avoid autocorrelation for a total of 10,000,000 MCMC iterations.
Lastly, we used Approximate Bayesian Computation implemented in DIYABC (Cornuet et al. 2014) to assess the relative support for competing models of population history. We developed four models based on outputs from genetic clustering and tests of genetic differentiation, with competing models differing in patterns and times since divergence from a common ancestral population (Fig. 2). The four models were: (M1) Simultaneous radiation of four populations from an ancestral population at time t A in the past; (M2) Ancestral divergence leading to MD and an unsampled ancestral population N W , with the latter subsequently differentiating into three populations (BW, PC, CC); (M3) an ancestral population diverged into MD and BW populations with subsequent colonization of the CC population from BW at t 2 , followed by the colonization of PC from BW at t 1 ; (M4) an ancestral population diverged forming MD and BW followed by the serial colonization of CC at t 2 and PC at t 1 . Prior selection, summary statistics and assessment of confidence in model choice are provided in the Supplemental File.

Mitochondrial sequencing and analysis
Given the relatively high microsatellite differentiation between MD and western populations (see "Results" section), we generated sequence data for the mitochondrial cytochrome b (cyt b) gene to improve the resolution of deeper historical patterns of divergence. Cyt b performs well at species level delimitation in bats including Eumops (Caraballo et al. 2020). We sequenced a subset of 49 samples: 13 MD, 11 BW, 17 CC, and 8 PC. Details on PCR and sequencing methods are provided in the Supplemental File. Sequences were aligned using Clustal Omega (https:// www. ebi. ac. uk/ Tools/ msa/ clust alo/) then imported into Mesquite vers. 3.7 (Maddison and Maddison 2021) for manual error correction, trimming, and to visually examine codon polymorphisms and confirm the absence of stop codons.
Due to the extent of haplotype divergence detected within Florida bonneted bats, and to improve our phylogeographic inference, we added available cyt b sequences of closely related Eumops lineages representing Cuba (N = 14), Jamaica (N = 3), and Mexico (N = 13) that have been shown to form a well-supported clade that includes E. floridanus (McDonough et al. 2008;Bartlett et al. 2013). We also included three previously sequenced Florida bonneted bats from near Fort Myers, Lee County FL, and from Miami-Dade County, and E. glaucinus (N = 6) from South America (See Supplemental File Table S4).
The number of unique haplotypes (h), polymorphic nucleotide sites (S), haplotype diversity (the probability that two randomly selected haplotypes will be identical, Hd; Nei 1987), and nucleotide diversity (the average number of nucleotide differences per site, π; Nei 1987) were calculated for each location and overall with DnaSP vers. 6.12.03 (Rozas et al. 2017). Pairwise differentiation was calculated as G ST (Nei 1973) and the pairwise average number of nucleotide differences among populations. Evolutionary relationships among haplotypes were visualized using the statistical parsimony network (Templeton et al. 1992;Clement et al. 2000) implemented in PopART vers. 1.7 (Leigh and Bryant 2015) which estimates the most parsimonious pairwise connections among haplotypes at ≥ 95% probability (Templeton et al. 1987). In addition, distance networks were created using the NeighborNet algorithm (Bryant and Moulton 2004) in SplitsTree vers. 4.11.3 (Huson and Bryant 2006). We used an uncorrected p-distance as input to identify similar haplotypes. The algorithm generates clusters of haplotypes and illustrates ambiguity of connections ('nets') where necessary.
Because we identified divergent haplotypes between western and MD bats (see "Results" section), we used DIYABC to test support for two competing models that contrast the matrilineal history of MD bats relative to the western populations of E. floridanus and that of the Caribbean bats. The first model (M t 1 ) assumed that the diversity of haplotypes and differentiation between MD and the other Florida populations reflected the retention of ancestral polymorphism from a Pleistocene bonneted bat population in Florida. The second model (M t 2 ) assumed that western E. floridanus had an independent population history from MD, and that the latter derived separately and more recently from a Caribbean ancestor ( Fig. 2; see Supplemental File for modeling details).

Sampling and marker assessment
We genotyped 125 Florida bonneted bats for our final dataset (BW = 60, PC = 15, CC = 31, MD = 19). Sample sizes were skewed primarily because of the ready access to breeding colonies in artificial roosts at BW, whereas samples from the other locations were from free-flying bats captured in nets or otherwise opportunistically sampled.
There was no evidence of null alleles present within any of the four regional samples. Nine of 22 loci were significant when all genotypes were pooled as a single population (L5989, L55690, TabrD15, L40155, L8400, L109865, L30125, L18018, L32421). Because null alleles were not detected in individual regional samples, and results of analyses of genetic structuring (see "Results" section) supported the independence of regional gene pools, we did not correct for null alleles. The number of alleles per locus ranged from 2 to 10.
Fisher's method for LD between locus pairs across all regional samples combined identified four pairs that were significant following sequential Bonferroni correction (all p adj < 0.0002; L55690 and L30125; L5989 and L117077; L38728 and L8400; L11764 and L39543). However, LD tests within individual regional samples revealed that significant values were each driven by a single significant LD test (i.e., 1 of 4 tests for each pair were less than p < 0.05); three significant tests were from PC and one from BW. Because LD was not present between loci in more than one regional sample, we treated loci as independent.
POWSIM simulations identified high power for our marker set to detect differentiation due to drift (i.e., in the absence of gene flow) after only a small number of generations. For even very large populations (N e = 5,000) power was above 0.8 after 50 generations. At N e = 500 power was near 1.0 after 10 generations of drift (see Supplemental File). These results suggested that the 22 loci had sufficient power to detect recent genetic divergence due to drift alone.
Two of 88 tests for departures from HWE were significant at α = 0.05 (PC, L91027, p = 0.043; CC, L11764, p = 0.037), but neither remained significant after Bonferroni correction of alpha (adjusted smallest α value = 0.0005). Global departures from HWE were significant (≤ 0.05) at 9 of 22 loci, with five loci remaining significant after Bonferroni correction (not shown). Evidence of population substructure (see below) and limited support for null alleles in conjunction with HWE test results suggested global departures from HWE were likely a product of population substructure (Crow and Kimura 1970).

Population structure and diversity
The number of genetic clusters detected via STRU CTU RE varied depending on the criterion used. For example, the Puechmaille (2016) metrics supported K = 4 across inclusion thresholds of 0.5, 0.6 and 0.7. At K = 4 each regional sample was clustered independently from one another (Fig. 3). Examining models where K = 5 to K = 10, MD, CC, and PC represented largely uniform, discrete clusters, whereas the remaining additional diagnosed clusters were admixed to various degrees among BW bats (Fig. 3). In contrast, the highest L(K) was observed for K = 8, whereas ∆K identified the greatest rate of change moving from K = 5 to K = 6 (Supplemental File, Fig. S1); however, ∆K values were small overall (all < 20; Supplemental File, Fig. S1).
Results from multivariate analysis mirrored those of STRU CTU RE. Specifically, each regional sample represented a largely non-overlapping scatter in our PCA, with the 1 st axis (12.34% of total variation) separating MD (east) from the western (BW + CC + PC) regions, and the 2 nd axis (7.14%) differentiating CC from northern sites (Fig. 4A). This trend was similar when comparing the 1 st and 3 rd (6.54%) principal components (Fig. 4B).
The results for DAPC resembled results from STRU CTU RE in that increases in K beyond 4 clusters added structure to bats from BW only (i.e., BW bats were assigned to > 1 cluster), whereas the group membership of bats from PC, CC and MD remained the same (not shown). The lowest BIC score (203.02) was at K = 6, followed by K = 5 (203.20), and K = 4 (204.44) (Supplemental File, Fig. S2). We present results from K = 4 for simplicity (Fig. 4C). At K = 4 (and K = 5 to K = 7), one bat genotyped from CC had a higher posterior probability of assignment to the BW cluster (p = 0.98), and all other bats were assigned to their regional . Results from multiple K models are shown to illustrate the three regional samples (PC, CC, MD) remained relatively homogeneous in their cluster membership, while additional clusters are apportioned among BW bats population of origin (all p > 0.99). Bats from PC and BW were least dissimilar, while PC + BW were most dissimilar to MD, with CC intermediate (Fig. 4C).
Genetic diversity (Table 1) measured as mean number of alleles ranged from 3.0 to 4.2 among regional populations; corrected for sample size this ranged from 3 to 3.8. Private alleles were found in all four populations: 3 in PC, 9 in CC, 10 in MD, and 12 in BW. Fixation indices (F IS ) were generally low across loci within each population. Locusspecific heterozygosity varied by marker (low: H e 0-0.10 to high: H e > 0.7), but was largely consistent across populations (Supplemental File, Table S6).
Results from the standard AMOVA found that most of the variation in microsatellite allele frequencies was detected within individuals. However, variation among individuals within populations was substantially lower (and  Table S4). Over all loci, 82% of variation in allele frequencies was observed within individuals (F IT = 0.18), 1.7% among individuals within regional populations (F IS = 0.02), and 16.3% among regional populations (F ST = 0.16). With western populations grouped to reflect potential hierarchical structure, 76% of variation in allele frequencies was within individuals (F IT = 0.239), 1.5% among individuals within regional populations (F IS = 0.02), 10.8% among regional populations nested within western and eastern (MD) groups (F SC = 0.12), and 11.6% representing western vs. eastern groups (F CT = 0.116). Regional genetic clusters were dissimilar based on pairwise F ST (range: 0.095 to 0.226), and differentiation based on D Jost (range: 0.08 to 0.218) ( Table 2). The highest levels of differentiation were all pairwise comparisons involving MD. Global values were significantly greater than zero for both F ST (0.178; 95% 0.166-0.192) and D Jost (0.270; 95% CI 0.253-0.288). Twenty-one loci had bootstrap values indicating 95% CI greater than zero for F ST , and 22 loci were significant for D Jost (Supplemental File, Table S7). From the corPlot examination of polymorphism bias, we found that both estimators (F ST and D Jost ) were positively correlated with the number of alleles at individual loci (Supplemental File, Fig. S3), suggesting that the mutation rate at these loci should not obscure their usefulness for inferring past demographic processes such as gene flow or population size.
Results from regression between linearized genetic distance and log geographic distance rejected the null  Table 2 Upper table: pairwise F ST (below diagonal) and D Jost (above diagonal) among regional genetic clusters of Florida bonneted bat based on 22 microsatellite loci (see Table 1 for a description of abbreviations) All values were significantly greater than zero (all P < 0.001). Lower table: Directional gene flow (Nm) estimates between regional samples. Rows reflect directional migration from the row into the column regional sample. Bold indicates Nm is significantly different from the alternative direction of migration hypothesis of no effect of distance in explaining differentiation among identified regional clusters (b = 0.143, 95% CI 0.085-0.281). Individual-level Mantel test also suggested that inter-individual Euclidean distances are positively correlated with geographic distance (r = 0.623). Results from regular and stratified permutations (MD vs western regional clusters) were both highly significant (p < 0.001).

Migration and gene flow
Tests to detect recent migrant ancestry using the GENS-BACK option in STRU CTU RE varied minorly depending on migration prior; all but two individuals were assigned to their sampling location with a probability > 0.9 (not shown). All recent migration rates estimated from BayesAss overlapped zero suggesting the recent migration was not significant (Supplemental File, Table S8). Longer-term gene flow calculations using divMigrate in diveRsity were lower than 1 Nm per generation and symmetrical for all estimates other than for the estimate from BW into CC (Nm = 1.0) which was significantly greater than Nm from CC into BW (Nm = 0.747; Table 2).

Demographic analyses
Assuming random mating, mean estimates of contemporary N e within regional samples were less than 100, with the exception of CC, which was above 100 only when rare alleles were considered (Table 3). Overall, N e estimates were higher when monogamy was assumed (Table 3), and under this scenario MD had an effective population size > 100. Confidence intervals on the N e estimates were narrow and similar between parametric and jackknife methods for BW and PC. For the parametric and jackknife CIs of MD, and for the jackknife CIs for CC, the upper 95% were infinity, suggesting a lack of power (i.e. small sample sizes) to estimate some CI with confidence. Tests for temporal changes in historic N e (i.e., over the last 500 generations) reconstructed using VarEff identified declines in each of the four populations over the past 100-200 generations (Supplemental File, Fig. S4). Declines were most pronounced in CC which showed a dramatic and consistent decline in N e over the last 500 generations, whereas MD and BW displayed an initial increase over the same period before declining.
For the DIYABC simulations, model choice and demographic estimates (ancestral N e and t) were consistent between models based on equal sexes and female bias ratios (median estimates had overlapping 95% CI for all parameters, so we report on only the equal sex ratios here). The ABC modeling of alternative demographic scenarios supported M2 (hierarchical structure with western populations sharing a most recent common ancestor, Fig. 2) based on the weighted logistic regression (P = 0.913, 95% CI 0.907-0.919). The second most supported model (M1, simultaneous divergence from an ancestral population) had posterior support of 0.075 (95% CI 0.069-0.081). The remaining two models had posterior probabilities and 95% CI below 0.01 (Supplemental File, Fig S5). Posterior-based error rate was 0.283. Prior Type I error (probability of rejection assuming hypothesis is true) for M2 was 0.079, and Type II error for M2 was 0.046, suggesting a small likelihood for selection of M2 over M1 assuming M2 is not the true best scenario. Although our primary objective for DIYABC was to assess the support for alternative demographic models, rather than to infer ancestral N e and t, we include these estimates in Supplemental File, Table S9.

Mitochondrial sequencing results
We aligned 934 bp of the cyt b gene (Genbank accession # OK165510-OK165558) for the 49 Florida sequenced samples in this study. There was a total of 16 polymorphic sites representing 7 haplotypes. Haplotype diversity was 0.740 overall and nucleotide diversity was 0.0057. Among all 49 sequences the mean number of nucleotide differences was 5.363, primarily reflecting comparisons between MD and the other sequences. Average nucleotide differences and pairwise G ST were highest between MD and the other three Florida populations (Table 4). The inclusion of the 39 published cyt b sequences reduced the alignment (n = 75) to 705 bp with 72 polymorphic sites for 31 total unique haplotypes. One haplotype from Cuba was shared with MD (haplotype 6), and all MD haplotypes were nested among Cuba and Jamaica haplotypes (Fig. 5). With the addition of previously published sequences there were a total of 19 private missense mutations amongst 14 haplotypes (Supplemental File, Table S8). The remaining missense mutations represented two codon changes (autapomorphies) that distinguished the 38 Florida bonneted bats collected in BW (and Lee Co.), PC, and CC (haplotypes 1-4, and 23-25) from all remaining haplotypes including all from MD (Supplemental File, Table S10).
For the combined Florida bonneted bat and Caribbean/ South American data set Hd was 0.911 overall and nucleotide diversity was 0.0154. MD was more similar (G ST ) to samples from Cuba (0.145), Jamaica (0.177), and Mexico (0.110) than to other Florida populations (BW = 0.247, PC = 0.295, CC = 0.292) ( Table 4).
Statistical parsimony analyses revealed that the MD haplotypes were closely related to the haplotypes from the Caribbean differing by 1 to 4 mutations and were nested within a cluster of Caribbean haplotypes (Fig. 5). In contrast, the haplotypes from BW, PC, and CC were relatively divergent from MD, differing by at least 8 mutational steps (Fig. 5) or 11 steps for the 934 bp alignment (not shown). The two previously published haplotypes from near BW (Lee Co.) were also divergent from MD (a minimum of 8 mutational steps) (Fig. 5).
The DIYABC modeling of the two scenarios, retention of ancestral polymorphism within Florida (M t 1 ) versus secondary colonization of MD from an ancestral population shared with the Caribbean (M t 2 ), showed strong support for secondary colonization (M t 1 , P = 0.017, 95% CI 0.015-0.019 vs. M t 2 , P = 0.983, 95% CI 0.981-0.985). Posterior-based error rate was 0.026 and both prior type I and type II error were 0.

Discussion
We present the first assessment of population genetic structure in the federally endangered Florida bonneted bat. We found strong population genetic structure among all regional sample locations with little evidence of recent connectivity Table 4 The average number of nucleotide differences (above diagonal) and Pairwise G ST (below diagonal), for mtDNA sequence data (934 bp) collected from Florida bonneted bats (see Table 1

Population structure and connectivity within Florida
Non-migratory bats typically display genetic structuring at smaller geographic scales relative to migratory species (Moussy et al. 2013), although patterns of structure vary among non-migratory types (Ruedi and McCracken 2009). We found significant and strong genetic fixation overall (i.e., global F ST = 0.166) and significant pairwise genetic among haplotypes. Vertical hatch lines represent mutational differences separating haplotypes. Numbers for both panels are haplotype identifiers. Note that haplotypes 5 and 6 are from Lee Co. FL, near BW differentiation between populations separated by < 250 km, which was surprising given the seeming absence of major biogeographic barriers. Our observations contrast with patterns from other bat genera (e.g., Myotis) which display non-significant levels of differentiation at comparable scales (Castella et al. 2001;Atterby et al. 2010;Laine et al. 2013;Miller-Butterworth et al. 2014;Tournayre et al. 2019). Higher levels of differentiation (F ST > 0.1) as was observed in Florida are more commonly associated with major biogeographic barriers (e.g., English Channel, Strait of Gibraltar; Castella et al. 2001;Tournayre et al. 2019); however, we are unaware of any fine-scale genetic studies of closely related species (e.g., Eumops spp.) for comparison. Tests of genetic differentiation among regional populations of Florida bonneted bats implies strong site fidelity despite long-distance foraging behavior (> 15 km) and large average home range sizes (50-200 km 2 ) (Webb 2018). The observed levels of divergence support the conclusion that these four regions represent individual demographic units and should be treated as 'distinct population segments' or 'management units' (Funk et al. 2012). Therefore, it will be important to account for the status and stability of all the regional populations when considering habitat conservation decisions and not just the overall status of the species throughout Florida. In addition, maintaining habitat connectivity across large scales will be important to ensure new colony formation is not impeded over time due to a lack of available habitat. The observed levels of genetic structure could be explained by a combination of the polygynous social system and demographic patterns (Ober et al. 2017b;Braun de Torrez et al. 2020). Within a harem-based social structure, females are unlikely to disperse long distances from their natal colony and dispersing juveniles are less likely to secure mates if they travel far from existing colonies. Both scenarios would limit gene flow and genetic connectivity between populations. In Florida bonneted bats, apparent survival of juveniles is low (Bailey et al. 2017b) implying that juveniles either have low survival or tend to disperse. Juveniles that disperse may either fail to establish new breeding colonies or, when successful, do so over small distances relative to the total range of the species.
Small estimates of N e , strong pairwise differentiation, large numbers of private alleles within each population, and short geographic distances separating populations support the interpretation that genetic drift has played a large role in shaping the dynamics of Florida bonneted bats. That we observed evidence of genetic drift and failed to detect evidence of dispersal supports the idea that western populations of Florida bonneted bats may represent small, recently (e.g., past few decades) established populations, as opposed to a scenario of long-term (e.g., > 100 s of years) allopatric isolation. Furthermore, in systems where genetic drift is prominent, spurious relationships between measures of differentiation and geographic distance can be erroneously attributed to processes such as isolation by distance (Prunier et al. 2017). Why regional populations are separated by significant gaps is unknown and requires further investigation but could in part be related to habitat availability/suitability (lack of roosting opportunities), the general spatial stochasticity of species with female defense or resource defense polygynous mating systems (Storz 1999;Parreira and Chikhi 2015) or some combination thereof. Importantly, our results suggest that genetic drift over short time periods (10 s of generations) are likely driving the level of differentiation observed among the western populations (BW, CC, PC), whereas historical factors must be accounted for in explaining the divergence between MD and western populations (see below).
Diversity did not vary strongly among regional populations (Table 4) and overall, each region appeared to fit closely to an equilibrium model of genetic variation. These observations suggest that regional populations are demographically stable, and despite their small effective population sizes, have not undergone significant inbreeding. Polygynous mating systems are expected to reduce N e relative to monogamy as fewer males contribute genes to subsequent generations (Storz et al. 2001). Consequently, we might expect our estimates of N e based on a random mating model to be artificially inflated given that Florida bonneted bats display a polygynous mating system. However, the effect of a polygynous mating system may be counteracted by the suspected life span of this species (~ 5 years, unpublished data), which would allow for more reproductive opportunities over multiple breeding seasons and would increase N e (Clutton-Brock 1988). This is in line with predictions that in highly social species, genotypic diversity can be maintained even in small populations (Parreira and Chikhi 2015). Data from BW suggests that turnover in dominant males within individual roosts and competition from outside males may be common (Ober et al. 2017b;Braun de Torrez et al. 2020), and the potential for cuckolding would further counter the downward impact of polygyny on N e. Taken together, our numbers suggest that low N e is not simply a result of the mating system but may reflect true small census population sizes.

Phylogeographic relationships within Florida and the Caribbean
Eumops in Florida have been proposed to be derived from a Caribbean origin (Baker and Genoways 1978;Morgan 1985) perhaps during one of the Pleistocene glacial stages when Florida's peninsular land mass would have increased by 25% to over 100% (Emslie 1998). Even at its peak gain in land mass, this still would have required colonization across the Florida Strait. Alternatively, range expansion from Mexico along the northern Gulf of Mexico, which facilitated the expansion of some other Neotropical species (Emslie 1998) has been proposed (Morgan 1991). The pattern of variation found in cyt b suggests a Caribbean origin, as western Florida haplotypes are more similar to Caribbean (and MD) haplotypes than they are to those identified from Mexico (Fig. 5).
Florida bonneted bats from MD and from western Florida have been demographically isolated for a substantial period. The presence of deep mtDNA divergence within Florida bonneted bats was initially suspected after close examination of the results presented in McDonough et al. (2008), where two bonneted bats from the west coast (Lee Co. < 20 km south of BW) were represented by a relatively long branch with high posterior support within a larger clade representing Eumops haplotypes from Cuba, Jamaica, Mexico, and one sample from Miami-Dade County (see Fig. 1 of McDonough et al. 2008). The level of divergence between western and MD+ Caribbean haplotypes (~ 1.5%) suggests that western bonneted bats may be descendants of populations represented by Rancholabrean fossils identified from eastern Florida (Nabholz et al. 2008). The similarity in cyt b between MD and Caribbean bats does not support the notion that the current MD population has had a long-term (i.e., early Holocene or older) presence in Florida. The latter scenario (longterm presence of MD haplotypes in Florida) would require incomplete lineage sorting, assuming all Florida bonneted bats share a most recent common ancestor, or alternatively, would suggest periodic maternal gene flow from the Caribbean to MD (mitochondrial introgression). Our ABC model comparison strongly supports the secondary colonization hypothesis over the retention of ancestral polymorphism in Florida (although we did not assess the possibility of gene flow). The addition of genomic data will improve our ability to resolve whether MD bats were recently introduced (Barbour 1936) or were more distantly established (Allen 1932;Koopman 1971).

Management implications
Given that MD and the western Florida bonneted bats likely have had mutually exclusive histories over the Holocene/ Pleistocene, they should be recognized as two Evolutionary Significant Units (ESUs). ESUs are defined here as populations that demonstrate independent evolutionary trajectories and restricted gene flow (Fraser and Bernatchez 2001). Ideally, diagnosis of ESUs would include both neutral and adaptive genetic loci (Funk et al. 2012); however in the absence of genomic data, combined mtDNA and nuclear microsatellites have been used to support the delineation of ESUs in various taxa (e.g., King et al. 2006;Austin et al. 2011;Carvalho et al. 2012;Puckett et al. 2021). While it is unknown whether the identified missense mutations have fitness effects, they do represent rare intraspecific variants that add support to the argument for managing as separate evolutionary units . Separate ESUs, for example, will each require sufficient critical habitat to ensure long term demographic resilience. These units will likely require ESU-specific conservation actions in anticipation of likely environmental shifts coming because of climate change. Management strategies, such as translocations of bats (i.e., remediation translocations) between ESUs would not be appropriate due to their genetic differentiation and potential ecological or adaptive differences. Additional research into the roosting and foraging requirements of bats within each ESU would provide insight into the degree of adaptive divergence between the two ESUs. Similarly, the three western distinct population segments may have different habitat preferences and responses to management actions. Additional research is needed to make appropriate population-level conservation recommendations across the species range and to facilitate movement among populations.
Although examining the taxonomic status of the bonneted bat was not an objective of this study, the cryptic variation uncovered in Florida merits discussion. Eumops in Florida was first described from Pleistocene fossils (Molossides floridanus, Allen 1932). Live specimens from the Miami area were subsequently recognized as being conspecific with E. glaucinus (Ray et al. 1963) and considered a subspecies (E. g. floridanus; Koopman 1971). The bat was elevated to species status by Timm and Genoways (2004) based on overall larger size and several significantly larger cranial measurements. They examined 47 E. floridanus specimens; however only 5 were from outside Miami Dade County. Their morphological analysis included 28 of the 47 specimens (see Table 2 in Timm and Genoways 2004). Thus, the most recent morphological description of E. floridanus relies disproportionately on the MD lineage. Further systematic study of the Florida bonneted bat should include a reevaluation of the morphology from across the range and include genome resequencing or single nucleotide polymorphism (SNP) markers to afford sufficient resolution for greater insight into the demographic history, taxonomy and conservation decisions (Funk et al. 2019).
These results provide novel insight into the dynamics and history of this endangered species. One limitation of our current study is that we collected samples only from regions of southern Florida where bonneted bats have been previously captured. Expanded sampling efforts covering areas between regions would facilitate an assessment of the geographic extent of distinct mitochondrial lineages. Samples from currently unsampled regions, including the natural areas in the proximity of Miami Dade County, would help to clarify connectivity across the range and refine our understanding of the distribution of genetically distinct lineages of bats. Similarly, broadening the extent of geographic sampling of Eumops throughout the Caribbean and Mexico/Central America may uncover additional cryptic genetic structuring that will help to clarify the population history of the two lineages in Florida. Genomic data will provide a clearer picture of the demographic history of Florida bonneted bats and estimate adaptive divergence among Florida and Caribbean populations. These additional insights may prove critical for evaluating potential conservation actions, including whether the translocation of bonneted bats will be an effective management tool.