Juveniles and adults were targeted for DNA extraction; postlarvae were excluded to ensure individuals collected had survived their initial migratory journey. Juveniles were expected in nursery areas and adults on spawning grounds. Only adults would be present on the spawning grounds as they return to spawn. DNA was extracted from the plucked abdominal muscle tissue with the DNeasy Blood and Tissue kit (Qiagen), following the manufacturer’s instructions. To ensure a sufficient amount of DNA had been obtained from an extraction for downstream ddRADseq library prep, DNA was quantified with the Qubit dsDNA High Sensitivity Analysis kit (ThermoFisher). Gel electrophoresis was used to confirm the presence of intact, high molecular weight DNA: DNA extractions were run through a 2% agarose gel for 90 min at 100 V, visualized with GelRed (Biotium). Only samples with more than 500 ng of unfragmented DNA were considered for ddRADseq library prep.

Of the 105 F. duorarum specimens that underwent DNA extraction, a subset were found to meet the criteria described above. Of these, 68 were chosen for next-generation sequencing library prep (∼10 samples per sampled region). Reduced representation libraries were prepared following the double digest Restriction-site Associated DNA sequencing (ddRADseq) method published by Peterson et al. (2012). Briefly, we began with a series of enzyme trials to determine the optimal enzyme combination and size selection range to provide adequate genomic coverage at adequate sequencing depth. At least 500 ng of extracted DNA was digested with SphI-HF and EcoRI-HF (New England Biolabs) for 3 h at 37°C. Enzymatic activity was stopped with a 30 min hold at 65°C. Custom barcode adapters (synthesized at Integrated DNA Technologies) were ligated to the double-digested fragments using T4 ligase (New England Biolabs). Following barcode adapter ligation, samples were pooled into nine samples of eight, uniquely barcoded libraries. Fragments between 270 and 330 bp, including adapter length, were size selected on a PippinPrep with a 1.5% Agarose Gel Cassette (Sage Science). To reduce the impact of PCR bias, each size-selected sample was subdivided into five parallel PCR amplification reactions and a negative control was used to ensure reagents were not contaminated. Using the Phusion Hi-Fidelity Polymerase (Thermo Scientific), the PCR reactions went for 10 cycles and incorporated i7 indices and Illumina adapters into every amplified fragment, allowing for pooling of all libraries into a single sample. This final sample was quality-checked on an Agilent BioAnalyzer 2100 (Agilent Technologies) immediately prior to sending it for sequencing with the Illumina HiSeq4000, SE150, at the University of Texas’ Genomic Sequencing and Analysis Facility.

Initial quality checks of the raw data were conducted with fastQC (Andrews, 2010) before data assembly began in STACKS v1.45 (Catchen et al., 2013) on Florida International University’s High Performance Computing Cluster (FIU HPCC). Given the risk of data assembly decisions resulting in a biased data set, recent literature was consulted before beginning the complex task of generating datasets from ddRADseq data (Mastretta-Yanes et al., 2015; Paris et al., 2017; Rochette and Catchen, 2017; O’Leary et al., 2018). Data assembly followed the recommended core pipeline for de novo data: process_radtags to demultiplex the reads, ustacks to align reads within each individual, cstacks to catalog these reads, sstacks to query putative loci against this catalog, and rxstacks to utilize population data to correct individual genotype calls. As any individual dataset, assembled according to the authors’ best judgment, can reflect biases stemming from assembly decisions, nine datasets were generated, differing in the maximum Hamming distance allowed between stacks (ustacks’ –M), the minimum depth required to designate a stack (ustacks’ –m), and the maximum Hamming distance allowed between sample loci (cstacks’ –n). The data assembly parameters for each dataset are presented in Table 1. These datasets are referred to as “batches” and reflect the parameter settings that generated them: “batch161” is the dataset assembled with a maximum Hamming distance of 1 allowed between stacks, a minimum stack depth of 6, and a maximum Hamming distance of 1 allowed between sample loci (–M 1 –m 6 –n 1).

