RAD-Seq Analysis and in situ Monitoring of Nassau Grouper Reveal Fine-Scale Population Structure and Origins of Aggregating Fish

Sherman, Krista D.; Paris, Josephine R.; King, Robert Andrew; Moore, Karen A.; Dahlgren, Craig P.; Knowles, Lindy C.; Stump, Kristine; Tyler, Charles R.; Stevens, Jamie R.

doi:10.3389/fmars.2020.00157

ORIGINAL RESEARCH article

Front. Mar. Sci., 19 March 2020

Sec. Marine Conservation and Sustainability

Volume 7 - 2020 | https://doi.org/10.3389/fmars.2020.00157

RAD-Seq Analysis and in situ Monitoring of Nassau Grouper Reveal Fine-Scale Population Structure and Origins of Aggregating Fish

$\r\nKrista D. Sherman,,,*$ Krista D. Sherman^1,2,3,4*

Josephine R. Paris¹

Robert Andrew King¹

Karen A. Moore¹

Craig P. Dahlgren^3,4

Lindy C. Knowles³

Kristine Stump^2,5

Charles R. Tyler¹

Jamie R. Stevens¹

¹Biosciences, University of Exeter, Exeter, United Kingdom
²Daniel P. Haerther Center for Conservation and Research, John G. Shedd Aquarium, Chicago, IL, United States
³Science and Policy, Bahamas National Trust, Nassau, Bahamas
⁴Perry Institute for Marine Science, Waitsfield, VT, United States
⁵Field Lab Consulting, LLC, Palmetto, FL, United States

Nassau grouper (Epinephelus striatus, Bloch 1792) are globally critically endangered and an important fishery species in The Commonwealth of The Bahamas (hereafter The Bahamas) and parts of the Caribbean, with an urgent need for better management and conservation. Here, we adopted a combined approach, integrating restriction-site-associated DNA sequencing (RAD-seq) and acoustic telemetry to establish country-wide demographic structure, diversity and connectivity, and the origins of Nassau grouper using an active fish spawning aggregation (FSA) in the central Bahamas. RAD-seq analysis of 94 Nassau grouper sampled from nine locations in The Bahamas generated a working dataset of 13,241 single nucleotide polymorphisms (SNPs). Similar levels of genetic diversity were found among sampled locations. Evidence of population sub-structuring across The Bahamas was demonstrated and supported by discriminate analysis of principal components (DAPCs) along with analyses of molecular variance (AMOVAs). Associated acoustic telemetry data indicated Nassau grouper tagged at an active FSA in the central Bahamas during the 2016–2017 spawning season migrated to the Exumas at the conclusion of the spawning period. Telemetry data suggest the likely origins of five individuals, which traveled one-way distances of up to 176 km from the FSA in the central Bahamas to two sites within a no-take marine protected area (MPA). Analyses of high-resolution SNP markers (including candidate loci under selection) illustrated patterns of spatial structure and genetic connectivity not reflected by telemetry data alone. Nassau grouper from Exuma and Long Island appear to have genetic signatures that differ from other islands and from the Hail Mary FSA. Collectively, these findings provide novel information on the intraspecific population dynamics of Nassau grouper within The Bahamian archipelago and within an active FSA.

Introduction

Approximately 25,000 of the species assessed worldwide have been classified as threatened and are at risk of extinction (IUCN, 2017), including 12% of grouper (Epinephelidae) species (Sadovy de Mitcheson et al., 2013). Changing this trajectory will require coordinated research and management approaches to develop a more robust understanding of population demographics and ecology at both small and large spatial scales. Current approaches to assess and manage fish species include analyses of fisheries dependent and independent data for stock assessments (Thurstan et al., 2015; Egerton et al., 2017; Pauly and Zeller, 2017), examination of habitat use and migratory behavior (Austin, 2007; Pittman et al., 2014; Aswani et al., 2015), and modeling predicted responses to climate change (Cheung et al., 2012; Crozier and Hutchings, 2014; Sunday et al., 2014). Designating and monitoring marine protected areas (MPAs) or reserves (Hughes et al., 2007; Mumby and Harborne, 2010; White, 2015), establishing stock enhancement programmes (Griffiths et al., 2011; Lorenzen et al., 2013; Abelson et al., 2016), recommending sustainable harvest regulations (Dann et al., 2013; Bozec et al., 2016; Hazen et al., 2016), and evaluating and incorporating the role of adaptation in survival and reproductive fitness (Eizaguirre and Baltazar-Soares, 2014; Delmore et al., 2015; Calosi et al., 2016) are other examples of valuable approaches for species management.

Nassau grouper (Epinephelus striatus, Bloch, 1792), exemplify a species requiring immediate management due to their IUCN critically endangered status, and economic, cultural, and ecological significance (Sadovy de Mitcheson and Colin, 2012; Carpenter et al., 2015; Sherman et al., 2016, 2018a). Nassau grouper are inherently vulnerable to population declines because of their reproductive and life history traits (Coleman et al., 2000; Dahlgren and Eggelston, 2001; Domeier, 2012). More specifically, the predictable occurrence of lunar-associated migrations to known fish spawning aggregations (FSAs) for annual reproduction with conspecifics increases their susceptibility to overexploitation (Domeier, 2012; Sadovy de Mitcheson and Erisman, 2012; Cheung et al., 2013). As a result, many historic Nassau grouper FSAs have disappeared and others have declined by orders of magnitude in the last few decades (Olsen and LaPlace, 1979; Sadovy de Mitcheson et al., 2013).

