The Bahamas is an archipelagic nation consisting of 700 islands and 3,000 cays located in the tropical Western Atlantic (~80.5 km/50 mi) southeast of the United States of America and north of Cuba (Buchan, 2000; Figure 1). A total of 464 Nassau grouper fin clip samples were collected from 13 islands (n = 407) and one active FSA in The Bahamas (n = 57) between August 2014–February 2017 (Figure 1, Table 1). Tissue samples from living fish (n = 186) were collected with permission from the Department of Marine Resources during 10 research cruises (see Acknowledgements). Nassau grouper were captured via fish traps or divers using hand-held nets in water depths ranging from 4.6 to 32 m. Fish were slowly brought to the surface and processed using methods described by Stump et al. (2017). Nassau grouper were Floy™ tagged, measured (standard length, SL and total length, TL to the nearest 0.1 centimeter), and a small (~5 × 5 mm) tissue sample was removed from the anal or dorsal fin. Each fish was weighed (to the nearest 0.1 kg) and allowed to recover. Fish were then individually escorted back to the reef from which they were captured and monitored for 1–2 min to ensure there were no adverse impacts due to handling. Fin clips from recently killed specimens (n = 280) were also opportunistically collected from local fishing docks and landing sites throughout The Bahamas. If the capture location for a fish was not specified or was ambiguous, these individuals were treated as unknown. Samples were preserved in 95–100% ethanol (n = 452) or DMSO (n = 14) prior to genomic DNA extraction.

Genomic DNA for microsatellite analysis was extracted from fin clip tissues (n = 464) using the HotSHOT method (Truett et al., 2000) and its quality assessed via NanoDrop Spectrophotometry (ND−1,000). Twelve species-specific and three cross-species grouper microsatellite loci: Est33a, Est416, Est360, Est92, Est265, Est376, Est267, Est420, Est338, Est340, Est290, Est262, EACD08, EACB6, and EACD02 developed by Bernard et al. (2012) and Jackson et al. (2012) were selected for this research. Polymerase Chain Reaction (PCR) amplification was performed using a BIO-RAD MyCycler Thermal Cycler® in a total of three multiplexes each containing 1 μl DNA, 3 μl of RNase-free water and 5 μl of Qiagen HotStarTaq plus MasterMix, and 1 μl of primer mixture under identical touchdown thermal cycling conditions (Supplementary Tables 1, 2; Hamilton and Tyler, 2008).

Amplified PCR products were visualized and scored with a Beckman Coulter CEQ^TM 8,000 sequencer (Fullerton, California, USA) using the manufacturer specified internal size standard. Genotypes with more than four missing loci (i.e., 10 individuals) were excluded from further analyses. Microsatellite loci were assessed for the presence of null alleles, large allele dropout, scoring error and stutter using Microchecker v.2.2.3 (van Oosterhout et al., 2004; Supplementary Table 3). The log-likelihood ratio statistic and probability test in GENEPOP v.4.2 was used to test for linkage disequilibrium and conformance to Hardy-Weinberg equilibrium (Rousset, 2008). To reveal Type I errors in both tests, the false discovery rate (FDR) correction was also applied (Storey and Tibshirani, 2003).

To address whether patterns of genetic composition were due to temporal and/or Wahlund effects (i.e., heterozygosity deficiency due to unknown/cryptic population sub-structuring), we computed a range of genetic diversity parameters for Nassau grouper by sampling year (2014–2016) and for sub-adult and adult fish from two locations (Eleuthera and Exuma) with relatively robust sample sizes (Supplementary Table 4). Specifically, observed (H_O) and expected heterozygosity (H_E) were calculated using GenAlEx v.6.503 (Peakall and Smouse, 2006, 2012). The rarefaction method in HP-Rare v. 1.1 (Kalinowski, 2005) was used to calculate allelic richness (A_R) and private allelic richness (PA_R) for 2014, 2015, 2016, sub-adult and adult fish with minimum sample sizes of 14, six, 10, 14, and 14 genes per locus respectively. Because these measures produced similar results (Supplementary Table 4), data was pooled for subsequent analyses. Genetic diversity statistics including gene diversity and allelic richness were determined using microsatellite profiles for 454 fish (Table 2). A_R and PA_R for the pooled data set was computed with a minimum sample size of 14 genes for each location. Significant differences in genetic diversity parameters (H_O, H_E, A_R, and PA_R) among sample locations were calculated using the non-parametric Kruskal-Wallis rank sum test in R v.3.4.0.