Quality filtering of the VCFs output from STACKS was accomplished with VCFtools (Danecek et al., 2011) on the FIU HPCC. First, the minimum read depth was set to 10×. Next, sites with ≥50% missing data were removed, followed by individuals with ≥90% missing data. The resulting VCF files were reformatted in PGDSpider v2.0.5.2 (Lischer and Excoffier, 2012) for analysis in BayeScan v2.1 (Foll and Gaggiotti, 2008; Foll et al., 2010; Fischer et al., 2011), which identifies loci which may be under natural selection, as well as GenAlEx v6.501 (Peakall and Smouse, 2006, 2012).

Pairwise measures between regions, including Nei’s unbiased genetic distances, which describe allelic differences assuming genetic drift and mutation are in equilibrium (Nei, 1972, 1987, and F_ST values, which quantifies the proportion of genetic variation explained by population structure (Wright, 1950), were calculated in the Excel data analysis suite, GenAlEx v6.501. GenAlEx was also used to identify private alleles within each region and conduct the Analyses of Molecular Variance (AMOVAs). The number of private alleles identified for every region were normalized by each region’s sample size (PA_norm). Pairwise F_ST values were calculated alongside the AMOVAs utilizing GenAlEx’s “AMOVA” option. Standard permutation was selected to calculate statistical significance of results over 999 permutations. Missing data were not imputed. Neighbor Joining (NJ) trees and Principal Component Analyses (PCAs) were constructed in the R package, adegenet (Jombart, 2008; Jombart and Ahmed, 2011). Three principal components (PCs) were calculated for each dataset, plotting the primary and secondary PCs with ggplot2 (Wickham, 2016). Ellipses, encompassing the 0.95 confidence levels, were added for each region. Finally, using the “dapc” command in adegenet, Discriminant Analyses of Principal Components (DAPCs) were built from the first three PCs for each dataset.

Population structure was tested in STRUCTURE v2.3.4 (Pritchard et al., 2000) with K taking values between 2 and 7, each tested 10× under the admixture model with allele frequencies correlated among populations. Initially, each analysis ran for 100,000 generations, and the first 25% were discarded as burn-in. Review of preliminary results found high agreement between replicates, indicating that this number of generations was sufficient to achieve convergence. After STRUCTURE analyses were complete, results were collated in STRUCTURE HARVESTER v0.6.94 (Earl and VonHoldt, 2012). Within STRUCTURE HARVESTER, the optimal K value was inferred using ad hoc posterior probability models (Pritchard et al., 2000) and the Evanno Method (Evanno et al., 2005). STRUCTURE plots were generated within the R package pophelper (Francis, 2017).

While Criales et al. (2015) presented a suite of models, for simplicity, here we focus only on the model that incorporates the larval behaviors of DVM and STST and describes the major and minor routes, as this was the only model that resulted in successful recruitment. Testing the hypotheses indicated by the modeling work of Criales et al. (2015) may be accomplished through patterns of unique haplotypes (private alleles), measures of genetic distance (Nei’s unbiased distance), population differentiation (F_ST), and components of genetic variance (AMOVA). It is important to note that statistical tests are performed on values calculated from pseudoreplicated datasets (batches), not fully independent data.

Expectations under the null hypothesis: Most of the survivorship research on F. duorarum has focused on recruitment success (Browder et al., 1999, 2002; Ehrhardt and Legault, 1999; Criales et al., 2006, 2007, 2015), describing a density-dependent trend (Ehrhardt et al., 2001). By definition, the spawning aggregation represents the highest population density of sexually mature, spawning shrimp. Adults found on nearshore nursery grounds have matured on those grounds or in nearby estuaries and will soon return to spawning grounds for their turn at spawning. Spawning aggregations hold greater genetic diversity than found on any one nursery ground when spawners come from several nursery locations. Under the null hypothesis, we expect the highest number of private alleles-per-individual (PA_norm) to come from sites representing a spawning aggregation. A t-test, assuming unequal variance, was utilized to statistically compare PA_norm for spawning (DT and MQ) vs. nursery (NBB, BB, SBB, EVG, and NEVG) regions. Most estuaries from which samples were collected for this study represent nursery areas, although young shrimp may move out of an estuary to avoid disruptive changes in conditions such as storms or cold snaps (e.g., see Tabb et al., 1962, pp. 26–27).

Finally, under the null, we expect little-to-no statistically significant pairwise population differentiation between regions; the vast majority of genetic variance should come from differences between individuals (F_IT). Pairwise F_ST values and AMOVA results will provide support in this regard.

Expectations under the alternative hypothesis: If the Dry Tortugas and the Marquesas support population-specific spawning aggregations, we expect statistically significant pairwise population differentiation between these sites, which was tested with an ANOVA comparing pairwise F_ST values by region type: spawning-spawning (DryTortugas-Marquesas), spawning-nursery (all region pairs containing DryTortugas or Marquesas), and nursery-nursery (all region pairs that do not contain DryTortugas or Marquesas). Moreover, while the majority of molecular variance may be attributable to variance among individuals (F_IT), F_ST should be greater than zero and statistically significant.

In additional to the statistical tests described, PCAs, DAPCs, and STRUCTURE results were evaluated for evidence of population structure. Any results, quantitative or qualitative, contradictory to both hypotheses will be considered as contradictions to the validity of the model presented by Criales et al. (2015), and the relative strength of such contradictions will be assessed in the context of the full study presented here.

Nine datasets were analyzed to better understand the robustness of results to data assembly decisions, however, results across datasets were highly similar. While all results are reported in the Supplementary Materials, only results from batch161, the dataset with the highest number of samples and SNPs (N = 57, SNPs = 799) are presented in-text. Because very few samples representing the SouthBiscayneBay region were retained following quality filtering (n = 3), these samples were removed for calculation of Nei’s unbiased distance calculation, estimation of pairwise F_ST between populations, and AMOVAs. The SouthBiscayneBay samples were included in PA_norm, PCAs, NJ trees, and STRUCTURE analyses.

Estimates of Nei’s unbiased genetics distance between all region-pairs, excluding SouthBiscayneBay, ranged from 0.003 to 0.006, with the highest value attributable to the comparison between BiscayneBay and North_of_Everglades (Figure 2A and Supplementary Figure 1). The shortest genetic distance was calculated for multiple region-pairs: North_of_BiscayneBay compared to either spawning region (DryTortugas and Marquesas), Everglades compared to either spawning region, and North_of_BiscayneBay compared to Everglades. All other region-pair distances fell between 0.004 and 0.005 (Supplementary Table 2).

Pairwise comparisons between regions, excluding SouthBiscayneBay, were also examined through estimates of F_ST (Figure 2B and Supplementary Figure 2), which ranged from 0.000 to 0.102. Many region-pairs returned null F_ST values: all comparisons including BiscayneBay or North_of_Everglades and a region associated with the nursery range of F. duorarum (Everglades and North_of_BiscayneBay), as well as the Marquesas-DryTortugas region pair. With the exception of the North_of_BiscayneBay-Marquesas region pair, all non-zero F_ST were characteristic of region pairs that included a spawning region, with the highest F_ST values calculated between DryTortugas and Everglades (Supplementary Table 3), though recall that DryTortugas-Everglades had a very low genetic distance.

Analyses of Molecular Variance across all nine datasets, excluding SouthBiscayneBay, yielded an average among-population variance value of 1.69% (standard deviation 0.71%, Table 3). The vast majority of molecular variance was attributable to differences among individuals (88.49%, standard deviation 1.68%) and the remainder came from differences within individuals. Overall average F_ST (the proportion of total genetic variance found within a population), F_IS (the proportion of genetic variance in a population which is found within an individual from that population), and F_IT (the proportion of total genetic variance found within an individual) reflected these values as well (0.017 ± 0.007, 0.900 ± 0.014, and 0.902 ± 0.014, respectively). Across the nine AMOVAs, F_ST was statistically significant in two cases, while F_IS and F_IT were statistically significant in every case.

Results from Principal Component Analysis (batch161 presented in Figure 3A; all batches presented in Supplementary Figure 3) and DAPCs (batch161 presented in Figure 3B, all batches presented in Supplementary Figure 4) include all samples, revealing a large, central cluster. However, samples from the DryTortugas are slightly shifted from the center. The NJ results (batch161 presented in Figure 3C; all batches presented in Supplementary Figure 5), which included samples from SouthBiscayneBay, show little structure. Across NJ trees, only two BiscayneBay samples are differentiated from the otherwise unstructured tree, but the other BiscayneBay samples do not reflect a larger separation of the region from the rest of the samples. To ensure these individuals did not represent contamination, we confirmed the taxonomic identification of these two samples as F. duorarum.

Further examining relationships between samples with the K-means clustering method STRUCTURE, the Evanno method was applied to identify the optimal K in each analysis. Across datasets, the Evanno method identified K = 2 as the optimal number of clusters within the data (Supplementary Figure 6). The two BiscayneBay samples differentiated in the NJ trees are clearly seen in the STRUCTURE plots as representing higher proportions of the minor cluster, otherwise all individuals appear highly similar, regardless of collection region (Figure 3D).

Normalized counts of private alleles within each region (PA_norm), including SouthBiscayneBay, ranged from 7.3 (standard deviation 0.7, Everglades) to 10.4 (standard deviation 1.1, BiscayneBay) private alleles per sampled individuals (Figure 4 and Supplementary Table 4). With the exception of BiscayneBay, spawning regions had higher normalized private allele counts (Marquesas = 10.2 ± 1.2, DryTortugas = 10.1 ± 0.6) than regions from the nursery range.

Expectations under the existing model were evaluated through several tests of significance (Table 4): to begin, we evaluated whether PA_norm differed significantly between spawning regions and nursery regions. A one-tailed, two-sample t-test assuming unequal variances between PA_norm of nursery regions (North_of_BiscayneBay, BiscayneBay, SouthBiscayneBay, Everglades, and North_of_Everglades) and spawning regions (DryTortugas and Marquesas) indicated significantly higher PA_norm in spawning regions (t_stat = –4.46; p = 2.23 × 10^–5). Next, we performed two single-factor ANOVAs to test whether Nei’s unbiased genetic distances or pairwise F_ST values differed significantly between types of region-pairs: spawning-spawning, spawning-nursery, and nursery-nursery (Table 4). The ANOVA analyzing Nei’s unbiased distances between region-pairs did not detect a statistically significant difference between region-pair types (F_stat = 1.95; p = 0.15). The ANOVA analyzing pairwise F_ST values, however, yielded a statistically significant result (F_stat = 42.83; p = 4.63 × 10^–15). We followed the ANOVAs with three two-tailed, paired t-tests comparing PA_norm, Nei’s unbiased genetic distances, and pairwise F_ST between all nursery regions and DryTortugas to all nursery regions and Marquesas. The normalized number of private alleles and Nei’s distances did not differ significantly by spawning region (PA_norm t_stat = 2.31, p = 0.91; Nei’s t_stat = 0.48, p = 0.63), while pairwise F_ST values were significantly higher between region-pairs including DryTortugas compared to region-pairs including Marquesas (t_stat = –8.22; p = 1.09 × 10^–9).

There is no paucity of potentially confounding variables when modeling current- and tide-mediated transport of dispersing larvae: the oversimplification of active swimming behaviors and the disparity between potential and realized dispersal has been described previously, including the biological importance of single individuals occasionally dispersing long distances (Shanks, 2009). However, biophysical modeling can be used in concert with genetic evidence to improve our understanding of the dynamic relationships between marine organisms and their environment (Liggins et al., 2013; Timm et al., 2020). Such an integrative approach has been utilized in studies of marine invertebrate populations (Dawson et al., 2005), including a recent study investigating the causes of population structure in an economically important decapod, the spiny lobster Panulirus argus (Truelove et al., 2017).

The biophysical oceanographic model developed by Criales et al. (2015) describes two migratory routes, which differ in their origin (Dry Tortugas and Marquesas vs. Dry Tortugas only), usage (many vs. few individuals, represented as particles), and recruitment success (majority of particles are successfully recruited to Florida Bay vs. few particles are successfully recruited). These differences have the potential to maintain intraspecific diversity via population differentiation. It is important to note that no model perfectly reflects reality; while the model developed by Criales et al. (2015) accounts for direction and velocity across water depth, this information is not discussed in the work. However, the model provides three questions that can be addressed with population genomics: Is there independent support for the model? Do the modeled spawning aggregations sufficiently explain the genomic results? Do we see evidence that the minor route sustains a differentiated population?

Next-generation sequencing data provided strong support for the existing model of larval transport: across analyses, samples collected from the Marquesas and Dry Tortugas were clearly part of a larger population present across the Florida Peninsula (see Figure 3). The presence of significantly more private alleles in the spawning regions compared to the nursery sites (Table 4) further supports the model of larvae originating from the Dry Tortugas and the Marquesas. It is worth explicitly addressing the two Biscayne Bay outliers identified throughout clustering analyses in Figure 3, which suggest recruitment to Biscayne Bay from a spawning aggregation that was not sampled in this study. In this regard, the existing model, which simulates spawning aggregations in the Marquesas and the Dry Tortugas, may not be complete and an additional spawning site contributes recruits to the region.

Under the null hypothesis, we expect one homogeneous population present throughout the study region. While cluster analyses (PCA, DAPC, and STRUCTURE) do not clearly delineate populations, we see some shifting of samples from the Dry Tortugas (Figure 3), and statistical tests of population differentiation (global and pairwise) indicate low levels of structure throughout (though these values are only rarely statistically significant). This structure provides evidence for the alternative hypothesis: the separate spawning aggregations and migratory routes (major and minor) support genetic structure in the pink shrimp population around the Florida Peninsula.

With few exceptions, significant pairwise population differentiation was highest and statistically significant when regions from the nursery range were compared to spawning regions (Figure 2, Supplementary Figure 2, and Supplementary Table 3). Examining pairwise population differentiation by region pair type (spawning-spawning vs. spawning-nursery vs. nursery-nursery) revealed significant differences (Table 4), with differentiation between the spawning-spawning pair < nursery-nursery pair < spawning-nursery pair. To a large extent, the Dry Tortugas seems to be driving this trend: analyses of population differentiation indicate the Marquesas region is better genetically connected to the nursery regions than the Dry Tortugas region is (Figure 2B and Supplementary Table 3). It should be noted that the highest significant pairwise population differentiation calculated in this study was relatively low, but low, statistically significant F_ST estimates are fairly common in the marine realm (Waples, 1998; Waples and Gaggiotti, 2006; Hauser and Carvalho, 2008; Therkildsen et al., 2013; Timm et al., 2020). This would also explain the lack of clear structuring in clustering analyses.

The presence of a differentiated population in the Dry Tortugas (hereafter referred to as the “Dry Tortugas subpopulation”) is unexpected. Recall that larvae are spawned offshore of the Dry Tortugas and the Marquesas. Larvae pass through a series of developmental stages as they migrate, taking the major or minor route (Figure 1), to estuarine nursery grounds around the Florida Peninsula where they complete their maturation into adults. Year-round, these adults migrate back to the spawning aggregations to reproduce, which should, theoretically, lead to sufficient mixing to result in a single, genetically homogeneous population. We suspect the maintenance of a Dry Tortugas subpopulation may be the result of geographic or temporal separation of spawning populations. By the geographic mechanism, the Dry Tortugas subpopulation spawns exclusively in the Dry Tortugas and solely utilizes the minor migratory route, while the larger population spawns in the Dry Tortugas and the Marquesas and utilizes the major migratory route. However, the lack of clearly defined genetic structure separating the Dry Tortugas subpopulation from the larger population suggests this geographic mechanism is not sufficient to explain the results presented here.

The population structure we identify may also be the result of a temporal mechanism: since the 1980s, Key West shrimpers have reported anecdotal evidence of two spawning surges annually for the past several decades (pers. comm.) and unpublished data of Robblee suggest two peaks in population abundance of juvenile pink shrimp in western Florida Bay (pers. comm.). Costello and Allen (1966) also remark on the seasonal nature of juvenile pink shrimp in the region. Without additional data, it is difficult to characterize this mechanism further; however, if adults of the Dry Tortugas subpopulation arrive at the Dry Tortugas spawning ground before or after the larger aggregation, they will only be able to reproduce with each other. Moreover, depending on the seasonal timing of this second spawning surge, larvae originating from the Dry Tortugas subpopulation may utilize the minor migratory route to Florida Bay, leading to higher mortality and lower recruitment success. Given the lack of a clearly distinguishable Dry Tortugas subpopulation in the clustering analyses, it may be that such a mechanism results in the differentiation of the Dry Tortugas subpopulation, with occasional gene flow between it and the larger population preventing strong genetic structuring.

Either mechanism, geographic or temporal, might be facilitated by local recruitment of the Dry Tortugas population to the Dry Tortugas or a region not represented by samples collected for this study. In line with the alternative hypothesis, we find evidence of an unsampled spawning aggregation contributing individuals to Biscayne Bay and the Everglades: both regions show low-but-significant differentiation from the Marquesas and the Dry Tortugas, but no differentiation between themselves. The Loop Current’s episodic influence may bring migrants into nearshore currents, bringing recruits to Biscayne Bay from the Caribbean (Saloman et al., 1968). Alternatively, migrants may be contributed from the Sanibel spawning aggregation. A previous mark-recapture study found that, while geographic ranges of stocks from the Dry Tortugas and Sanibel overlap in nursery grounds, there is only evidence of weak, one-way migration of Sanibel stocks to the Dry Tortugas (Costello and Allen, 1966). Such separation between spawning grounds could provide a basis for population differentiation. Interestingly, this mark-recapture study did not find any evidence of shrimp migration between Biscayne Bay and the Sanibel grounds, nor between Biscayne Bay and the Dry Tortugas; indeed, all individuals marked and released within Biscayne Bay were only ever recovered from Biscayne Bay. The results presented here contradict this study, finding gene flow between the Dry Tortugas and Biscayne Bay (though Biscayne Bay may also receive recruits from a spawning aggregation that was not sampled in the current study).

Spawning-recruitment relationships of pink shrimp in south Florida appear to be more complex than previously believed and additional research is needed to investigate the mechanisms we hypothesize here. Representative sampling of F. duorarum from Sanibel, Cuba, and the Bahamas would be needed to further investigate the relative support for these potential sources of postlarval migrants. Anecdotal evidence of spawning surges, and the role this may play in the population structure of pink shrimp around the Florida Peninsula, would require a longitudinal study to better understand this mechanism.

The fisheries supported by F. duorarum contribute to economies internationally (Sheridan, 1996; Ramírez-Rodríguez et al., 2003; Hart et al., 2012), and the continued exploitation of this natural resource is critically dependent on the stability and sustainability of the species in the Gulf of Mexico and around the Florida Peninsula. One crucial factor contributing to species stability is successful larval recruitment: the movement of larval and postlarval individuals from spawning aggregations to nursery grounds.

Our results support the biophysical oceanographic model developed by Criales et al. (2015), which indicates a major route, traversed by larvae from the Dry Tortugas and the Marquesas, and a minor route, which only resulted in successful recruitment when larvae originated from the Dry Tortugas. Moreover, we find evidence that samples from the Dry Tortugas represent a differentiated population. Co-located with this region is a pink shrimp fishery (Klima et al., 1987; Hart et al., 2012), which harvests mature and young adult shrimp year-round on the lower southwest Florida shelf (Ehrhardt and Legault, 1999; Browder et al., 2002), perhaps with important implications for intraspecific genetic diversity: individuals harvested near the Dry Tortugas may represent the subpopulation indicated by our analyses. The removal of these individuals could undermine the subpopulation’s stability by reducing the density of juveniles and subsequently decreasing recruitment success (Ehrhardt et al., 2001).

Additional work is needed to further characterize the role of these two spawning grounds and migration routes, particularly by including individuals collected from the Sanibel grounds and the Caribbean. Such research will assist in determining whether the species should be managed as a single stock or if more complex management is required. Enhancing our understanding of larval recruitment success in F. duorarum will ultimately improve the long-term sustainability of these fisheries while protecting diversity within the species.