Investigations of Nassau grouper FSAs have described the dynamic nature of these migrations – particularly, intraspecific differences in the timing of migrations (in relation to the winter solstice) and distances traveled to and from FSAs (Smith, 1972; Colin, 1992; Bolden, 2000; Starr et al., 2007; Dahlgren C. P. et al., 2016). Additionally, acoustic telemetry studies suggest that adults exhibit site fidelity to specific FSAs (Heppell et al., 2009; Kadison et al., 2010; Dahlgren C. P. et al., 2016) and tend to home to their respective reefs (Bolden, 2000; Dahlgren C. P. et al., 2016). However, FSAs are comprised of hundreds to thousands of individuals that migrate from multiple locations (Domeier, 2012), and larvae are dispersed with ocean currents, which can facilitate genetic mixing of different stocks during the spawning season (Colin, 1992, 2012; Molloy et al., 2012). Indeed, previous research has demonstrated how variability in migratory behavior can influence genetic connectivity (Dann et al., 2013; Delmore et al., 2015; Moore et al., 2017). The Bahamian archipelago represents an important area for Nassau grouper (Dahlgren C. et al., 2016; Sherman et al., 2016), yet data on the current status of Nassau grouper, its spawning stocks, patterns of gene flow, population structure and the origins of migrating adults within the country are limited (Dahlgren C. P. et al., 2016; Sherman et al., 2016, 2017; Stump et al., 2017). Such information is critical in order to effectively manage this critically endangered species.

In recent decades, conservation biology has increasingly incorporated genetic tools to assist with fisheries management (Reiss et al., 2009; Flanagan et al., 2017; Paris et al., 2018). Several reviews have highlighted the advantages and limitations of molecular approaches for both marine and freshwater species (e.g., Selkoe et al., 2016; Bernatchez et al., 2017; Lowe et al., 2017). In particular, highly connected marine species, characterized by weak genetic subdivision and large effective population sizes (N_e), pose difficulties for population assignment and ascribing management units (MUs) (Allendorf et al., 2010), in comparison to freshwater species (e.g., Vidal and García-Marín, 2011; Whiteley et al., 2011). For Nassau grouper, most assessments of genetic differentiation have been based on relatively small numbers of nuclear markers, i.e., mtDNA or microsatellites (Jackson et al., 2014; Bernard et al., 2016; Sherman et al., 2017), which may provide insufficient resolution to delineate fine-scale population structure.

Restriction-site-associated DNA sequencing (RAD-seq) can be used to identify and genotype large numbers of single nucleotide polymorphisms (SNPs). Increased coverage throughout the genome using thousands of SNP markers have substantially improved the capacity to unravel genetic diversity and divergence in species with both strong (Perrier et al., 2013) and weak (Larson et al., 2014; Benestan et al., 2015) population structure. While SNPs have been identified and utilized in Caribbean-wide assessments of genetic subdivision for Nassau grouper (Jackson et al., 2014), such an approach has not been used to explore fine-scale patterns of genetic variation in The Bahamas. Recent research based on microsatellite data showed no geographic population structure, suggesting that contemporary genetic connectivity may be higher than previously thought (Sherman et al., 2017). Better characterization of the spatiotemporal population dynamics of Nassau grouper through genome-wide assessments of genetic diversity and in situ investigations of FSAs and spawning migration patterns are critical to develop effective conservation management strategies (e.g., optimizing MPA networks). Emerging research has begun to integrate genetics and telemetry data, providing complimentary approaches to understand migration, spawning behavior and to elucidate patterns of connectivity for aquatic species [e.g., Arctic ringed seals (Pusa hispida, Schreber, 1775) (Martinez-Bakker et al., 2013); lake sturgeon (Acipenser fulvescens) Rafinesque, 1817 (Donofrio et al., 2018); Arctic char (Salvelinus alpinus) Linnaeus, 1758 (Moore et al., 2017)].

In this study, we used both RAD-seq and acoustic telemetry data to provide a more detailed insight into the spatial population dynamics and intraspecific genetic population structure of Nassau grouper in The Bahamas. The specific objectives were to (1) use SNPs to investigate fine-scale genetic variation throughout The Bahamas, (2) combine population genomics and acoustic telemetry to determine the origins of Nassau grouper migrating to Hail Mary, an active FSA in the central Bahamas and (3) investigate the role of selection on the genetic population structure of Nassau grouper across The Bahamas. These data are required to assess spatial and temporal changes in the health of Nassau grouper, inform national fisheries management strategies and direct future research and monitoring (Sherman et al., 2018a).

Materials and Methods

Sampling, DNA Extraction, and Quantification for RAD-Seq

Nassau grouper fin clips (n = 96) were collected according to methods described by Sherman et al. (2017). Samples originated from eight islands (Abaco, Andros, Eleuthera, Exuma, Great Inagua, Long Island, New Providence, and Ragged Island) and the Hail Mary FSA in The Bahamas (Figure 1). Fin clips were preserved in 95–100% ethanol prior to genomic DNA extraction. Genomic DNA was isolated using the Qiagen DNeasy Blood and Tissue Kit^® according to the manufacturer’s instructions. DNA was quality assessed via agarose gel electrophoresis and quantified using the dsDNA BR Qubit Quant-iT assay (Invitrogen).

FIGURE 1

Figure 1. Sampling locations of Nassau grouper (blue circles) and prevailing currents (black arrows) showing the extent of the Exuma Cays Land and Sea Park (rectangle).