Ascertainment bias was not a consideration in the comparison of our diversity results with those of Jackson A. M. et al., 2014. The current study had three loci in common with that of Jackson A. M. et al., 2014: EACB6, EACD02, and EACD08; these loci were originally designed for another grouper species (Gulf coney-Hyporthodus acanthistius) and successfully cross-amplified for Nassau grouper (Jackson et al., 2012). In both our study and that of Jackson A. M. et al., 2014 these loci were used with the same species (Nassau grouper).

Locations with sample sizes of ≥10 individuals were used in assessments of genetic differentiation, corresponding to a total of 421 Nassau grouper sampled from Abaco in the north to Great Inagua in the south. Significant differences in genetic differentiation (i.e., global and pairwise F_ST values; Weir and Cockerham, 1984) were computed using 10,000 permutations of the data set in Microsatellite Analyzer (Dieringer and Schlötterer, 2003). To examine the extent of genetic variance among sampled locations, GenAlEx v.6.503 was used to perform an analysis of molecular variance (AMOVA). Statistical significance was assessed using 999 permutations of the data. Principal coordinate analysis (PoCA) plots were generated from codominant genetic distance matrices for all samples (n = 454) initially and then for locations with sample sizes of ≥10 per location (n = 421) in GenAlEx v.6.503. A discriminate analysis of principal components (DAPC) was also constructed using the adegenet package in R to expose any fine-scale genetic partitioning (Jombart et al., 2010).

To visualize genetic population structure, the admixture ancestry model using correlated allele frequencies in STRUCTURE v.2.3.4 (Pritchard et al., 2000) was implemented with a pre- and post-burn-in of 75,000 and 150,000 Monte Carlo Markov Chain iterations respectively over 10 independent runs for inferred populations, K, ranging from 1 to 17. STRUCTURE analysis was re-run using a pre- and post-burn-in of 250,000 and 1,000,000 MCMCs respectively for K 1–10. The most likely partition of the dataset for both STRUCTURE runs was determined using the delta K (ΔK) statistic (Supplementary Figures 1, 2; Evanno et al., 2005) and was calculated using the web-based Structure Harvester v. 0.6.94 programme (Earl and vonHoldt, 2012). The PopHelper web application v.1.0.10 (with a ΔK of 2) was used to create a Structure plot of aligned and merged repeats for the 10 independent runs of K 1–10. Because marine species often show weak patterns of genetic differentiation, the locprior model (Hubisz et al., 2009) was also run with parameters identical to those used for the admixture model for K 1–10. Because the ΔK was the same for both models, only the results of the more conservative admixture model (which assumes no prior information about structure) are reported.

Finally, an assessment of whether genetic structure could be attributed to isolation-by-distance was performed using a Mantel test in GenAlEx v.6.503. Matrices of unbiased pairwise genetic distances F_ST/(1−F_ST) were compared with the logarithm of geographic distances between sites (km). Google Earth Pro was used to estimate the shortest straight line distance between islands. Significance was based on 999 permutations of the data.

