ASSESSING HORSESHOE CRAB (LIMULUS POLYPHEMUS) POPULATION STRUCTURE WITHIN SOUTHERN NEW ENGLAND

The Atlantic horseshoe crab, Limulus polyphemus, is a commercially, ecologically, and economically important species. Current management practices of this species may not be well informed enough to avoid jeopardizing the future health of these animals. This thesis argues that population differences, as defined by morphometrics, behavior, and genomics, may be visible over smaller geographic scales than current fisheries management observes. Specifically, this work focuses on three states in southern New England: Rhode Island, Massachusetts, and Connecticut. Ten sampling locations were observed over two years, 2020-2021, from which over 500 crabs were sampled. Width and weight data were collected to assess whether size differs by location using nonparametric approaches. Tagging data from the US Fish and Wildlife Service was analyzed to assess whether small, localized movement patterns or broad range geographic movement was more prevalent throughout the range. Tissue samples were processed to extract genetic information (single nucleotide polymorphisms) to inform upon adaptive and migratory traits across the range. Morphometric data identified that 36% of pairwise comparisons were significant. Tagging data showed 69% of recaptured crabs were caught in the same water body they were originally released. Genomics tools suggested that outlier loci, more so than neutral loci, were driving the population structuring observed. Cumulatively, these results suggested that population differences can be observed over a smaller scale than currently employed for fisheries management.


LIST OF TABLES
and 18% of male and female, respective, prosomal width comparisons exhibited different mean ranks. Of those comparisons, 69% and 63% occurred across state lines for males and females, respectively. The same test was applied to weight. A total of 44% (45% interstate) of pairwise comparisons of male weight data and a total of 16% (100% interstate) of pairwise comparisons for female weight differed in mean ranks.
Overall, tagging data analysis identified the 68% of crab recaptures occurred in the same water body the crab was initially tagged in; however, these results are confounded by the large range in the sample size of recaptures by location (min n=1, max n=1449). The results in this study agree with previous work in Massachusetts that identified significant differences in horseshoe crab size over a short geographic scale.
However, this previous study found higher percentages of differences between sampling sites using fewer crabs and a parametric approach (ANOVA). The tagging results do not agree with many previous New England studies which identified very few crabs exiting the water body they were tagged in. All of the studies I analyzed used acoustic telemetry and observed fewer than 50 crabs. While the results in this study are not entirely consistent with previous work in the same geographic area, they do provide a baseline for future work that needs to be completed to inform sustainable management, most of which can be done through volunteer networks and educating the community of citizen scientists.