RAD Library Development and Sequencing

Libraries were prepared using a modification of the original single-digest RAD protocol (Baird et al., 2008). RAD-seq libraries were prepared by the Exeter Sequencing Service (University of Exeter, United Kingdom) using the SbfI (New England Biolabs) enzyme. Libraries were sequenced on two lanes (48 individuals per lane) of Illumina HiSeq 2500 using single-end sequencing. Additional details on library preparation and sequencing methodology can be found in Supplementary Material 1.

SNP Discovery and Genotyping

FastQC was used to inspect the quality of raw sequence reads (Andrews, 2010). Stacks v. 1.42 (Catchen et al., 2011, 2013) was used for RAD-seq data processing, SNP discovery and standard population genetic statistics. Raw sequence data were de-multiplexed and cleaned (e.g., quality ≤ 10) and truncated to 96 bp using process_radtags. Following recommendations by Paris et al. (2017), de novo assembly parameters were optimized (m: 3–7, M: 1–7, and n: 2–4; see Supplementary Figure 1).

The denovo_map pipeline was run using optimal parameters of m5, M3, and n3. Two samples – one from Andros (AN026) and one from Hail Mary (LI239) – were removed from subsequent analyses due to low coverage (Supplementary Tables 1, 2). The program, cstacks was used to rebuild a catalog of loci using individuals with coverage depths of ≥ 15x (n = 59) to reduce the likelihood of genotyping recurrent sequencing error. Loci from all individuals were then matched back to this catalog using sstacks to generate genotypes for each individual. populations was used to filter loci (p = 5, r = 0.8, min-maf = 0.1, max_obs_het = 0.7), resulting in a final dataset consisting of 10,031 polymorphic loci (and ∼13 k SNPs), which was exported as GENEPOP (Raymond and Rousset, 1995) and vcf (Danecek et al., 2011) file formats. Further data conversion was implemented using PGD Spider v. 2.0.1.0 (Liscjher and Excoffier, 2012).

Loci were assessed for conformity to Hardy-Weinberg equilibrium (HWE) using GENEPOP v.4.2 (Rousset, 2008), and multiple sample testing was corrected using the false discovery rate (FDR) (Storey and Tibshirani, 2003). To address linkage disequilibrium (LD) between SNPs and to account for potential ascertainment bias only a single random SNP was selected from each RAD-tag (–write_single_snp in populations).

Population Diversity, Differentiation, and Structure

Standard diversity measures, including allele frequencies, expected and observed heterozygosity (H_E and H_O), the inbreeding F-statistic (F_IS), numbers of polymorphic loci, and global and pairwise values of F_ST (Weir and Cockerham, 1984) were computed in GenoDive v. 2.0b27 (Meirmans and Van Tienderen, 2004).

For analyses sensitive to outlier loci, a neutral dataset containing one random SNP from each locus was used (10,031 loci; see section Results). Assessments of genetic population structure were performed with a discriminate analysis of principal components (DAPCs) using the adegenet package in R v. 1.0.143 (Jombart, 2008; Jombart et al., 2010). The number of principal components employed to create each DAPC was based on inspections of the optimal α-score and cross-validation outputs (Jombart et al., 2010). To avoid over-fitting the data, the number of principal components retained for analysis was based on the lower of the two outputs. DAPCs were constructed using the optimal α-score value (40 PCs) and three discriminant functions.

Results of the sub-structure suggested by DAPC analysis was used to develop and test hypotheses for analyses of molecular variance (AMOVAs). The infinite alleles model (IAM) in GenoDive was selected to perform three AMOVAs based on 10,000 iterations of the neutral SNP dataset to test whether genetic variation was greater between (1) The Bahamas and Exuma, (2) The Bahamas and Long Island, or (3) The Bahamas versus Exuma and Long Island.

Identifying Loci Under Selection

Because F_ST values were low and non-significant (see section Results), two environmental association tests: latent factor mixed models (LFMMs, Frichot and François, 2015) and redundancy analysis (RDA, Forester et al., 2018) were employed to screen for potentially adaptive loci. We obtained environmental data for eight sample sites (excluding the Hail Mary spawning aggregation) from the Bio-ORACLE database (Tyberghein et al., 2012). Data for each site [salinity (SAL), minimum and maximum sea surface temperature (SSTmin and SSTmax), current (CURRmin and CURRmax) and primary productivity (PPmin and PPmax)] were extracted using the R package sdmpredictors (Bosch et al., 2018). Environmental variables were selected based on their significance to the timing of spawning migrations, spawning, larval distribution and survival for this species (Colin, 1992, 2012; Eggleston, 1995). Pearson’s r correlations were used to check for relationships between the measured environmental variables.

LFMM v 1.5 was used to detect correlations between genotypes and environmental variables (Frichot et al., 2013). LFMM analysis was run for 10,000 iterations and a burn-in of 5,000 reps, as suggested by Frichot et al. (2013). The number of latent factors was set at K = 3, as suggested by estimating the minimal cross-entropy criterion of the snmf function (iterating through K1 to K20) in the LEA package in R, and also based on the results of the DAPC analysis (see section Results). Prior to analyses, missing data were imputed using K = 3 based on 10 repetitions of the dataset. The Gibbs sampler was run five independent times on the dataset to estimate the appropriate average LFMM parameters (z-score values). p-values were adjusted using an estimation of the genomic inflation factor (λ) and a FDR of 0.05 to test for significance.