Estimates of contemporary effective population size, N_e based on linkage disequilibrium were assessed using the random mating model with minimum allele frequencies of 0.05, 0.02, 0.01, and 0.005, as implemented in LDNE v.1.31 (Waples and Do, 2008). LDNE analyses were performed using sub-adult and adult fish from Eleuthera and Exuma and for the entire dataset. Additionally, we used STRUCTURE outputs for ΔK = 2 to assign individuals to the most likely genetic cluster and re-analyzed contemporary N_e using these two populations. M-ratio was used to test whether anthropogenic activities (e.g., overfishing) have led to declines in population size of Nassau grouper throughout the archipelago and within an active FSA. We calculated the M-ratio, the ratio of the number of microsatellite alleles compared to allele size range (Garza and Williamson, 2001). Critical M (M_c) values were computed using a pre-bottleneck N_e ranging from 10,000 to 15,000; corresponding to theta (θ) values between 5.6 and 120 and a mutation rate of μ = 5 × 10⁻⁴. We also tested mutation rates of 5.57 × 10⁻⁴, 2.0 × 10⁻³, and 1.5 × 10⁻⁴ derived from Common carp (Cyprinus carpio, Yue et al., 2007), Broadnosed pipefish (Syngnathus typhle, Jones et al., 1999) and Zebrafish (Danio rerio, Shimoda et al., 1999). The proportion of stepwise mutations was set to 0.78 and the mean size of non-stepwise mutation (δ_g) to 3.1 following the suggestions of Peery et al. (2012).

Temporal changes in N_e were analyzed using the VarEff package in R to determine when historical population declines occurred (Nikolic and Chevalet, 2014). The VarEff model uses the coalescent method and approximate likelihoods to derive posterior distributions of N_e over a specified past generation time. Initially, this analysis was performed assuming a two-phase model, T, with a commonly used microsatellite mutation rate, μ = 5 × 10⁻⁴ (Estoup et al., 2002), and burn-in of 10,000 over the past 1,000 generations. Additional parameters included setting the number batch to 50,000, length and space batch to 10, with an acceptance rate of 0.25 and diagonale of 0.5, following recommendations from Nikolic and Chevalet (2014). The model was also run using mutation rates of 5.57 × 10⁻⁴ (C. carpio), 2.0 × 10⁻³ (S. typhle) and 1.5 × 10⁻⁴ (D. rerio) respectively, under the same parameters. The procedure was re-run modifying the number of times N_e changed (JMAX 2 and 3) and assumed prior values of N_e (NBAR 10,000 and 15,000) for each of the mutation rates listed above (Supplementary Figures 3–6).

Following the removal of 10 individuals with missing data for four or more loci, genotypes from 454 Nassau grouper were retained for analysis. Evidence of null alleles and homozygote excess or heterozygote deficiencies were detected in eight loci, but these were inconsistent across sample locations (Supplementary Table 3). Three occurrences of stuttering and potential scoring error were detected for loci EACB6 and Est360. However, no large allele dropout was detected, and stuttering and scoring error were negligible across populations (Supplementary Table 3). Tests for linkage disequilibrium revealed associations between 10 pairs of loci for fish collected from unknown locations within The Bahamas and Abaco, Eleuthera and Exuma. Samples at three locations, Abaco, Exuma and Long Island, showed possible deviations from Hardy Weinberg equilibrium. Locus EACB6 exhibited significant deviations from HWE after FDR correction across ≥50% of populations and 10 instances of homozygote excess, which can be indicative of null alleles. Accordingly, genotypes for this locus were excluded from subsequent analyses.