INTRODUCTION
The Atlantic horseshoe crab, Limulus polyphemus, is the only extant member of the Family Limulidae found in the Atlantic Ocean. This species spans from northern Maine to the Gulf of Mexico and the Yucatan Peninsula with the epicenter of biomass, and research, residing in the Delaware Bay Shuster and Sekiguchi 2009;Luo et al. 2020). The history of horseshoe crab harvest starts with the animals being used as fertilizer in the mid-1800s  To sustainably manage any species, especially one described as vulnerable by the IUCN red list, special attention must be given to accurately defining stock units as each unit may require different tools (Begg et al. 1999, Smith et al 2016

Field Sampling
A total of eleven sampling locations were selected based on historical horseshoe crab activity and suitable spawning habitat to represent the Southern New England region for this study ( Figure 1). This research relied heavily on Rhode Island, Massachusetts, and Connecticut's respective state environmental agency staff for guidance on sampling site selection. While crabs have been observed opportunistically utilizing a variety of substrates for spawning, most activity occurs in protected sandy bays, estuaries, and ponds , Smith et al 2002, Landi et al. 2015.
Eight locations were selected in Rhode Island and one location was selected from both Massachusetts and Connecticut. Ultimately, due to a lack of crabs, a second Massachusetts location was selected to supplement the first. For the entirety of this study, crabs from both Massachusetts sampling locations were treated as a single population.
Sampling occurred by hand collection between May and August over a two- year period (2020 and 2021) focused around the new and full moons. This approach ensured that adult crabs were being sampled as spawning activity typically coincides with major moon phases (Brockmann 1990). Crabs were also sampled by the Rhode Island Department of Environmental Management (RI DEM) and The Nature Conservancy staff opportunistically as encountered during seine surveys taking place at the sampling locations. The goal was to survey at least 30 crabs, split evenly between males and female, in each location in each of the two years. Once captured, crabs were separated by sex, measured across the widest point of the prosoma ( Figure   2), weighed, inspected for injuries, and affixed with a uniquely numbered US Fish and Wildlife Service (USFWS) horseshoe crab monitoring survey disc tag. For future genetic research, a segment of walking leg was clipped, placed in 5mL of DNA/RNA shield, and stored at -80˚C (Zymo Research). When possible, water temperature, salinity, and dissolved oxygen measures were collected at each sampling event.

Tagging Data Filtering
Horseshoe crabs have a terminal molt phase, relatively long lifespan, and few predators in their adult stage, making this species the perfect candidate for long-term tagging studies (Loveland 2002).All historical tagging records that contained a crab tagged and/or recaptured in Rhode Island were obtained from the USFWS in November 2021. The tagging program and subsequent data set originated in 1999 and has since continued to grow by providing tags to agencies along the United States' east coast for a multitude of long-term monitoring and research projects. This program relies heavily on citizen scientists and members of the public to report recaptures.
Once received, the data were filtered to standardize release and recapture locations to a scale fitting this project ( Figure 3). It should be noted that the standardized location of "Little Narragansett Bay" includes sections of both Connecticut and Rhode Island waters. For the purpose of this study, these crabs were treated as native to Rhode Island. Additionally, a tagging location of "Narragansett Bay" was added to represent all previous Narragansett Bay tagging events and tagging that took place at Conimicut point as part of this research. After standardizing, the complete data set was filtered to remove recaptures that took place within the same year the tag was affixed to the crab and multiple sightings of the same crab in the same month, year, and location. These methods were used to reduce inflation introduced by multiple people reporting the same crab, reflect the interannual time scale we aimed to observe, and capture any departures a crab made from their original tagging location.
A column for linear distance travelled and a column for direction of travel were added to the data set to estimate average bearings and to quantify significant movement.
Consistent with other studies, this work defines significant movement as linear travel greater than 3.5km . All data filtering was performed in R Studio(R Core Team 2022) using tidyverse (Wickham et al. 2019).

Analyses
Morphometric data, pooled by year, were analyzed using two approaches to determine if size differed among locations. The Kruskal-Wallis test was used as a nonparametric alternative to a one-way analysis of variance (ANOVA) to determine if size variables, grouped by location, originate from the same distribution (Kruskal and Wallis 1952). Pairwise Wilcoxon comparisons were then completed to determine which sampling locations differed from each other (Wilcoxon 1945). All p-values obtained from Wilcoxon tests completed in this work were corrected using the Bonferroni method to decrease chances of false discovery (Bonferroni 1935). This analysis was completed for both prosomal width and weight. Given the sexual dimorphism present in crabs, all morphometric analyses were performed separating for each sex (Botton and Loveland, 1992). Analysis of variance (ANOVA) could not be used with this dataset given the data did not meet the assumption of common variance in all populations (Eisenhart, 1947). Summary statistics of significant interstate comparisons will be utilized to further investigate the opposing state trends used in the 2019 ASMFC stock assessment.
Summary statistics were conducted to determine what proportion of crabs travelled to a location other than the water body they were tagged in. Tagging and recapture data also did not satisfy the ANOVA assumption of equal variance so again the Kruskal-Wallis and Wilcoxon tests were used to determine if distance travelled differs by tagging location. Recapture rates were low in some areas, so this analysis was not separated by sex. Additionally, this part of the analysis was only completed for crabs tagged in Rhode Island as this work is focused on characterizing Rhode Island crab behavior.

Morphometrics and Environmental Attributes
In 2020 and 2021, respectively, a total of 293 and 280 crabs were sampled and tagged (Table 1). Connecticut was the only site not sampled in 2020 (due to . The average male carapace width is 198.0mm and the average male weight is 1.1kg across all samples (Figure 4, Figure 5). The average female carapace width is 261.2mm and the average female weight is 2.8kg across all samples ( Figure 6, Figure 7). Kruskal-Wallis tests rejected the null hypothesis that width is consistent across sampling locations for both males and females (p <0.0001). Pairwise Wilcoxon tests showed significant differences in prosomal width between sites in 29% of male comparisons and 18% of female comparisons. Of these significant differences, 69% and 63% occurred across state lines for males and females, respectively ( The average water temperature during the 2020 sampling season was 17.7˚C (n=22) and 20.1˚C during 2021 (n=28). In 2020, the mean salinity (n=20) was 28.0ppt and the mean level of dissolved oxygen (n=20) was 10.3mg/L. In 2021, the mean salinity (n=15) was 26.9ppt and the mean level of dissolved oxygen (n= 16) was 8.7mg/L. Measures of salinity and dissolved oxygen were not obtained for Block

Island (2020 and 2021) and Massachusetts (2021). Mean environmental variable
values by year and location can be found in Table 4.

Tagging
Since 1999, a total of 6,106 crabs have been tagged in Rhode Island (41% female, 59% male, < 0.01% unknown). There have been 1,489 unique crabs recaptured (24% of releases) with 1,303 (21% of releases) occurring in Rhode Island. These metrics do not include crabs that have been recaptured more than once. The recapture rate for this study was higher than comparable studies by at least 10% (James-Pirri et al. 2005; James-Pirri 2010). The average time between release and recapture is 2.7 years with two crabs exceeding 20 years between release and recapture.
A total of 531 crabs (31%) were recaptured in a different location from where they were tagged (Table 5). The location with the highest proportion of crabs recaptured outside of their tagging location is Narrow River (100%); the location with the lowest proportion is Ninigret Pond (0%). The direction travelled by most crabs who left their tagging site was northeast (44%). The tagging location with crabs exhibiting the greatest mean distance of travel was Narragansett Bay (41km); the shortest mean travel distance greater than zero was 7.01km (Quonnie Pond). All mean travel distances were greater than the value that describes significant movement (3.5km).
There are four significant pairwise Wilcoxon results when examining differences in distance travelled by location (Block Island v Narragansett Bay, Block Island v Narrow River, Block Island v Winnapaug Pond, and Little Narragansett Bay v Narragansett Bay) ( Table 6). The number of recaptures by location varied widely (n=1-1449) ( Table 5). Another important study to examine in comparison to this study is Shuster's (1982) study, which places the largest horseshoe crabs in South Carolina and the smallest crabs at the northern and southern extremes of the species' range For males, this finding was consistent with the current study with the smallest median crab size (181.5mm) being observed in Massachusetts (most northern location) and the largest median crab size (~210mm) being observed in Connecticut and Ninigret Pond (most southern locations) (Figure 4). For females, however, this study observed the inverse.

DISCUSSION
The smallest median female crab size (245.mm) was observed in Quonnie Pond and the largest mean female crab size (277.5mm) was observed in Massachusetts ( Figure   6). Possible causes for this difference could be the higher statistical power provided by ANOVA analysis in the 1982 study or changes in size distribution in the past 40 years.
Tagging results, specifically the proportion of crabs recaptured in the same water body they were tagged (68%), support hypothesis 2: crabs are statistically more likely to be recaptured at their original tagging location than any other sampling location. .. Given the variance in estimates and low recapture rates, extrapolating the findings beyond the samples collected may not be warranted until supporting data can add more certainty to these findings. In contrast to the New England analyses performed by , James-Pirri (2010), and , all of which found few to zero instances of crabs exiting the water body they were tagged in, this work found instances of departure in 87% of tagging locations (

MANUSCRIPT II
Discerning horseshoe crab (Limulus polyphemus) population structure in southern New England using next generation sequencing.
Natalie Ameral 1,2* and Jonathan Puritz 1 To further distinguish the subtle signals present in these data it is necessary to expand the environmental variables present and increase the number of individuals per location.

INTRODUCTION
The Atlantic horseshoe crab (Limulus polyphemus), an invertebrate species with origins possibly as early as the Paleozoic era, finds itself approaching heightened concern in a conservation context. From fertilizer to bait to endotoxin detector, this previously widely unregulated species has served a variety of commercial purposes and now faces considerable population decline in recent years Berkson & Shuster, 1999;ASMFC 2019). Beyond commercial significance, the horseshoe crab plays a key ecological role as a food source for migratory shore birds, many of which are at risk of being endangered (Tsipoura & Burger, 1999). Along the east coast of the United States, the migratory bird species of greatest concern are red knot (Calidris canutus), ruddy turnstones (Arenaria interpres), and sanderlings (Calidris alba).
To assuage these pressures, the Atlantic States Marine Fisheries Commission (ASMFC) adopted the first horseshoe crab fisheries management plan (FMP) in 1998.
As with many FMP's, the goal for the horseshoe crab plan is to protect and conserve the species for current and future generations of humans and wildlife alike. A common theme throughout the history of this FMP and related stock assessments is the need for well-defined stock units (ASMFC, 1998; ASMFC, 2019).
The most recent coastwide genetic studies (King et al., 2005;King et al., 2015) identified 13 microsatellite DNA markers to analyze for population structuring using 1,684 crabs from Maine to Mexico. To summarize, this work found high regional FST values (82% pairwise FST significantly > 0), considerable allelic diversity and differentiation, and low gene flow and concluded that eight units would be appropriate for management (King et al., 2015). An extensive review of horseshoe crab conservation status summarized these eight regions down to six (Smith et al., 2017).

Additionally, Johnson's (2016) work analyzing 22 microsatellites using crabs in the
Cape Cod, Massachusetts area suggest embayment level structuring via PCA and GST.
Currently, ASMFC uses four regions for fishery management based on similarities in trends (from fishery independent and dependent catch data) and genetics information from King et al. in 2005 (Table 1). horseshoe crab population genetic structure occurs over finer scales (King et al., 2015;Johnson, 2016).
assessed by Benestan (2019), over 35% of recent genomic studies aimed at evaluating management units of exploited marine species uncovered a "mismatch" between population structure and current stock units. Another layer of complexity unique to the horseshoe crab fishery is that crabs are returned to the sea live after being bled for biomedical use; this means that associated mortality with this practice is significantly less than 100%. As it is difficult to enforce crabs be returned to their location of harvest, this species may be at risk of a fate similar to southern European Atlantic Salmon (Ayllon, 2006). This region used stocking programs to revive dwindling salmon populations which ultimately led to decreased genetic diversity of wild fish.
Returning crabs to a location other than the harvest site and restocking discrete water bodies with non-native fish are both examples of human admixture.
To further investigate population structuring that could inform upon the best management tools to implement in the New England region, the work presented here focuses on genetic analysis of southern Massachusetts and Rhode Island. Here, we used quaddRAD sequencing to identify single nucleotide polymorphisms (SNPs) present in crabs from nine sampling locations. SNPs were then separated into putatively adaptive and neutral datasets to analyze drivers of adaptation, structure, and estimated migration rates. This work hypothesized that population-genetic breaks will be discernable by sampling location given the population structuring identified by King et al. (2015) and Johnson (2016) and the higher power afforded by SNPs versus microsatellites. Ultimately, the goal of this research is to provide genetically informed guidance on whether current stock units in the Northeast are appropriately defined or if the species could be more appropriately managed under alternate designations.

Sample collection and library preparation
Adult horseshoe crabs were collected in 2020 from ten sampling locations (Manuscript 1 Figure 1 After quantification and qualification, we followed a modification of the double-digest restriction-site associated DNA protocol (Peterson et al., 2012) called quaddRAD (Franchini et al., 2017), which allows for higher-level multiplexing and detection of PCR duplicates using unique identifying barcodes. First, restriction enzymes were tested to find the optimal combination. Enzymes PstI, EcoRI, MspI, and SphI were tested for efficiency and selected based on estimated number of fragments produced (20,0000 -30,0000) and fragment size using Agilent TapeStation platform; PstI and EcoRI were chosen after analysis.

Bioinformatics and genotyping
Raw sequencing reads were deduplicated using the module clone_filter and demultiplexed using the process_radtags module of Stacks (Catchen et al., 2011(Catchen et al., , 2013) before being processed with the dDocent bioinformatic pipeline (Puritz et al., 2014). Briefly, dDocent removed low quality bases and adapters, mapped reads to the Resulting SNPs were first filtered using VCFtools to remove individuals that did not sequence well and include only variants genotyped in 50% of samples with a minimum quality score above 30 and a minor allele count of at least 3 (dDocent.com). Remaining SNPs were then filtered to keep SNPs with minor allele frequency of 0.05, filter populations by specific call rates, apply an allele balance between 0.25 and 0.75, remove sites with reads from forward and reverse strands, apply a mapping quality ratio between reference and alternate alleles of >0.9 and <1.05, apply a filter for properly paired status for reads supporting reference or alternate alleles, remove loci that did not have quality scores 2 times depth (Heng, 2011), filter loci above a mean depth of 280, and apply Hardy-Weinberg equilibrium.
In addition to accounting for all SNP elimination, these steps were included to reduce false positive calls (O'Leary et al., 2018).

Outlier Detection
The resulting genotypes were scanned for loci that have statistically higher than background FST values (outlier loci), looking for small portions of the genome that have larger differences in allele frequencies than expected by random (also known as genome scans based on population differentiation) (Hoban et al. 2016). The three methods used for outlier detection were pcadapt v.4.3.3 (Privé et al., 2020), OutFLANK v.0.2. (Whitlock & Lotterhos, 2015) and BayeScan v.2.1 (Foll & Gaggiotti, 2008). Pcadapt performs principal component analysis (PCA) on a centered and scaled genotype matrix then computes correlation values (p) between SNP's and PC's (Luu et al. ,2017;Privé et al., 2020). To select outliers, we applied a standard alpha level of 0.05 to the computed q-values (false discovery rate) of the p-values.
OutFLANK uses a VCF file of SNPs with minor allele frequency (MAF) >0.05 and a proxy for neutral FST, developed from the real dataset, to detect FST outliers with a threshold q-value of 0.1 (Whitlock & Lotterhos, 2015). OutFLANK has been shown to have a lower false positive rate than other outlier selection methods across a range of demographic histories and selection regimes (Whitlock & Lotterhos, 2015). BayeScan uses a Bayesian (posterior predictive distribution) approach to predict loci deviation (outliers) from Hardy Weinberg proportions using FST values allocated into population specific and locus specific segments (Foll & Gaggiotti, 2008). BayeScan was run using a thinning interval of 50, burn-in length 50,000, and 30 pilot runs. The threshold q-value used for detecting outliers with BayeScan was also 0.1. A dataset of putatively neutral loci was created by removing all identified outliers from the total SNP dataset.

Summary statistics and population structure
The outlier, neutral, and combined datasets were analyzed using the R package hierfstat to compute summary statistics and generate global and pairwise FST values (Goudet, 2004). Global FST examines allele frequency variation in an entire population (as a ratio of between-population variance (S) to total variance (T) while pairwise FST examines the subpopulation component of population structure (Wright, 1965;Reynolds et al., 1983;Weir & Cockerham, 1984). Heatmaps of pairwise FST values were generated to visualize population structure.

Redundancy analysis (RDA) is an extension of the multiple matrix regression
and is used to model multivariate response data (Legendre and Gallagher, 2001;Benestan et al., 2016). RDA, computed using adegenet (Jombart, 2008;Jombart & Ahmed, 2011) in R, was used to test for significant relationships between allele frequencies, environmental variables, and geographic distance in the outlier and neutral datasets. If water temperature, salinity, and dissolved oxygen measures were not collected in the field, the data were supplemented with temporally comparable data from RI DEM. The geographic distance was tested as distance-based Moran's eigenvector maps (dbMEMs) calculated using the R package, adespatial (Dray et al., 2022). Computing dbMEMs spatially decomposes and summarizes a study area into independent vectors that can be more easily used to describe spatial structure (Borcard & Legendre, 2002). Multicollinearity checks were run prior to identifying significant variables for use in RDA.
To discern genetic groupings while minimizing between-individual variation, we used discriminant analysis of principal components (DAPC) on the neutral and outlier datasets (Jombart et al. 2010). A thinned neutral loci dataset (5,413 SNPs) was also used to examine patterns of genetic diversity and infer small-scale patterns of genetic connectivity using estimated effect migration surfaces (EEMS) (Petkova, Novembre, & Stephens 2016). This program allows the user to calculate and visualize spatial variation in gene flow and genomic diversity by identifying geographic regions where genetic similarity decays more rapidly than expected. With this approach, the expectation is always isolation by distance and the program seeks to identify where this expectation is not supported. EEMS was run using three deme sizes (1000, 1200, and 1500) with a burn-in of 100,000 and MCMC length 500,000 iterations. The Reemsplots R package was used to average the results of the three runs and visualize effective migration rate (m) and effective diversity (q).

Bioinformatics and Outlier Detection
A total of over three million raw sequence reads resulted from Illumina sequencing.

Summary statistics and population structure
Expected heterozygosity, He, had a small range from a minimum of 0.3125 (outlier dataset) to a maximum of 0.3269 (neutral dataset). Observed heterozygosity, HO, followed a similar pattern with a minimum of 0.2977 (outlier) and a maximum of 0.3157 (neutral) ( Table 2). Global FST for the outlier dataset was 0.0072, an order of magnitude greater than both the neutral (0.0007) and combined (0.0008) data sets.
Mean pairwise FST followed the same pattern with less than a 0.0001 difference between global FST and mean pairwise FST for corresponding datasets (Table 2) Multiple runs of DAPC for the outlier and neutral datasets did not identify any significant groupings. With forced parameters, non-overlapping groupings could be discerned; however, as Jombart and Collins (2015) note, it is more important that clusters usefully describe the data, not that they merely exist.
All three estimated effective migration surface (EEMS) runs using 1000, 1200, and 1500 demes exhibited similar visual results meaning that averaging the three outputs would be appropriate. The averaged model identified two regions, Block Island Sound to the southern Rhode Island coast and Buzzards Bay, that displayed estimated migration (log(m)) > 0 and one region, Narragansett Bay, that displayed estimated migration (log(m)) < 0. The broad area centrally located on Figure 7 is attributed to a lack of data in this region. The probability plot reflects high confidence (95%) in a condensed area of Block Island Sound and the Conimicut Point sampling location ( Figure 8). The model did not identify any locations where mean effective diversity (q) differs from expected diversity (0.0) (Figure 9). shows there may be a barrier to gene flow.

Summary statistics and FST values
Although minute, lower observed than expected heterozygosity, as seen in both the neutral and outlier datasets, signifies that sampling localities are deviating from Hardy-Weinberg equilibrium, and this could be caused by recent barriers to migration.
For marine invertebrates, common barriers to gene flow include loss of habitat, patchy spawning habitat, ocean currents, and commercial harvest (Walsh et al., 2005;Cheng et al., 2020;Vu et al., 2021). Guided by other studies of marine invertebrates, neither the global nor pairwise FST values for the neutral data set would suggest high levels of population structure (Uthicke & Benzie, 2003;Lal et al., 2016). As this is the first genomic survey of horseshoe crab populations, it is still important interpret results relative to each other, especially in the context of fisheries management. The absence of significant differences in confidence intervals in the neutral data set suggests that there were significant historical levels of gene flow, that there is contemporary migration between populations, and/or that the populations are evolutionarily related (Hart & Marko, 2010). This result opposes typical behavior of horseshoe crabs which includes strong site fidelity and year-round occupation of small water bodies James-Pirri, 2010) and reports of low gene flow using microsatellites (King et al., 2105

Redundancy analysis and estimated effective migration surfaces
FST values informed upon whether each data set exhibited population structuring while RDA aimed to elucidate the drivers of both neutral and outlier population structure.
The inability of RDA to find an environmental or geographical correlate to describe the neutral data set is not surprising given the low FST results; however, the low FST does not immediately preclude the presence of a significant variable thus the test is still important to perform. The significant relationship between db-MEM3 and the outlier loci suggests that there is some localized adaptation present, but the environmental variable driving the selection is not contained in our dataset. Possible environmental variables to explore, given their status as drivers of selection in other marine invertebrates, are air temperature and pollution (Denis-Roy et al., 2020; Nielson et al., 2020). Another possibility is that the dataset could contain the environmental variable driving selection, but the number of observations for each variable may be too low or the accuracy too low to represent the trends in each location. Several previous studies that have found sea surface temperature and dissolved oxygen to be significant drivers of marine invertebrate population structure (Benestan et al., 2016;Bernatchez et al. 2018;Nielson et al. 2020) support attempting to compile a more robust dataset including the variables tested in this study.
Estimated effective migration surfaces (EEMS) are used to approximate possible historic migration patterns that would result in the diversity observed today with a baseline expectation that isolation by distance is occurring. The areas highlighted in blue on the EEMS map, which indicate greater than expected migration (or departure from isolation by distance), are consistent with FST values between sites ( Figure 7). The orange area of the map, located in upper Narragansett Bay and indicative of lower-than-expected migration, offers insight into a possible barrier to gene flow. The dynamic currents, mixing, or submarine topography of Narragansett Bay are all possible explanations of this relative isolation (Weisberg & Sturges, 1976;Spaulding & Swanson, 2008). A lack of confidence in the blue migration polygon in Massachusetts could be improved with the addition of more samples as this region had relatively low representatives (N=17).

Implications for fisheries management
This work illustrates that the strongest signal of population structuring is present in the outlier dataset. While overall gene flow is present, as demonstrated by the neutral loci, there may be selection pressures acting on putatively adaptive loci. This discovery is important especially regarding the return of bled crabs to their location of harvest. As we have seen, there is evidence of migration of individuals between sites. Negative impacts to populations can arise when individuals are returned to a site that does not match their fitness level as provided by localized adaptation (Schwenk & Spaak, 1995;Molofsky et al., 2014). Transplanting many unfit individuals over many generations could lead to a gene pool with less adaptive alleles present (decreased evolutionary potential); ultimately resulting in local extinctions (Bouzat, 2010;Markert et al., 2010). In addition to admixture, overharvest for commercial use could likely precipitate localized and total extinction. These concerns are exacerbated by climate change, the only solutions to which, besides human interference, are dispersal, phenotypic plasticity, and adaptation (Kirk & Freeland, 2011).
This work hypothesized that population-genetic breaks will be discernable by embayment. While the results presented here do not entirely support the hypothesized fine-scale population structuring, there is ample evidence that structuring, via putatively adaptive loci, is far more partitioned than current ASMFC management units. Notably, there are significant differences observed in adaptive allele frequencies over a 50km range of discrete coastal ponds and the current range of the management unit is over 500km of coastline. Identifying the drivers of selection, as well as increasing sampling sites and individuals per site, is crucial to perform before proposing management action. Possible measures to employ include stricter enforcement of returning crabs to their harvest location (to ensure that human-aided admixture is avoided) and location-specific harvest quotas. Ultimately, conserving this species means continued medical advancement, job security for commercial fisheries workers, and supporting fragile ecosystems.