RDA, conditioned on site latitude, was performed using vegan (Oksanen et al., 2019) to identify candidate SNP loci associated with environmental variables. Significance of the full model and of each constrained axis was tested using an analysis of variance (ANOVA) and a marginal ANOVA, respectively with 999 permutations of the data. We used a ± 3 standard deviation cut-off from the mean SNP loading for each axis, as suggested by Forester et al. (2018) to optimize the number of correct positive detections while reducing false positives and identify SNPs significantly associated with environmental variables.

Finally, in order to investigate the functional attributes of identified outlier loci from LFMM and RDA, homology to the genome of the red spotted grouper (Epinephelus akaara, Temminck and Schlegel, 1842) (Ge et al., 2019) was assessed. RAD tags were aligned to the genome using minimap2 (Li, 2018), first directly to the gene.cds.fasta for direct hits to coding sequence (CDS). Remaining RAD tags were aligned to the genome.fasta and the annotation gff was used to guide hits to the nearest CDS. Resultant alignments for LFMM and RDA hits can be found in Supplementary Tables 4, 5.

Tagging and Telemetry of Nassau Grouper

During December 2014 – December 2016, Nassau grouper collected around the Exumas, Long Island and the Hail Mary FSA, were externally tagged (n = 103) using color-coded Floy^TM tags and fish of an appropriate size, i.e., ≥ 54 cm total length, TL (n = 45) were also surgically implanted with Vemco^TM V13 or V13P (13 mm × 36 mm, 6.5 g in water; Vemco, Ltd., Nova Scotia, Canada) acoustic transmitters and released as per procedures outlined by Stump et al. (2017). Specifically, 43 fish were acoustically tagged at the Hail Mary FSA, one fish was tagged in the Exuma Cays Land and Sea Park (ECLSP) and one individual in Long Island (Table 4 and Supplementary Table 6). Genetic analysis was also performed on a subset of acoustically tagged Nassau grouper (Tables 1, 4 and Supplementary Table 6). Eleven acoustic receivers (VR2Ws) stationed in the Exumas (n = 3), Long Island (n = 3), and at the Hail Mary FSA (n = 5) were used to record data on Nassau grouper migrations to and from the Hail Mary FSA over consecutive spawning seasons between December 2014 and January 2017 (Figures 1, 2).

TABLE 1

Table 1. Summary indices of genetic diversity per location based on 10,031 loci.

FIGURE 2

Figure 2. Map of acoustic receivers (VR2Ws) [denoted by green circles] deployed around the Exumas and Long Island, Bahamas. Receiver codes correspond to: SCHA, Shroud Cay; COCU, Conch Cut; GETO, Georgetown; HM, Hail Mary; NOPO, North Point; NECA, Newton’s Cay; POPO, Poseidon’s Point. SHCA and COCU are located within the boundaries no-take Exuma Cays Land and Sea Park (rectangular shape outlined in black). The inset map shows the Exuma Sound study area (black box) in relation to the rest of The Bahamas, Florida and the Caribbean.

Results

SNP Discovery and Genotyping

Sequencing of 96 Nassau grouper samples produced a total of 159,195,960 raw RAD-tags with an average of 1,658,291 (± 671,515 SD) RAD-tags across all samples (Supplementary Table 1). The number of RAD-tags generated for individual Nassau grouper ranged from 272,483 (LI239) to 3,653,394 (LI227). Following filtering and quality control, 114,174,309 RAD-tags (∼72%) were retained, with a mean of 1,189,316 (± 478,376 SD) RAD-tags across samples (Supplementary Table 1). After parameter optimization, a total of 5,139,297 assembled loci (min = 20,298; max = 63,594), 610,951 polymorphic loci (min = 1,382; max = 8,687), 869,675 polymorphic SNPs (min = 2,489; max = 12,139), and 78,615 fixed SNPs (min = 19; max = 8,216) were obtained (Supplementary Table 2). Average numbers of assembled loci, polymorphic loci, SNPs and fixed SNPs were 54,673 (± 8,317 SD), 6,499 (± 1,695 SD), 9,252 (± 2,173 SD), and 836 (± 1,253 SD) respectively (Supplementary Table 2). Approximately 35% of loci were polymorphic and contained an average of 1.3 SNPs per locus. After filtering in populations, the final dataset contained 28,636 assembled loci, of which 10,031 were polymorphic, containing a total of 13,241 SNPs.

Genetic Diversity, Differentiation, and Structure

A total of 262 loci deviated from HWE after FDR correction (p < 0.05). Measures of genetic diversity and allele polymorphism were generally consistent across sample locations (Table 1). Mean genome-wide diversity across all locations was H_E = 0.3627 (range 0.340–0.373; Table 1) and H_O = 0.332. H_O was lowest in Abaco (0.284) and highest in Exuma (0.363; Table 1). Inbreeding coefficients were highest in Abaco (F_IS = 0.127, G_IS = 0.165) and lowest in Exuma (F_IS = 0.024, G_IS 0.027; Table 1).

Global F_ST was low and non-significant (Global F_ST = 0.0003; p = 0.325). Genetic differentiation based on pairwise F_ST values ranged from -0.003 to 0.003 and was greatest between Exuma and Great Inagua (0.003), and Hail Mary and Ragged Island (0.003). Pairwise F_ST values were significant (p < 0.05) between only three population pairs (Supplementary Table 3).