Microsatellite loci were highly polymorphic, with the number of alleles per locus ranging from 21 (Est33a) to 42 (EACD02; Supplementary Table 5). A total of 394 alleles were detected, averaging 28.14 across all 14 loci. Measures of genetic diversity and allele size range varied among loci and across sample locations (Table 2, Supplementary Tables 4, 5). Gene diversity by sampling year and life stage ranged from H_E = 0.746–0.873 and from 0.838 to 0.877 respectively (Supplementary Table 4). Estimates of allelic richness were variable, ranging from A_R = 4.4–7.84 (Supplementary Table 4). Inbreeding coefficients also varied by sampling year and life stage, with higher F_IS values observed for sub-adults vs. adults (Supplementary Table 4). Allelic richness for the pooled dataset (adjusted for sample size) was similar among locations, with values ranging from 6.99 to 7.95 (Table 2). Private allelic richness varied from 0.26 to 0.89. Mean allelic richness was highest in Long Island (7.95) and lowest in Eleuthera (6.99). Private alleles were identified in all locations, with more private alleles present in Exuma (mean PA_R of 0.98) specimens (Table 2). Heterozygosity, H_E was generally high across the sampling area (ranging from 0.753 to 0.873), with the exception of Bimini (0.321), which contained the smallest samples size (Table 2). Mean H_O was lowest in Bimini (0.643) and highest in New Providence (0.872). Long Island had the highest mean inbreeding coefficient, (F_IS = 0.117) and New Providence had the lowest (F_IS = −0.028). Kruskal-Wallis tests showed no significant differences in H_E, H_O, A_R, and PA_R among locations (p > 0.05).

PCoA, DAPC, and STRUCTURE analyses revealed no distinct clustering or definitive geographic structure in the genetic composition of Bahamian Nassau grouper when examined by island or individual (Figures 2, 3). However, there was a low but significant difference in the overall F_ST value across samples and loci (Global F_ST = 0.00236, p = 0.0001; Table 3). Genetic differentiation based on pairwise F_ST values ranged from −0.00097 to 0.01194 and was greatest between Eleuthera and Long Island. Pairwise F_ST values were significant between only 8 population pairs (Table 3). The Mantel test revealed no significant isolation-by-distance (rxy = −0.208, p = 0.255; Figure 4). AMOVA results indicated that most genetic variability, 90%, occurs within individuals, with 10% among individuals (Supplementary Table 6).

Contemporary linkage disequilibrium based estimates of N_e varied among locations and in association with the allele frequency used (Table 4, Supplementary Table 7). N_e also varied by sampling year, life stage and population (Supplementary Table 7). Values based on a critical allele frequency of 0.02 are reported (Table 4) following recommendations by Waples and Do (2010). N_e was lowest for New Providence (−888.3), but this was likely due to the small sample size. To investigate this, we randomly sub-sampled data (n = 10) from Abaco, Eleuthera and Andros respectively and re-ran LDNE analysis, which also produced negative N_e estimates (Supplementary Table 8). Excluding New Providence, N_e estimates were lowest for Abaco (~375) and highest for Andros (~2,978).

The VarEff modeling approach using μ = 5 × 10⁻⁴ revealed historically very large and stable effective population sizes for Nassau grouper within The Bahamas (Figure 5). The upper estimates of median N_e ranged from ~111,000 in Long Island to ~315,000 in Andros. N_e estimates for The Bahamas peaked at ~227,000 and began to gradually decline around 400 generations ago. Similar patterns of population decline were identified for all islands in The Bahamas and Hail Mary around the same period, with a severe decline occurring within the last 150 generations (Figure 5). The upper estimates for Nassau grouper within the last 100 generations were variable, ranging from ~135,742 in Great Inagua to 5,024 Long Island (Figure 5). Lower and upper contemporary N_e estimates during the last 50 generations varied from ~830 to 1,514 in Great Inagua and ~2,672 to 4,539 in Exuma. The upper estimate for The Bahamas during this timeframe was ~2,016 (Figure 5).

Tests using faster and slower mutation rates shifted the timing of genetic bottlenecks (Supplementary Figures 3–6). Under a faster rate of mutation (μ = 2 × 10⁻³), bottlenecks would have occurred within the last 50 generations, whereas under the slower mutation rate (μ = 1.5 × 10⁻⁴), declines occurred around 300–400 generations ago (Supplementary Figure 6). All model simulations consistently revealed a pronounced decline in Nassau grouper within the past 400 generations for all islands (Supplementary Figures 3–6). The only exceptions were Exuma and Long Island, which under the slower mutation rate, showed population declines within the last 450 generations (Supplementary Figure 5).

