Population Genomics of the Commercially Important Gulf of Mexico Pink Shrimp Farfantepenaeus duorarum (Burkenroad, 1939) Support Models of Juvenile Transport Around the Florida Peninsula

The Gulf of Mexico pink shrimp, Farfantepenaeus duorarum, supports large fisheries in the United States and Mexico, with nearly 7,000 tons harvested from the region in 2016. Given the commercial importance of this species, management is critical: in 1997, the southern Gulf of Mexico pink shrimp fishery was declared collapsed and mitigation strategies went into effect, with recovery efforts lasting over a decade. Fisheries management can be informed and improved through a better understanding of how factors associated with early life history impact genetic diversity and population structure in the recruited population. Farfantepenaeus duorarum are short-lived, but highly fecund, and display high variability in recruitment patterns. To date, modeling the impacts of ecological, physical, and behavioral factors on juvenile settlement has focused on recruitment of larval individuals of F. duorarum to nursery grounds in Florida Bay. Here, we articulate testable hypotheses stemming from a recent model of larval transport and evaluate support for each with a population genomics approach, generating reduced representation library sequencing data for F. duorarum collected from seven regions around the Florida Peninsula. Our research represents the first and most molecular data-rich study of population structure in F. duorarum in the Gulf and reveals evidence of a differentiated population in the Dry Tortugas. Our approach largely validates a model of larval transport, allowing us to make management-informative inferences about the impacts of spawning location and recruitment patterns on intraspecific genetic diversity. Such inferences improve our understanding of the roles of non-genetic factors in generating and maintaining genetic diversity in a commercially important penaeid shrimp species.