DAPC analysis revealed three genetic clusters of Bahamian Nassau grouper (Figure 3). The analysis separated fish from Long Island and Exuma from the other locations (Figure 3). Based on the AMOVA results, most genetic variability (∼97%) occurs within individuals (Table 2). However, of the three AMOVA models, the third model, which separated Exuma and Long Island versus the rest of The Bahamas, explained most of the genetic variance (p = 0.055), compared with the other two models comparing the rest of The Bahamas and Exuma (p = 0.111) or the rest of The Bahamas and Long Island (p = 0.334; Table 2).

FIGURE 3

Figure 3. Discriminate analysis of principal components (DAPC) analysis of Nassau grouper based on 10,031 SNP loci illustrating separations of Exuma and Long Island from other sampled locations within The Bahamas. The DAPC was constructed using 40 principal components (PCs) and 3 discriminate functions. The percent variation explained by each axis appears in brackets.

TABLE 2

Table 2. Analysis of Molecular Variance (AMOVA) for the three models tested to examine genetic variation in 94 Nassau grouper sampled from eight islands and one spawning site in The Bahamas.

Identifying Loci Under Selection

Preliminary analysis showed that some of the environmental variables were highly correlated (Supplementary Figure 2). We therefore performed an initial RDA only on five environmental variables – maximum sea surface temperature, minimum sea surface temperature, minimum primary productivity, maximum current and salinity (SSTmax, SSTmin, PPmin, CURRmax and SAL). The resultant variance inflation factors (VIF) suggested collinearity between SSTmin and other variables. We therefore removed this variable and performed the final RDA on four variables (SSTmax, PPmin, CURRmax, and SAL). The first two axes explained 32 and 27.53% of the variation, respectively with a total of 16 SNPs associated with environmental variables – four with CURRmax, four with salinity, six with PPmin and two with SSTmax (Table 3 and Figure 4). However, results of this analysis should be treated with caution as ANOVAs showed both the full model and each axis to be non-significant.

TABLE 3

Table 3. Environmentally linked loci under putative selection based on latent factor mixed model (LFMM) and redundancy (RDA) analysis.

TABLE 4

Table 4. Details of acoustic telemetry data including transmitter identification (ID), tag date, detection information [i.e., dates, times, total number of detections across spawning seasons (shaded in progressively darker gray colors)] for Nassau grouper recorded on receivers in the Exuma Sound.

FIGURE 4

Figure 4. Redundancy analysis triplot of association between Nassau grouper single nucleotide polymorphism (SNP) frequency and four environmental variables: maximum current velocity (CURRmax), minimum primary productivity (PPmin), maximum sea surface temperature (SSTmax) and salinity (SAL) for eight sampling locations (AB, Abaco; AN, Andros; EL, Eleuthera; EX, Exuma; GI, Great Inagua; LI, Long Island; NP, New Providence; RI, Ragged Island). SNPs associated with each environmental variable are color coded (CURRmax – purple, PPmin – green, SSTmax – white and SAL – blue). The bottom and left axes relate to site and environmental variable scores while the top and right axes relate to SNPs scores.

LFMM analysis showed that seven loci were associated with two or more environmental variables (Table 3). A total of 67 (non-overlapping) hits associated with 58 genes was returned when the dataset was aligned with available red spotted grouper sequences for outlier SNP loci detected from LFMM analysis (51 loci) and RDA (16 loci) (Ge et al., 2019; Supplementary Tables 4, 5).

Tagging and Telemetry

Externally tagged Nassau grouper ranged in size from 21.2 to 73.0 cm TL, with an average of 51.97 (± 14.0 SD) cm TL (Supplementary Table 6). Local fishers reported only five tag returns from Nassau grouper. Acoustically tagged fish ranged in size from 53.1 to 73.0 cm TL, averaging 63.84 (± 4.2 cm TL) and 4.5 (± 1.1 SD) kg in weight (Supplementary Table 6). Hail Mary VR2Ws recorded detections from 43 Nassau grouper, 18 of which returned to the site during subsequent spawning seasons (Table 4). Of the 16 Nassau grouper tagged at Hail Mary in December 2014, one was detected between September and December 2015, two were detected 3 days after the November 25th, 2015 full moon, three were detected within a week around the December 25th, 2015 full moon, and one was detected 2 days before the January 4th, 2015 full moon (Table 4).

For Nassau grouper tagged in 2015, one individual that was tagged and released at Parrotfish Reef, ECLSP migrated to the Hail Mary FSA during the December 2015 spawning season, and another fish returned to the FSA for the week of the January 23rd, 2016 full moon (Table 4). One Nassau grouper, also tagged in December 2015, arrived at Hail Mary on the day of the full moon (December 13th, 2016) and remained within detection range at the FSA for approximately 4 h prior to its departure. Additionally, seven fish tagged in December 2016 were detected during the week of the January 12th, 2017 full moon (Table 4).

Five Nassau grouper that were acoustically tagged at the Hail Mary FSA in December 2016 were detected on two VR2Ws in the ECLSP a few days after the December 2016 full moon. One individual remained within detection range of Hail Mary VR2Ws for 5 days prior to departing for the ECLSP. Migrating Nassau grouper covered one-way distances of ∼137 km from Hail Mary to Conch Cut (COCU) in the south and ∼176 km to Shroud Cay (SHCA) in the north of the ECLSP during the 2016–2017 spawning season (Table 4). Three Nassau grouper that were tagged in December 2016 at Hail Mary were detected on the Georgetown (GETO) receiver in Exuma during the 2016–2017 spawning season (Table 4). Two of these fish were detected within 1–2 days after the December 13th, 2016 full moon and the other fish was detected on January 4th, 2017 (Table 4). Migrating Nassau grouper traveled ∼15.7 km (9.8 mi) from Hail Mary to GETO.