Dating when bottlenecks occurred varied in accordance with the generation time used. We explored this using generation times ranging from 4 to 9 years. These times encompass the earliest age of sexual maturity, 4 years (Sadovy and Colin, 1995), the suggested generation time calculated for an unexploited Nassau grouper FSA, 9 years (Sadovy and Eklund, 1999), and the approximate age at which Bahamian Nassau grouper undertake their first spawning migration, 7 years (Dahlgren et al., 2016a). Comparisons of historical vs. contemporary derived VarEff N_e estimates for The Bahamas indicate a dramatic reduction in N_e over the last 400 generations, which would have occurred between ~1,600–3,600 years ago (corresponding to generation times of 4 and 9 years respectively). M-ratio estimates also varied in relation to the mutation rate used. Generally, under the fastest mutation rate, no bottlenecks were detected, however, under more realistic mutation rates, M-ratio values were below critical M values for both pre-bottleneck thetas (Supplementary Table 9). Overall, VarEff outputs were consistent with M-ratio analyses, showing historic bottlenecks on several islands and throughout The Bahamas (Table 4, Figure 5).

Results from microsatellite analysis of genotyped Nassau grouper revealed weak genetic differentiation, lacking definitive geographical population structure. Patterns of weak or no genetic subdivision at varying spatial and temporal scales using microsatellites and other molecular markers have been reported for a number of Epinephelids, e.g., Camouflage grouper, E. polyphekadion (Rhodes et al., 2003); Hawaiian grouper, E. quernus (Rivera et al., 2004); Red grouper, E. morio (Zatcoff et al., 2004); Dusky grouper, E. marginatus (Maggio et al., 2006; Schunter et al., 2011; Buchholz-Sørensen and Vella, 2016) and Goliath grouper, E. itajara (Silva-Oliveira et al., 2008). We expected to see a stronger signal of genetic differentiation separating the northern and southern Bahamas given that the Exumas were purported to present a barrier to gene flow (Jackson A. M. et al., 2014). However, Nassau grouper sampled throughout the archipelago were overlapping in DAPC analysis and no geographical clustering was apparent in either PCoA or Bayesian-based Structure analyses. Moreover, a closer examination of fish sampled north to south across (~200 km) in the Exuma chain, revealed considerable genetic mixing, which is inconsistent with the hypothesis of reduced gene flow. In contrast, low estimates of genetic differentiation and similarly high levels of genetic diversity suggest moderate to extensive gene flow or connectivity throughout The Bahamas.

There are several plausible explanations for the patterns of genetic variation observed. Firstly, pelagic larval dispersal (PLD) and recruitment success are affected by many factors including interspecific competition, predation, larval behavior, oceanographic processes, environmental conditions (e.g., water temperature) and habitat suitability (Colin, 1992, 2012; Shenker et al., 1993; Choat, 2012; Hamner and Largier, 2012; Kough et al., 2016). The stochastic nature of these biotic and abiotic factors helps to determine population structure and demographic connectivity of adults. For example, in bonefishes (Albula spp.), the convergence of larval dispersal and ecological processes has been shown to influence genetic population structuring (Wallace, 2015; Wallace and Tringali, 2016). Nassau grouper have a protracted PLD of up to 50 days, with a mean of 42 days (Tucker and Woodard, 1991; Colin et al., 1997). In the Caribbean, in situ monitoring of sub-surface currents along with modeled trajectories of ocean currents and larval dispersal of Nassau grouper have shown varied recruitment patterns including advection from spawning sites and retention via eddies (Colin, 1992; Heppell et al., 2009, 2011). However, behavioral traits and ecology of larval Nassau grouper in tandem with local oceanographic-driven current patterns have not been explicitly investigated for The Bahamas.