INTRODUCTION
The Gulf pink shrimp, Farfantepenaeus duorarum (Burkenroad, 1939) supports multiple, international fisheries along its described geographic range, representing millions of dollars of economic activity (Sheridan, 1996;Ramírez-Rodríguez et al., 2003;Hart et al., 2012). Over 7,000 tons of pink shrimp were harvested across fisheries in the Gulf of Mexico in 2016, the last year for which such data are available (Hart, 2017). Given the economic and social influence of the large-scale fishing effort directed at Farfantepenaeid species in the Gulf, proper management is critical to the sustained stability of the species and protection of economic interests in the region: all of the species within the Farfantepenaeus group are targeted by fisheries to some extent (see Timm et al., 2019 for more information).
Management of fished species requires understanding the biology and ecology of the organism, including assessments of intraspecific biodiversity and the evolutionary processes that drive it (Bernatchez, 1995). Management of F. duorarum by Mexico and the United States of America makes such insight particularly crucial: shrimp fisheries have supported regional Mexican economies for decades, and pink shrimp have contributed substantially to these fisheries, with 90% of fished shrimp in the 1990s being F. duorarum (Arreguín-Sánchez et al., 2008). In the late 1990s, however, the F. duorarum fishery in the southern Gulf of Mexico was declared collapsed (Arreguín-Sánchez et al., 1997). Investigation of possible underlying causes of the collapse found evidence for decreased stockrecruitment (Arreguín-Sánchez et al., 1997, and efforts were undertaken to promote recovery (Arreguín-Sánchez et al., 2008). Such events have occurred in United States fisheries as well, resulting in the closure of the northern brown shrimp (F. aztecus) fishery along the Texas coast in the 1980s (Klima et al., 1987). The co-occurrence of several, economically important species of Farfantepenaeus along the coasts of the Gulf of Mexico further complicate management. Specifically, juvenile individuals of F. brasiliensis and F. duorarum look very similar, and the ability to confidently identify juvenile individuals taxonomically by reproductive structure morphology (Pérez-Farfante, 1988) is nearly impossible (Ditty and Alvarado Bremer, 2011;Teodoro et al., 2016). A recent study found cryptic diversity within F. brasiliensis, identifying two distinct populations (one occupying United States coasts and the other present along the coasts of South America). The study called for additional efforts to better understand population structure and evolutionary history within managed species . A break in species composition exists between the Gulf of Mexico and the greater Atlantic; divided by prevailing environmental features (Avise, 1992;Young et al., 2002). Studies focused on genetic diversity and population connectivity in species that span this break (such as F. duorarum) might prove particularly informative.
Life history can be significant in determining the composition and structure of adult assemblages, especially in species with complicated development cycles. Adults of F. duorarum spawn year-round in aggregations offshore of the Dry Tortugas and the Marquesas on the southwest Florida shelf (Cummings, 1961;Roberts, 1986). There is a distinct spawning aggregation on the Sanibel grounds as well, and, despite geographic overlap between Sanibel and Dry Tortugas nursery grounds, a division between shrimp originating from these two spawning grounds has been noted near Indian Key (i.e., between Sanibel and Dry Tortugas; Costello and Allen, 1966;Robblee et al., 1999): shrimp emanating from Sanibel nursery grounds only rarely migrate into the Dry Tortugas trawling grounds south of Indian Key and vice versa. After hatching, larvae rapidly progress through 11 developmental stages [nauplii (5), protozoea (3), and mysis (3)] in approximately 15 days (Dobkin, 1961). During this time, larval individuals exhibit a vertical migration pattern, alternating between deeper waters and surface waters (Rothlisberg, 1982;Rothlisberg et al., 1995Rothlisberg et al., , 1996Condie et al., 1999). For the first 15 days of development, vertical migration is triggered by light [diel vertical migration (DVM)]. During the subsequent 15 days, as individuals pass through several postlarval stages (3-6 stages; Ewald, 1965), vertical migration is timed to tidal movement [selective tidal-stream transport (STST)], allowing postlarvae to take advantage of tidal movement toward nursery grounds and avoid tidal movement in the opposite direction (Forward and Tankersley, 2001;Queiroga and Blanton, 2005). A modeling study by Criales et al. (2015) suggests that these two behaviors, DVM and STST, facilitate movement from spawning grounds toward primary nursery grounds in Florida Bay and mangrove estuaries along the southwest coast (Tabb et al., 1962;Costello and Allen, 1966;Browder and Robblee, 2009), where they grow through the juvenile stage, returning to spawning grounds as young adults. Environmental factors such as salinity and temperature on their nursery grounds affect their rate of growth and mortality (Browder et al., 1999(Browder et al., , 2002Ehrhardt and Legault, 1999), potentially influencing recruitment to the offshore fishery.
Two routes have been proposed for larval/postlarval migration (Figure 1): larvae may drift east and northeast along the Florida Current, to enter Florida Bay through the Florida Keys (Munro et al., 1968;Criales et al., 2003). The other route posits that larvae move more directly across the southwest Florida shelf, entering Florida Bay at its northwest side (Jones et al., 1970;Criales et al., 2006). Recently, Criales et al. (2015) found support for both suggested migration routes with a biophysical model utilizing Lagrangian modeling to display larval-to-postlarval behaviors, receiving output from a physical oceanographic model providing the drivers. The modeling system supported an investigation into the influence of spawning location, larval traits, and oceanographic features (tides, winds, and currents) on larval transport. Virtual larvae were released near the water column's surface from the Dry Tortugas and the Marquesas areas, mimicking a combination of DVM and STST behavior, and allowed to be transported for 28-30 days according to current speeds and directions and larval position in the water column (i.e., bottom vs. middle to surface). Finally, a benthic habitat module reflected larval aggregations on the offshore spawning grounds and suitable settlement habitats near the coast. The biophysical and physical oceanographic model developed by Criales et al. (2015) indicated that recruitment success was largely determined by season and spawning ground: generally, larvae FIGURE 1 | Individuals of Farfantepenaeus duorarum were collected from seven regions around the Florida Peninsula and Florida Keys, including two spawning aggregations, DryTortugas and Marquesas (yellow and underlined). The major (thick lines) and minor (thin line) migratory routes described by Criales et al. (2015) are depicted between the spawning aggregations and the major nursery grounds in Florida Bay (Photograph of F. duorarum from wikimedia commons. Base map from Google Earth). simulated from the Marquesas were several times more likely to reach nursery habitat than those simulated from the Dry Tortugas, and summer simulations consistently resulted in higher larval settlement compared to winter simulations. Simulated larvae were most likely to settle in nursery habitat when they were released from the Marquesas in the summer, migrating eastnortheast across the southwest Florida shelf. When simulated larvae originated from the Dry Tortugas, they were likely to become entrained in the Florida Current, exiting the Gulf of Mexico entirely and entering the greater Atlantic. The few simulated larvae released from the Dry Tortugas that successfully reached Florida Bay did so through both hypothesized routes, while those simulated larvae successfully recruited to the Florida Bay recruitment area from the Marquesas never migrated through the Florida Keys. These results provide expectations of population dynamics that can be tested with molecular methods.
The model of larval transport and migration developed by Criales et al. (2015) leads to testable, if relatively qualitative, hypotheses. Under the null hypothesis, all pink shrimp around the Florida Peninsula represent a single, genetically homogeneous population, originating from spawning aggregations offshore of the Dry Tortugas and the Marquesas, traveling either migratory route (Figure 1), and reaching adulthood on nursery grounds around the Florida Peninsula. From a fishery management perspective, this would be the simplest conclusion: lacking differentiated intraspecific diversity, all fisheries targeting the species can be managed as one. The alternative hypothesis, however, posits that the two spawning aggregations and different migratory routes to the nursery grounds support at least two genetically differentiated populations. If the alternative hypothesis holds, the Dry Tortugas and the Marquesas represent separate spawning aggregations to some extent, maintaining at least two distinct populations (these may be characterized by spawning aggregation, i.e., Dry Tortugas vs. Marquesas, or migratory route, i.e., the moretraveled east-northeast "major" route across the southwest Florida shelf vs. the less-traveled south-southeast "minor" route through the Florida Keys), and more complex management strategies would be needed to protect both populations during the stock-recruitment phase.
A better understanding of these two routes, major and minor, is of primary concern to researchers focused on sustainable fishing of pink shrimp (Browder et al., 1999(Browder et al., , 2002Ehrhardt and Legault, 1999;Criales et al., 2000Criales et al., , 2003Criales et al., , 2006Criales et al., , 2007Criales et al., , 2010Criales et al., , 2015Ehrhardt et al., 2001;Ogburn et al., 2013). The major route, which traverses the southwest Florida shelf, crosses through a regional fishery operating year-round near the Dry Tortugas and Key West (Klima et al., 1987;Upton et al., 1992;Hart et al., 2012), catching both fully mature and young adult shrimp (Ehrhardt and Legault, 1999;Browder et al., 2002). The co-localization of these large, highly productive pink shrimp fisheries with spawning grounds and out-migrating larvae makes an understanding of population dynamics in the region especially important to longterm species sustainability. Here, we utilize a next-generation sequencing method, double digest Restriction-site Associated sequencing (ddRADseq) to investigate the fine-scale population structure of F. duorarum in the eastern Gulf of Mexico. Our overall objective is to characterize diversity and connectivity in terms of the larval migration and transport within the area for the purpose of informing and improving fishery management. To accomplish this, we: (1) validate the biophysical oceanographic modeling results of Criales et al. (2015) with an independent data type (ddRADseq data); (2) investigate any evidence of population differentiation within F. duorarum in the region, including whether postlarvae recruited to Biscayne Bay originate from the Dry Tortugas; and (3) contextualize the population genomics results in terms of fisheries management.

MATERIALS AND METHODS
Because the migratory routes between pink shrimp spawning aggregations and nursery habitat span a relatively small geographic range, our sampling effort targeted proximal locations around the Florida Peninsula. Over 100 postlarval, juvenile, and adult specimens of Farfantepenaeus were collected from several sites representing seven regions around the Florida Peninsula between 2011 and 2015 (Figure 1): New Smyrna ("North_of_BiscayneBay" or "NBB"), Hobie Beach, Bear Cut, and South Virginia Key ("BiscayneBay" or "BB"), NOAA sampling stations 2.1-2.3 and 7.1-7.3 ("SouthBiscayneBay" or "SBB"), Bradley Key ("Everglades" or "EVG"), Pumpkin Bay, Estero Bay, Fakahatchee Bay, and Pine Island Sound ("North_of_Everglades" or "NEVG"), Fort Jefferson to Key West, which sampled across the Marquesas spawning ground ("Marquesas" or "MQ"), and the Dry Tortugas ("DryTortugas" or "DT"). The majority of samples collected from nursery habitats around the Florida Peninsula were acquired by Jackson as part of a collaboration between the Ecosystems Investigations Unit of the Southeast Fisheries Science Center (SEFSC) in Miami. South Biscayne Bay samples were collected as part of a nearshore southwestern Biscayne Bay monitoring project. Because some of the samples, primarily those representing spawning aggregations, were obtained from shrimping vessels, exact collection coordinates were not obtained. Sampled specimens were frozen after collection and shipped to the Ecosystems Investigations Lab at SEFSC for taxonomic identification, specifically focused on the diagnostic characters associated with reproductive morphology (gonopore, thelyca, and petasmata; see Pérez-Farfante, 1969, 1970, 1988Pérez-Farfante and Kensley, 1997). After identification to species, 105 frozen individuals identified as F. duorarum or likely to be F. duorarum (labeled F. sp.) were transferred to the CRUSTOMICS Lab in North Miami, Florida, where each was given a unique voucher ID in the Florida International University Crustacean Collection (FICC). The ID and all metadata associated with collection were entered into the FICC database. Samples were thawed and muscle tissue was plucked from each specimen by lifting the integument of the second abdominal segment and removing a few milligrams of tissue, using care to avoid puncturing the digestive tract. Tissue was stored at -20 • C in 70% EtOH. The intact whole-specimens were preserved in 70% EtOH and stored in the FICC. All specimens included in the study presented here, including all relevant metadata, are presented in Supplementary Table 1.

DNA Extraction and Next-Generation Sequencing Library Preparation
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.

Data Assembly and Quality Filtering
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 TABLE 1 | Details of data assembly in STACKS v1.45 are provided below, including flags and settings used at every step of the pipeline.

Population Genomics Analyses
Pairwise measures between regions, including Nei's unbiased genetic distances, which describe allelic differences assuming genetic drift and mutation are in equilibrium (Nei, 1972(Nei, , 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).

Validating the Existing Biophysical Oceanographic Model
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(Browder et al., , 2002Ehrhardt and Legault, 1999;Criales et al., 2006Criales et al., , 2007Criales et al., , 2015, describing a densitydependent 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-perindividual (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.

RESULTS
The preparation of next-generation sequencing libraries occurred for 68 individuals collected from 19 sites representing seven regions. Over 117 million SR150 raw reads were returned from the Illumina HiSeq4000. Demultiplexed data were submitted to the NCBI SRA database under BioProject PRNJA554161 and are also publicly available through the Gulf of Mexico Research Initiative's Information and Data Cooperative (doi: 10.7266/n7hhnq-kh83; Timm, 2019). Nine parameterizations of STACKS yielded nine data assemblies (batches, see Table 1) with 11,971-20,820 single nucleotide polymorphisms (SNPs). Additional quality filtering was executed in vcftools: when minimum read depth was set to 10×, 4,025-13,267 SNPs remained; applying a missing data filter (<90% individual missingness and <50% missing SNP data allowed) resulted in 740-800 high-confidence SNPs. BayeScan identified no loci under selection. See Table 2 for a detailed report of this information.
The sample sizes across regions included in the research presented here could be considered low compared to traditional population genetics studies of microsatellites or multilocus datasets. However, reduced representation library (RRL) approaches, such as ddRADseq, generate vastly more data, sampled from across the genome of each individual, and this increase in genomic data for each individual empowers the detection of fine-scale population structure with substantially fewer samples (Willing et al., 2012;Jeffries et al., 2016;Nazareno et al., 2017). North of Biscayne Bay Marquesas 10 9 10 10 10 9 10 10 10 Dry Tortugas 10 10 10 10 10 10 10 10 10 Information about the numbers of reads and SNPs passing each step of data assembly and filtering are reported in the upper section of the table, including the number of raw reads, the number of SNPs assembled within STACKS, and the number of SNPs that passed quality filtering (including minimum read depth of 10×, site missingness of 50%, individual missingness of 90% allowed, and removal of sites under selection). The lower section of the table reports final sample sizes for each region in the datasets that were analyzed in this study.

Population Genomics Analyses
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 amongpopulation 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  Supplementary Figures 1, 2 (Nei's and F ST , respectively). Line color indicates the type of region pair: nursery region to nursery region (red), nursery region to spawning region (purple), or spawning region to spawning region (blue). Line width indicates genetic distance with greater genetic distance or higher F ST illustrated with a narrower connecting line. Please note that these lines solely represent pairwise values, not movement of individual shrimp between regions. Spawning regions are labeled in yellow (Base map from Google Earth).  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).

Validating the Existing Biophysical Oceanographic Model
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 ).

DISCUSSION
The study presented here used next-generation sequencing data to inform management strategies by characterizing the population dynamics of F. duorarum around the Florida Peninsula, with specific focus on the role of migration from spawning aggregations to nursery grounds. Much of this work was motivated by the biophysical oceanographic model of FIGURE 4 | The number of private alleles, normalized by each region's sample size (referred to as PA norm in-text), are presented for each region. Detailed results for each dataset are presented in Supplementary Table 4. Box color indicates the type of region: nursery region (white) or spawning region (blue). Outliers are red. In all cases, the significance test, comparison of interest (Samples), test statistic value, and p-value are presented. *Indicates a statistically significant test (where p < α when α = 0.05).
larval transport from spawning aggregations offshore of the Dry Tortugas and Marquesas to nursery grounds in Florida Bay (Criales et al., 2015). The model supported two migration routes from spawning regions to nursery grounds: the major route crosses the southwest Florida shelf in a fairly direct eastnortheast path (Munro et al., 1968;Criales et al., 2003); the minor route involves downstream transport along the Florida Current, bringing larvae east-northeast with the Current and then breaking with the Florida Current to move west-northwest toward Florida Bay through the passes in the Middle and Lower Florida Keys (Jones et al., 1970;Criales et al., 2006). These two routes have the potential to sustain population differentiation within the species, representing overlooked biodiversity. Independent analysis of next-generation sequencing data revealed some population differentiation associated with the Dry Tortugas. With some caveats, the work presented here provides strong support for the model of larval migration and recruitment developed by Criales et al. (2015).

Utilizing Population Genomics Data to Validate a Biophysical Oceanographic Model
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.

Evidence of Population Structure in the Study 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 < nurserynursery 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.

Relevance to Fisheries Management
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.

AUTHOR CONTRIBUTIONS
TJ collected samples for inclusion in this study. JB and HB-G provided financial support for the research. LT performed the lab work and data analysis, generating an early version of this manuscript, including all tables, figures, and Supplementary Materials. All authors contributed to the study design and final manuscript.