Nassau grouper were not detected on the Poseidon Point (POPO) receiver on the eastern side of Long Island. However, one individual was detected on the Newton Cay (NECA) VR2W for 27 consecutive days beginning on November 28th, 2015 and then again for another 4 days around the time of the January 2016 full moon (Table 4). Additionally, tag return information from fishers suggests that one Nassau grouper migrated from Hail Mary to Clarence Town, located in the southeastern end of Long Island.

Discussion

RAD-seq was used to discover and genotype SNPs for Nassau grouper, enabling assessments of genomic diversity and population divergence within The Bahamas. We also sought, for the first time, to investigate the role of selection on genetic population structure of Nassau grouper. Analyses of high-resolution SNP markers (including candidate loci under selection) coupled with acoustic telemetry revealed intraspecific population structure of Bahamian Nassau grouper, enabled us to determine the origins of individuals migrating to the Hail Mary FSA, and suggest that multiple processes influence the spatial dynamics of this species.

Genetic Diversity, Differentiation, and Structure

Global estimates of genetic differentiation (Global F_ST = 0.0003; p = 0.325) were lower than that reported by Jackson et al. (2014) using SNPs (F_ST = 0.002; p = 0.014) and microsatellite markers (F_ST = 0.002; p = 0.004), respectively. This, however, is not surprising given that Nassau grouper in the Jackson et al. (2014) study were sampled over an extensive geographical area, encompassing The Bahamas and parts of the Caribbean. Genome-wide assessments of diversity and differentiation with large numbers of SNPs have been shown to yield more accurate estimations of patterns of genetic variation (Helyar et al., 2011; Fischer et al., 2017; Hodel et al., 2017). Differences have been reported in gene diversity values between two studies for Nassau grouper using microsatellites and SNP markers (Jackson et al., 2014; Bernard et al., 2016). In the present study, however, findings from RAD-seq data concur with microsatellite data reported by Sherman et al. (2017), with both techniques demonstrating similar levels of genetic diversity in Bahamian Nassau grouper.

DAPC analysis is multivariate approach that attempts to simultaneously maximize variability between populations, while minimizing variation within populations. As a result, DAPC analysis has become an increasingly used method for revealing fine-scale patterns of population structure, especially when population differentiation estimates are low (Jombart et al., 2010; Kalinowski, 2010). Multivariate DAPC analysis enabled detection of three distinct genetic clusters (Figure 3), offering an apparently finer level of the resolution than that obtained by microsatellite analysis, which showed no structure (Sherman et al., 2017). The observed genomic divergences of Exuma and Long Island from the rest of The Bahamas (despite exhibiting similar levels of heterozygosity) suggest subtle differences in gene flow may be occurring among Nassau grouper from these islands.

In the absence of overt physical barriers to gene flow, probable explanations for the observed patterns of genetic connectivity include, (1) varied migration patterns and adult dispersal capabilities (Bolden, 2000; Dahlgren C. P. et al., 2016); (2) intraspecific differences in larval survivorship, dispersal, and recruitment (Eggleston, 1995; Kough and Paris, 2015); (3) low densities of Nassau grouper (0.13-0.25 fish/100 m²) outside the ECLSP (Dahlgren C. et al., 2016; Sherman et al., 2018b); and (4) divergent selective pressures between islands (e.g., Delmore et al., 2015; Moore et al., 2017).

As a transient spawner (Domeier, 2012), most of the reproductive effort for Nassau grouper is limited to a relatively narrow temporal window, which is further influenced by environmental and oceanographic processes (e.g., water temperature, currents, and the lunar cycle) (Colin, 1992, 2012). Through biophysical modeling, Kough and Paris (2015) demonstrated that infrequent spawners often exhibit reduced larval connectivity networks. Therefore, heavy exploitation of Nassau grouper outside the MPA coupled with reduced larval connectivity networks may potentially limit the relative contribution of gene flow throughout parts of The Bahamas. While residents of the ECLSP experience negligible fishing pressure compared to the rest of The Bahamas (Dahlgren C. et al., 2016), telemetry data in the present study, as well as previous research (Bolden, 2000; Dahlgren C. P. et al., 2016), support connectivity through adult migrations for spawning, which are variable in space and time.

Genetic connectivity in marine species is often complicated by geographical separation, migration, larval dispersal, life history characteristics (e.g., age at maturity, ontogeny, reproductive strategies), and oceanographic conditions (Rivera et al., 2011; Hemmer-Hansen et al., 2014; Selkoe et al., 2016). Future research is required to clarify patterns of gene flow and provide additional insights into the mechanisms influencing source-sink dynamics for Nassau grouper throughout The Bahamas. Such data will have important implications for maintaining genetic connectivity within Bahamian MPA networks, including the protection of specific FSAs, and in recommending fishery quotas. However, in the interim, stronger enforcement of national fisheries regulations particularly around the Exumas and Long Island are warranted to protect both genetic stocks.

Identifying Loci Under Selection

Examination of loci under selection has been employed to explain genetic population structuring in marine species (Overgaard Therkildsen et al., 2013; Benestan et al., 2015; Forester et al., 2018), which in turn is fundamental for understanding potential adaptive responses and associated relationships between genetic diversity. These data are therefore useful from both evolutionary and conservation management perspectives, providing the associated functions are known (Funk et al., 2012). To this end, we performed RDA of outlier SNPs to visualize weather patterns of genetic structure could be attributed to loci under selective processes (Figure 4 and Table 3). Additionally, LFMM analysis was used to detect correlations between SNPs and the environmental variables considered in the study (Supplementary Figure 2 and Table 3).