Secondly, adult Nassau grouper can undertake long migrations (>300 km) to participate in annual transient FSAs (Bolden, 2000; Domeier, 2012; Dahlgren et al., 2016a; Stump et al., 2017). Telemetry studies have demonstrated that adults exhibit philopatry to specific FSAs (Starr et al., 2007; Dahlgren et al., 2016a). It is likely therefore, that the population dynamics of Nassau grouper is influenced by both adult migration to and from FSAs via continental shelf margins (Bolden, 2000; Stump et al., 2017), FSA site fidelity (Sala et al., 2001; Starr et al., 2007; Dahlgren et al., 2016a), larval dispersal (Colin, 1992; Shenker et al., 1993) and ocean currents (Colin et al., 1997; Paris, 2009; Paris et al., 2013). Ocean circulation patterns in The Bahamas show predominately northerly flowing currents (Figure 1); however, these patterns are not static and exhibit considerable complexity (Gunn and Watt, 1982; Colin, 1992, 1995). Strong winds and cross-continental shelf currents driven by winter cold fronts have been shown to facilitate larval dispersal from FSAs to recruitment habitats (e.g., Shenker et al., 1993; Colin, 1995; but see Colin et al., 1997). Additionally, mesoscale gyres have been documented in the Exuma Sound, enabling entrainment of larvae through to recruitment stage within this area (Colin, 1995). Collectively, the processes of adult migration and larval dispersal have the potential to increase gene flow, thus contributing to the apparent maintenance of genetic diversity despite ongoing anthropogenic pressures (e.g., FSA fishing).