In the present study, pairwise F_ST with neutral loci showed no significant F_ST values (Supplementary Table 3), and there were no clear geographic patterns related to gene diversity for outlier loci linked to environmental variables (Figure 4). RDA results showed correlations between outlier SNPs and maximum sea surface temperature, maximum current, minimum primary productivity, and salinity, with maximum sea surface temperature and maximum current influencing the northern and southern extremes of The Bahamas (i.e., Abaco and Great Inagua, Figure 4). RDA and LFMM outlier SNPs possibly under selection were connected to a spectrum of important biological and molecular functions (Supplementary Tables 4, 5). Gene ontology indicates that loci of interest were related to chaperone and developmental proteins, oxidative stress, osmoregulation, metabolic processes, muscle differentiation, and function (Smithers et al., 2000; Chaube et al., 2017; Dardenne et al., 2017). Indeed, several of these gene families and transporter proteins (e.g., slc, ATPase, cpne, ptpn) have been shown to be connected to osmoregulatory and physiological processes in fish (DeFaveri et al., 2011; Tine, 2017).

In other research (e.g., Hess et al., 2013; DeFaveri and Merilä, 2014; Xu et al., 2017), loci under potential divergent selection were used to explain observed spatial patterns. For example, using SNP data obtained from RAD-seq, Hess et al. (2013) were able to link adaptive loci with environmental conditions that helped to explain demographic population structure in Pacific lamprey (Entosphenus tridentatus, Richardson, 1836). By examining patterns of gene flow, dispersal and putative loci under adaptive selection, Moore et al. (2017) found that Arctic char (S. alpinus) altered their migratory patterns to reduce exposure to unfavorable environmental conditions outside of the reproductive season.

Mean monthly water temperatures in The Bahamas range from 24°C in the winter to 29°C in the summer, but rarely drop below 26°C in the southern Bahamas during winter months¹. However, elevated water temperatures in nursery habitats (exposed to tidal fluctuations) and on the Bahama Banks are common during summer months, where maximum temperatures can exceed 40°C (Murchie et al., 2011; Shultz et al., 2014). Nassau grouper migrate during winter months, when water temperatures are cooler (Colin, 1992). Any extremes in water temperature are likely to be experienced by juvenile and sub-adult fish transitioning from nursery habitats to reefs (Eggleston et al., 1998). However, long distance migrations are likely to be more energetically costly, and thus may explain the presence of outlier loci associated with osmoregulation, metabolic processes and muscle function, which were also linked to environmental variables (Table 3 and Supplementary Tables 4, 5). There are, however, limitations with the different methods used for detecting loci under selection and even moderate gene flow may mask signals of divergence (Davey et al., 2011; Narum and Hess, 2011). Development and annotation of a genome for Nassau grouper would allow for rigorous exploration of genes associated with reproduction, physiology and other important biological functions, as well as correlations between fitness traits and adaptive variation.

Tagging and Telemetry

Round-trip migration distances covered by Nassau grouper in The Bahamas can exceed 300 km (Dahlgren C. P. et al., 2016), thus surpassing distances reported for fish in the Cayman Islands and Belize (e.g., Whaylen et al., 2004; Semmens et al., 2007; Starr et al., 2007). The migratory patterns displayed by Nassau grouper in this study are congruent with previous research in The Bahamas in terms of routes and distances (∼15.7–176 km) traveled (Bolden, 2000; Dahlgren C. P. et al., 2016; Stump et al., 2017). Adult Nassau grouper migrate along the shelf edge in groups to and from FSAs within 1–2 weeks of the full moon (Dahlgren C. P. et al., 2016; Stump et al., 2017). Moreover, Dahlgren C. P. et al. (2016) demonstrated that some fish from the ECLSP utilize FSAs around Long Island and highlighted temporal differences in migratory patterns, which coincided with the timing of the full moon. In the present study, of 45 acoustically tagged fish, 22% were detected by receivers traveling between Hail Mary and the Exumas (Figures 1, 2 and Table 4), and 14% entered or exited the ECLSP, a no-take marine protected area (MPA). Approximately 37% of Nassau grouper from the Exuma Cays were sampled within the ECLSP (Supplementary Table 5). Acoustic telemetry revealed this MPA as the likely origin of five Nassau grouper, representing 5 of 19 individuals (26%) tagged at the Hail Mary FSA in December 2016 (Supplementary Table 5). These results not only corroborate patterns of long-distance adult dispersal during the winter spawning season (Dahlgren C. P. et al., 2016), but also highlight the relative contribution of the ECLSP to the Hail Mary FSA.

However, the discrepancy between genetic and telemetry data suggest that fish may use other FSAs around Long Island or originate from other islands. Genetic data may not reflect this pattern due to low sample sizes (i.e., < 20) for several islands and limited spatial coverage across certain islands (e.g., Abaco, Grand Bahama, and Long Island). Future Nassau grouper sampling should aim to increase both geographic coverage and sample sizes. More detailed tagging data via satellite tagging or a more extensive network of acoustic receivers placed around migration routes, in tandem with additional genetic analysis of other potential reef sites of origin, may be necessary to fully understand spawning migration patterns around this part of The Bahamas. Identifying where all active FSAs are, the exact locations Nassau grouper migrate from to participate in annual spawning events, and patterns of larval dispersal and recruitment will be essential to the long-term protection and sustainable exploitation of this species.

Conclusion

Here, we aimed to determine whether a genomic approach, in conjunction with acoustic telemetry, could resolve intraspecific population structure and provide additional insight into the demographics of Nassau grouper in The Bahamas. RAD-seq analysis using more than 13,000 SNPs yielded new genomic data for the species. Key findings from this research show that genomic patterns of diversity in Nassau grouper exhibit a high degree of uniformity throughout The Bahamas. DAPC analyses revealed fine-scale genetic clustering – separating Exuma and Long Island from the rest of The Bahamas.

Patterns of cryptic divergence were further supported by AMOVAs, showing differences in the relative contribution of putative loci under selection in parts of The Bahamas. Acoustic telemetry data were used to confirm migratory corridors and track northwest movements of Nassau grouper from an active FSA into the Exumas and into the ECLSP. Long-distance movement patterns displayed by acoustically tagged fish provide additional evidence for a mechanism of genetic exchange, which is likely to be ubiquitous throughout The Bahamas. However, genomic analysis of Nassau grouper sampled at Hail Mary, did not reflect genetic signatures of fish from Exuma and Long Island, highlighting discrepancies between the genetic and telemetry data. Despite evidence of some tagged Nassau grouper moving to and from ECLSP and Hail Mary, genetic data suggest Hail Mary fish may also originate from other islands. MPAs like ECLSP that do not account for genetic connectivity and/or have not been explicitly designed to address fisheries management objectives are also unlikely to be completely effective conservation management tools. Collectively, these results demonstrate the critical importance of understanding the reproductive output of active FSAs and connectivity for well-designed MPAs and MPA networks, as they can be valuable tools for species conservation.

Using RAD-seq, new information on putative loci associated with divergent selection in Nassau grouper has been derived. It appears that selection may also be an important driver in shaping the contemporary genetic population structure of the species and these data certainly warrant further exploration to identify potential processes influencing selection. These findings may have implications for long-term survival in light of future anthropogenic disturbances (e.g., forecasted increases to water temperatures, which may impact the timing, and energetic costs of spawning migrations). Moreover, the use of SNP data will facilitate regional comparisons of genetic diversity, differentiation and structure once such data become available from other countries where Nassau grouper populations still exist. In the future, Nassau grouper SNP data may have applications beyond demographic inferences, such as examining the extent of illegal unreported and unregulated (IUU) fishing for fisheries management (e.g., Martinsohn and Ogden, 2009) and population assignment (e.g., Larson et al., 2014).

Overall, data from this study build upon previous Nassau grouper research in The Bahamas. It demonstrates the power of RAD-seq to reveal complex spatial patterns and provide novel insight into population dynamics of a critically endangered marine species. Furthermore, the application of both in situ methodologies (i.e., acoustic telemetry) and molecular research highlight the benefits and importance of integrative approaches for advancing understanding the dynamics of aggregating species.

Data Availability Statement

RAD-seq data are available in the European Nucleotide Archive (ENA) via study accession number PRJEB36904. Sample accession numbers as ERX3958594–ERX3958689.

Ethics Statement

Established ethical procedures were reviewed by Shedd Aquarium’s Internal and External Research Committee and were adhered to during animal handling.

Author Contributions

KDS conducted fieldwork, performed genomic DNA extractions and quality testing, prepared samples for RAD-tag library development, analyzed the data, and wrote the manuscript. JP assisted with bioinformatics and contributed to the manuscript. RK provided laboratory support, performed RDA analysis, and contributed to the manuscript. KM developed RAD-tag libraries and completed sequencing. CD, LK, and KS were involved with fieldwork and commented on the manuscript. LK also produced Figures 1–2. JS and CT helped in project design and contributed to improving the manuscript. The final version was approved by all authors.

Funding

RAD-sequencing was funded by the Save Our Seas Foundation (small Grant No. 365) and Lyford Cay Foundation via the Shirley Oakes Butler Charitable Trust scholarship to KDS. Fieldwork was funded by Shedd Aquarium. Open access publication fees were covered by the University of Exeter.

Conflict of Interest

KS was employed by company Field Lab Consulting, LLC.

The remaining 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.

Acknowledgments

We would like to thank the Exeter Sequencing Service and Computational core facilities team at the University of Exeter. Fin clips were obtained with permission from The Bahamas Department of Marine Resources under scientific research permit (MAMR/FIS/17 – Sherman) to KDS and follow UK Home Office guidelines. Fieldwork and telemetry research were completed under permits issued to KDS (MAMR/FIS/17 – Sherman) and KS (MAMR/FIS/17 – Stump). Finally, we are grateful to research volunteers from the Shedd Aquarium, Bahamas National Trust, Atlantis, Cape Eleuthera Institute, Bahamas Reef Environment Educational Foundation, The Nature Conservancy, University of Miami, T. Kamerman, S. Buhler-Hensel, D. Rottman, D. Delaney, N. Truelove, B. Crown, and R. Fluck who assisted with fieldwork between December 2014–May 2017 and to Keith Pamper and the crew of the R/V Coral Reef II for logistical field support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.00157/full#supplementary-material

Footnotes

^ https://www.seatemperature.org/central-america/bahamas/nassau.htm

References

Abelson, A., Obolski, U., Regoniel, P., and Hadany, L. (2016). Restocking herbivorous fish populations as a social-ecological restoration tool in coral reefs. Front. Mar. Sci. 3:138. doi: 10.3389/fmars.2016.00138