However, the ability to detect genetic differentiation can be masked by high levels of polymorphism, genetic diversity and large N_e, which can act to retard mutation-drift equilibrium (O'Reilly et al., 2004; Hellberg, 2009). In this study, both patterns of genetic partitioning (K = 2) and a high degree of genetic diversity were found. The measures of genetic diversity reported here (H_E = 0.753–0.873) using neutral microsatellite markers are comparable to those found for Nassau grouper in the Caribbean (H_E = 0.608–0.966, Jackson A. M. et al., 2014; H_E = 0.32–0.91, Bernard et al., 2016) and to other grouper species (H_E = 0.703–0.762, Schunter et al., 2011; H_E = 0.04–0.886, Buchholz-Sørensen and Vella, 2016). However, given the extent of overfishing and other anthropogenic activities, which have been shown to collapse FSAs (Aguilar-Perera, 2006; Sadovy De Mitcheson et al., 2008) and alter essential marine habitats (Buchan, 2000; Gardner et al., 2003; Pandolfi et al., 2003; Jackson J. B. C et al., 2014), the high levels of genetic diversity observed were somewhat surprising. Indeed, genetic diversity measures for Hail Mary, an active but fished Nassau grouper FSA, were similar to those for the rest of The Bahamas, suggesting no loss in genetic diversity. In other exploited grouper species, e.g., E. itajara, molecular analyses employing microsatellites have provided evidence for reduced genetic variability (Silva-Oliveira et al., 2008). However, this is not always the case, as demonstrated by Zatcoff et al. (2004) for E. morio and Scamp, Mycteroperca phenax; in their study, both species were genetically homogeneous and genetically diverse, despite intense fishing pressure.

Typically, high heterozygosity and allelic richness is beneficial because, potentially, it confers a fitness advantage and has implications for both short and long-term adaptations to both natural and anthropogenic stressors (Bouzat, 2010; Bernatchez, 2016). We found no significant difference in allelic richness in fish sampled from across the islands, but these values were low and may indicate possible loses of genetic diversity (Pinsky and Palumbi, 2014). The allelic richness results reported here are analogous to other published data for Nassau grouper (Bernard et al., 2016) and E. marginatus (Schunter et al., 2011), but higher than values reported from E. quernus (Rivera et al., 2011). In populations that undergo a bottleneck, rare or unique alleles may be lost, deleterious alleles may become fixed and inbreeding is more likely to occur (Garza and Williamson, 2001). In our study, coefficients of inbreeding were high in Abaco (F_IS = 0. 098) and Long Island (F_IS = 0.117), suggesting that inbreeding of Nassau grouper may be occurring at these locations. However, it is important to highlight that samples from both these islands showed possible deviations from Hardy-Weinberg Equilibrium. Additionally, sub-adults comprised the majority of Nassau grouper from Abaco and most of these fish were caught around the same area (Table 1). This is likely due to Wahlund effects (from the incorporation of mixed cohorts and/or population sub-structuring), providing some explanation for one of the highest F_IS values found in this study (Table 2, Supplementary Table 4). Similar patterns may not have been observed if Nassau grouper were sampled over a broader geographic area and across different life stages. In contrast, both sub-adults and adults were sampled over a wider area in Long Island, although mostly around the northern part of the island.

In most instances, LDNE derived N_e estimates varied in accordance with the critical allele frequency used, with the smallest allele frequency (0.005) and smaller sample sizes generating the lowest estimates for N_e (Table 4, Supplementary Table 7). Analysing the data by sampling year, life stage, and population greatly reduced sample sizes and produced mostly negative results, which is either an artifact of sampling or suggest that N_e is indistinguishable from infinity (Waples and Do, 2010; MacBeth et al., 2013). Recently, MacBeth et al. (2013) compared empirical and simulated data of various population sizes to estimate N_e and examine the accuracy of these estimates for Spanish mackerel (Scomberomorus commerson), demonstrating how sample size and the presence of conspecific migrants from genetically differentiated populations could result in reduced or negative N_e estimates.

In the present study, contemporary values of N_e were generally low, but higher estimates were observed in Andros, Hail Mary, and Great Inagua (Table 4). M-ratio and VarEff analyses provided strong evidence in support of both recent and historic population declines. Recent bottlenecks were detected for Abaco, Eleuthera, Exuma and Hail Mary fish despite showing high levels of genetic diversity consistent with the rest of The Bahamas. This is particularly interesting for Hail Mary, which was also the location with the second highest N_e estimate (Table 4). Fishing pressure on Nassau grouper around Eleuthera and Exuma, and FSA fishing at Hail Mary are likely to have contributed to the observed bottlenecks. The majority of contemporary LDNE-based N_e values we derived are smaller than those reported (using different microsatellite markers) for E. marginatus—an endangered grouper species with similar life history characteristics to Nassau grouper (Schunter et al., 2011). However, the smallest contemporary N_e value observed (Abaco: 375) was similar to those reported for other teleosts (e.g., Santer and Slinger sea bream, Cheimerius nufar and Chrysoblephus puniceus; Coscia et al., 2016).

We recognize that most Nassau grouper used in this study were sampled over consecutive years and across life stages with a few exceptions (Table 1), thus potentially violating assumptions regarding discrete generations (Waples, 2006; Waples et al., 2014). Variability in N_e estimates has been linked to differences in reproductive success and/or the methods used in their computation (Hare et al., 2011; Waples et al., 2013, 2014; Coscia et al., 2016; Waples and Anderson, 2017). For example, Waples et al. (2014) demonstrated that N_e estimates based on a variety of age-structured samples are likely to be underestimated due to Wahlund effects. In the present study, both LDNE and VarEff produced similar estimates of contemporary N_e, providing additional confidence in modeled results. However, based on potential biases associated with N_e estimates, results should be interpreted with caution.

The timing of contemporary N_e losses (i.e., within the last 100 generations) coincides with known anthropogenic impacts, e.g., fishing, habitat degradation to Bahamian marine resources (Smith, 1972; Lang et al., 1988; Buchan, 2000; Jackson J. B. C et al., 2014). Historical estimates of Nassau grouper N_e were very large. Jue (2006) also reported a large historical N_e of ~30,000 for Gag (Mycteroperca microlepi), another exploited aggregating grouper species. M-ratio values for Bahamian Nassau grouper were below the critical M (0.68) recommended by Garza and Williamson (2001), highlighting the severity of population declines. In contrast, Bernard et al. (2016) reported a weakly bottlenecked Nassau grouper population from an FSA in the USVI.

More remarkably, if we examine historical N_e reductions in The Bahamas, assuming a generation time of 9 years, the first signs of population decline would have occurred 3,600 year ago. Although Florida is within close proximity to The Bahamas, during the 1800s most of the population resided in northern Florida (Smith, 2005). Consequently, it seems unlikely that the native Indian tribes residing there 3,600 years ago would have traveled to The Bahamas to fish with access to ample marine resources around the shallow Floridian coast. The earliest known inhabitants to populate the Bahamian archipelago were the Lucayan or Arawak Indians (900–1500 AD) who survived by engaging in low technology subsistence fishing and farming, before their eradication in the 1500s (Buchan, 2000). A decline in N_e 3,600 years ago pre-dates the existence of human occupancy in The Bahamas and implies that historical population bottlenecks were likely due to natural as opposed to anthropogenic events.

According to marine geological records, four major glacial events happened during the Pleistocene (2.6 Ma to 11,700 years. ago), resulting in massive (120 m) fluctuations in sea level, which impacted The Bahamas (Sealey, 1994). This was punctuated by smaller variations in sea level (5–6 m) above current conditions during inter-glacial periods (White et al., 1998) that also would have affected marine organisms. Unfortunately, no historical genetic material was available to establish baseline genetic variability and cross-validate modeled estimates of temporal changes in the effective population sizes and contemporary genetic diversity values of Nassau grouper. However, previous studies have shown that the demography, distribution and genetic composition of species can be strongly influenced by geological processes and associated eustatic sea level changes (e.g., Hewitt, 2000; Pellissier et al., 2014; Brüniche-Olsen et al., 2017). For example, an analysis of >6,000 fish species showed strong correlations between current species richness and distribution patterns with historic climatic events that negatively impacted reef habitats (Pellissier et al., 2014). This finding was particularly pronounced for Caribbean reefs and strongly associated reef species, which experienced drastic reductions in water temperatures during past glacial events (White et al., 1998; Hewitt, 2000; Pellissier et al., 2014). Although speculative, it is possible that historical bottlenecks experienced by Bahamian Nassau grouper ~3,600 years ago could have been associated with severe climatic perturbations. Collectively, these results suggest that both climate-driven biological processes and human exploitation have impacted the genetic composition of Nassau grouper.

The concept of management units (MUs) is based upon populations with genetically distinct population dynamics and demographics, resulting in management strategies tailored to preserve genetic variability within units (Palsbøll et al., 2007). Two putative genetic clusters were revealed by STRUCTURE analysis, but these do not appear related to geography. Possible explanations for this finding include panmixia, temporal differences in spawning between two sympatric stocks, and secondary contact from previously isolated Nassau grouper populations (e.g., Reid et al., 2016). However, our study was not designed to investigate these hypotheses, and to address it would require expanding spatial and temporal sampling of Nassau grouper throughout The Bahamas and the Caribbean. In the absence of demographic genetic structure, N_e is a viable alternative to use for conservation management. Prior to this study, no estimates of N_e were available for Nassau grouper. Genetic analyses have been successfully used to demonstrate the genetic consequences of exploitation through reductions in N_e (e.g., Hauser et al., 2002; Hoarau et al., 2005; Poulsen et al., 2006), a finding we observed in the present study. The protection of FSAs has been identified as a biodiversity target for The Bahamas (Moultrie, 2012) and has been incorporated into Marxan analyses for the establishment of new MPAs (Moultrie, 2012; Moss and Moultrie, 2014). Given the substantial reduction in N_e and the low contemporary values of allelic richness, strengthened national, regional and international relationships are recommended to tackle illegal fishing, to eliminate FSA targeted fishing and to establish MPA networks to protect marine habitats that are of critical importance to rebuild Nassau grouper stocks.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer RH and handling Editor declared their shared affiliation.