Rangewide Population Structure of the Clearnose Skate Rangewide Population Structure of the Clearnose Skate

Skates (family Rajidae) are benthic elasmobranchs that are highly vulnerable to incidental ﬁ shery bycatch, are discarded at sea, and are poorly accounted for in catch records. Many aspects of skate life history, such as population structure, are not well understood. Without this knowledge, indiscriminate removal may have deleterious effects on scienti ﬁ c, conservation, and management efforts. The Clearnose Skate Rostroraja eglanteria is seasonally migratory and widely distributed in the coastal waters of the eastern United States and in the northeastern Gulf of Mexico. This study used molecular techniques to assess the population structure of Clearnose Skate for use as a biological reference point for further research and management. Specimens were collected from 2014 to 2019 by ﬁ sheries-independent surveys. High-throughput genotyping-by-sequencing was used to identify single nucleotide polymorphisms, resulting in two data sets: one consisting of 8,914 loci (outlier and neutral) and the other comprised of 30 outlier loci. Results from all analyses and using both data sets indicated a high level of genetic differentiation between specimens from the Gulf of Mexico and specimens from the U.S East Coast. Using the outlier data set, a low but signi ﬁ cant level of genetic differentiation was also found among specimens from the U.S. East Coast, with a subtle break near the North Carolina and South Carolina border. Genetic differences along the U.S. East Coast were spatially autocorrelated, indicating a latitudinal genetic gradient. The level of observed genetic differentiation between the Gulf of Mexico and the U.S. East Coast is likely due to physical barriers such as Florida and the Gulf Stream current, while the subtle structure along the U.S. East Coast is likely attributable to isolation caused by dispersal limitations and local temperature preferences. The results from this investigation of Clearnose Skate population structure can be used to better monitor and manage this vulnerable elasmobranch.

. In instances when harvest is tracked, skates are often placed into nonspecific skate groups (e.g., "skate unidentified") instead of speciesspecific listings. This method is problematic as grouped abundance may remain stable over many years when, in fact, local depletion of one or more species may be masked by the increased abundance of another (Dulvy et al. 2000). This can have deleterious effects on scientific, conservation, and management efforts (Marandel et al. 2018).
Similarly, understanding population structure is vital for characterizing other aspects of life history that are important for appropriate management, including demographic connectivity, reproductive patterns, segregation, migration, local adaptation, and phenotypic variation. Understanding each of these specifics is key for assessing biodiversity, conservation and risk assessment (Bräutigam et al. 2015;Domingues et al. 2017), biomass and harvest limits (Palsbøll et al. 2006;Dulvy et al. 2008;Ying et al. 2011), and phylogeny (Palumbi 1994;Natoli et al. 2003). Falsely assuming a single, panmictic population and mismanaging fishing pressure could result in unequal harvest rates, loss of genetic diversity, or stock collapse (Dulvy et al. 2000;Hueter et al. 2004;Laikre et al. 2005;Ying et al. 2011;Spies and Punt 2015).
Clearnose Skate Rostroraja eglanteria are flat, benthic, elasmobranchs and members of the family Rajidae (McEachran 2002). Clearnose Skate are distributed in the coastal waters of the western North Atlantic Ocean from Massachusetts to Florida and in the northeastern Gulf of Mexico. This species' wide latitudinal distribution is made possible by their range of temperature, salinity, oxygen concentration, and depth tolerances (Fitz and Daiber 1963;Schwartz 1996Schwartz , 2000Packer et al. 2003;Hogan et al. 2013;Schwieterman et al. 2019). There is evidence to suggest that their growth follows Bergmann's rule (Mayr 1956;Lindsey 1966). The reported maximum total length of Clearnose Skate from Massachusetts, North Carolina, and South Carolina typically ranges from 700 to 745 mm, 570 to 687 mm, and 568 to 640 mm, respectively (Bigelow and Schroeder 1953;Schwartz 1996). The length at maturity is estimated to be between 560 and 650 mm from New England skates (Sosebee 2004) and between 535 and 568 mm from South Carolina and northern Florida skates. Aging Clearnose Skate by counting vertebral bands has been challenging, but age-length calculations estimate maximum age to be 6-8 years (Fitz and Daiber 1963;Schwartz 1996;Gelsleichter 1998) and age at maturity to be 5-6 years (Northeast Fisheries Science Center 2000).
Clearnose Skate exhibit seasonal inshore-offshore migrations and are not present in all areas of their range during all seasons. Clearnose Skate in the Gulf of Mexico overwinter inshore and migrate offshore in the spring (Luer and Gilbert 1985), whereas Clearnose Skate between North Carolina and Delaware exhibit the opposite behavior by overwintering offshore and migrating inshore and northward during the spring (Packer et al. 2003). Some individuals visit Cape Fear River, North Carolina, between February and November, with an occasional absence in June and July when temperatures are highest (Schwartz 2000). It has been noted that Clearnose Skate between South Carolina and eastern Florida remain in coastal waters year-round (Bigelow and Schroeder 1953).
The reproductive patterns of Clearnose Skate are not well understood. Like other rajids, Clearnose Skate are oviparous and deposit numerous egg cases (colloquially referred to as mermaid's purses), each containing a single embryo (Last et al. 2016). The timing of reproduction and egg deposition coincide with regional migration patterns and temperature cues (Carrier et al. 2010). Clearnose Skate in the Gulf of Mexico mate in the winter and deposit eggs in the spring (Luer and Gilbert 1985;Rasmussen et al. 1999). For Clearnose Skate between Delaware and New York, egg deposition also occurs in the spring, although it is unknown where or when mating takes place (Breder and Nichols 1937;Breder and Atz 1938;Fitz and Daiber 1963). To add further uncertainty, female Clearnose Skate have been documented to store sperm in the shell gland of the oviduct for at least 3 months (Luer and Gilbert 1985). Some skate species have been observed to repeatedly utilize specific habitats as nurseries for depositing eggs, termed "regional philopatry," which may include their own place of emergence, termed "natal philopatry" (Mayr 1963;Quattrini et al. 2009;Hoff 2010;Hunt et al. 2011;Amsler et al. 2015;Chapman et al. 2015;Flowers et al. 2021), but it is unknown whether Clearnose Skate exhibit philopatry of any kind.
Efforts to manage Clearnose Skate have been sparse. They are not currently managed by the Southeast Fishery Management Council or the Gulf States Fishery Management Council but are managed by the New England Fishery Management Council (NEFMC) as part of a single management unit that spans from Maine to Cape Hatteras, North Carolina. This management unit includes six other co-occurring species of skate as a species complex. Inaccurate species identification prohibits the collection of accurate discard or landings records and biological data. The most recent stock assessment indicated that Clearnose Skate are not overfished and overfishing is not occurring (Northeast Fisheries Science Center 2007). The primary source of mortality comes as bycatch in the United States bottom trawl and dredge fisheries (Northeast Fisheries Science Center 2007) to which they are most susceptible during seasonal migration and congregation. It is assumed by the NEFMC that Clearnose Skate consist of a panmictic population, though lack of information regarding population structure has been consistently noted as a research 2 need and impediment to accurate management (Bonfil 1994;Northeast Fisheries Science Center 2007;Northeast Fisheries Management Council 2018;Anderson et al. 2020). Understanding this species' population structure would allow for subsequent delineation of management units, which is essential for preserving genetic diversity and abundance when faced with discard mortality and unequal fishing pressure across management councils. Additionally, if Clearnose Skate could be managed within the NEFMC skate complex more precisely, targeted skate species would potentially benefit from an increased total allowable catch (Hogan et al. 2013), thus ensuring the economic viability of skate fisheries in the United States (Northeast Fisheries Management Council 2020). The observed regional differences in Clearnose Skate growth, migration, and reproduction suggest some level of divergence from panmixia. The purpose of this study was to use genotyping-by-sequencing and single nucleotide polymorphisms (SNPs) to investigate the population structure and genetic diversity of the Clearnose Skate, thereby testing the null hypothesis that this species consists of a single, panmictic population. Specimens collected from throughout the Clearnose Skate range were examined along both broadscale and small-scale geographic delimitations.

Specimen Collection, DNA Extraction, and Filtering
Specimen collection. with an otter trawl. The area covered by these surveys was hierarchically stratified using two geographic scales. For broadscale analysis, the range of specimen collection was delimited into three study regions bounded by biogeographic barriers ( Figure 1): (1) a "northern range" study region (NOR) consisted of the NEAMAP and VIMS sea scallop surveys, which begin at Cape Hatteras and extend northward, (2) a "southern range" study region (SOU) consisted of South Carolina Department of Natural  Table 1 available in the online version of this article). Fin, muscle, or heart tissues were collected from skate specimens either at sea by survey crews or from frozen survey specimens that were shipped to VIMS, and the tissue was stored in 95% solution of ethanol at −20°C until processing.
Extraction of DNA.-All tissues were first soaked in Longmire lysis buffer between 2 and 12 h at room temperature (Longmire et al. 1997) followed by DNA extraction using either Qiagen DNeasy Blood and Tissue kit (Qiagen Sciences, Germantown, Maryland) or Puramag carboxylated magnetic beads (MCLAB, South San Francisco, California) using standard protocols with minor modifications. To achieve increased DNA concentration for a subset of specimens with poor quality tissue, two, three, or four extractions per specimen were performed, eluted in 50 μL of elution buffer (Qiagen), and subsequently combined and concentrated using the Zymo Research DNA Clean & Concentrator kit (Zymo Research, Irvine California). The DNA quality was assessed on 1% agarose gels that included GelRed Nucleic Acid Stain (Biotum, Fremont, California) and a 1Kb Plus ladder (Invitrogen, Carlsbad, California). The DNA concentration and purity were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts).
Next generation sequencing.-Extractions yielding highmolecular-weight DNA and a concentration of at least 50 ng/µL were sent to Diversity Arrays Technology (DArT Pty, Canberra, Australia) for high-throughput genotypingby-sequencing using DArTseq (Sansaloni et al. 2011). Briefly, restriction endonucleases were used to reduce the complexity of the genome followed by hybridization to microarrays and sequencing on an Illumina HiSeq 2500 (Illumina, San Diego, California). Then, SNPs were identified using a proprietary data pipeline developed by DArT PL (DArTsoft14), where fragments were aligned to the Little Skate Leucoraja erinacea reference genome (Wang et al. 2012; http://www.diversityarrays.com).
Filtering of SNPs.-The resulting matrix of identified SNPs and corresponding metrics for each specimen were downloaded from the Diversity Arrays Website and subsequently underwent further quality control filtering steps. We used the "dartR" package (version 1.3.5; Gruber et al. 2019) in R 3.6.1 (R Core Team 2019) to remove lowquality loci and individuals as follows: (1) monomorphic loci, (2) read depth using a lower bound of 5 and upper bound of 30, (3) reproducibility average of 1.00 across technical replicates, (4) a call rate by locus threshold of 0.98, (5) a minor allele frequency threshold of 0.01, (6) hamming distances with a threshold of 0.25, (7) the R package "radiator" (version 1.1.6; Gosselin et al. 2020) used to identify loci out of conformance to the expectations of Hardy-Weinberg equilibrium in at least two populations and with a P <0.05, (8) secondary SNPs, and (9) call rate by individual using a threshold of 0.95. This resulting data set is hereby referred to as the "full data set."

Genetic Variation and Population Structure
Principal component analysis, heterozygosity, and relatedness.-Principal component analysis (PCA) plots were created to visually examine genetic variation and clustering patterns among all sampled individuals in the R package "adegenet" (version 2.3.1; Jombart et al. 2008;Jombart and Ahmed 2011). Individuals that were plotted far from their respective clusters as outliers were examined for evidence of relatedness, improper species identifications, or cross contamination resulting from errors in laboratory procedures by estimating individual heterozygosity levels in "radiator." Outlier individuals were also examined for potential familial relationships using relatedness values (r). The r-values were calculated in the Coancestry software (Wang 2011) using the triadic likelihood estimator (Wang 2007) and were then compared with values of four simulated relationship categories (full sibling, halfsibling, parent-offspring, and unrelated). Ten simulations were conducted using the R package "related" (Pew et al. 2015), each with an allele frequency matrix calculated using 100 randomly selected loci and 100 simulated relationship pairs. The resulting simulated genotypes were analyzed using the triadic likelihood estimator with an error rate of 0.001 and were used to create box plots to illustrate the range of expected relatedness values recovered for each family relationship category. Calculated rvalues for each specimen were compared with the simulated results and specimens with relatedness values consistent with parent-offspring, full sibling, or half-sibling. Specimens with values similar to or above these categories were removed to remove potential bias.
Values of F ST .-We estimated the relative levels of differentiation, using Wright's F-statistics (Weir and Cockerham 1984), among samples at both broad (study region) and small (state) geographic scales in the program GenoDive version 3.03 (Meirmans 2020). Significance was assessed using an initial α = 0.05, and false discovery rate (FDR) corrections were applied to provide more conservative estimates of significance (Benjamini and Hochberg 1995).
Analysis of molecular variance.-Hierarchical groupings of the genetic data were explored by conducting an analysis of molecular variance (AMOVA; Excoffier and Lischer 2010), which tests grouping scenario hypotheses based on results from the PCA and F ST results. Tests were run in Genodive using the infinite allele model, and significance was assessed using 10,000 permutations of the data.
Identification of population clusters.-Identification of the number of population clusters (i.e., population groups) and an individual's cluster assignment was examined using two methods. First, discriminant analysis of principal components (DAPC; Jombart et al. 2010) was conducted using "adegenet." The number of retained principal components was determined using two approaches: a visual assessment of the graphs showing the variance explained by the PCA and by the function "xvalDapc" that uses cross-validation procedures to optimize the number of principal components retained. The optimal number of clusters, K, was determined using the Bayesian information criterion (BIC), and individuals were probabilistically assigned to each of the identified clusters. The second method used Bayesian clustering implemented in the program STRUCTURE (Pritchard et al. 2000). The likelihood of each K was estimated using the admixture model, sampling locations as prior information, and correlated allele frequencies. Five replicates of simulations of K = 1-5 were run using a burn-in of 250,000 Markov chain−Monte Carlo (MCMC) iterations followed by an additional 500,000 MCMC iterations. Results from each simulation were imported in the R package POPHELPER version 2.3.0 (Francis 2017), tabulated, and visualized using bar plots, and the likely K was determined using ΔK and maximum likelihood criteria based on the Evanno method (Evanno et al. 2005). The program STRUCTURE has a tendency to fail to detect subtle population differentiation in cases of hierarchical population structure or when the data is obscured by strong genetic differentiation (Evanno et al. 2005;Janes et al. 2017;Roycroft et al. 2019); therefore, a second data set without the GOM samples was used to determine K for the NOR and SOU samples as above. After the optimal K was selected for each data set (with and without the GOM samples), individual assignment to each K was determined by running STRUCTURE with only the selected optimal K, a 500,000 MCMC iteration burnin, and subsequent 2,000,000 MCMC iterations. All results were imported in POPHELPER, and individual population assignment was visualized as STRUCTURE bar plots.

Spatial Autocorrelation
Spatial patterns of genetic variation were assessed using spatial principal component analysis (sPCA; Jombart et al. 2008) in "adegenet." The sPCA method is similar to that of a PCA in that it identifies genetic variation, but it additionally calculates and incorporates Moran's index, I, a measure of spatial autocorrelation (Moran 1948(Moran , 1950). An individual's genotype was compared with that of its neighbors and examined for whether the genotype was more similar or more dissimilar than what would be expected in a random spatial distribution of the allele frequencies. Neighbors were determined using a Euclidean distance-based network with a minimum distance of zero and a maximum distance of six. Missing allele frequencies were replaced with the mean allele frequency, and the presence of global structure (i.e., spatially autocorrelated) and local structure (i.e., ecologically driven) was tested using 999 permutations of the data.

Genetic Diversity and N e within Identified Populations
Once putative genetic populations were identified, summary statistics including expected heterozygosity (H e ) and observed heterozygosity (H o ) were calculated in GenoDive. Additionally, an inbreeding coefficient, F IT , was calculated for each specimen in Coancestry. For each putative population, contemporary effective population size, N e (Wright 1931), was calculated using NeEstimator version 2.1 (Do et al. 2014) with the linkage disequilibrium method and for three allele frequency thresholds: 0.01, 0.02, and 0.05.

Outlier Loci
Outlier loci are markers that exhibit patterns of variation that are divergent from the rest of the genome and can be used to identify population structure related to local adaptation (Luikart et al. 2003;Russello et al. 2012). Outlier loci were identified using the dartR wrapper Out-FLANK (Whitlock and Lotterhos 2015) with a q-value threshold of 0.05, left and right trim fractions of 0.05, and minimum heterozygosity at 0.10. The identified outlier loci were saved as a separate data set and were used to evaluate population structure using PCA, F ST , DAPC, STRUC-TURE, AMOVA, and sPCA methods as described above.

Specimen Collection and Data Filtering
Clearnose Skate DNA was extracted from a total of 197 specimens across the three

Population Structure
Principal component analysis, heterozygosity, and relatedness.-An initial PCA plot was used to visually explore the genetic differentiation among Clearnose Skate samples using the full data set (Supplementary Figure 1 available in the online version of this article). Using broadscale groupings, the NOR and SOU specimens formed a single large cluster and the GOM specimens formed a second cluster. The two clusters were highly separated from each other along PC1, which suggested a high level of genetic differentiation. The PCA plots also showed three pairs of individuals, two from the GOM region and one from the NOR region, that were plotted as outliers. Individual heterozygosity was calculated and examined for evidence that these outliers were the result of cross contamination. The range of individual heterozygosity among all samples (0.111-0.253) was not consistent with evidence of cross contamination or laboratory procedure errors, and no observed individual heterozygosity was anomalously high. To discriminate whether these outliers represented duplicated, misidentified, or related samples, relatedness values were calculated and compared to simulated values. Based on the genetic differentiation evident in the PCA plots, r-values were calculated for GOM individuals separately from the NOR and SOU individuals (Supplementary Figure 2). The estimated means and interquartile ranges of each family relationship category for NOR and SOU specimens were as follows: siblings = 0.5603 (0.4791-0.6518), half siblings = 0.3023 (0.1949-0.4113), parent-offspring = 0.5688 (0.5000-0.6341), and unrelated = 0.0777 (0.0-0.1250). The estimated means and interquartile ranges for the GOM specimens were as follows: full sibling = 0.5677 (0.4735-0.6644), half siblings = 0.3131 (0.1741-0.4361), parent-offspring = 0.5888 (0.5017-0.6704), and unrelated = 0.01036 (0.0-0.1607). The estimated r-values of the two GOM-GOM sample pairs were 0.79 and 0.72, which is greater than the third quartile of simulated values for any relationship category but less than values that would be expected if samples were duplicates. The relatedness value of one NOR-NOR sample pair was estimated to be 0.47, consistent with either a full sibling or a parent-offspring pair. To avoid bias that may result from the sampling of family members or duplicated samples, one individual from each of those three pairs was deleted from the data set, thereby reducing the final number of individuals to 194 (NOR = 127, SOU = 45, GOM = 22; Figure 1). The PCA was recalculated using the final data set, and outlier specimens were no longer visible. For this, PC1 explained 4.73% of the variation in the data, and PC2 explained 0.74% of the variation in the data (Figure 2A). To further examine patterns among the NOR and SOU specimens, the GOM specimens were removed from the data set and the PCA was recalculated. For this subset, PC1 explained 0.85% and PC2 explained 0.84% of the variation of the data ( Figure 2B). Specimens from both study regions formed a single cluster in the PCA plot, and there was no suggestion of differentiation. To look for evidence of clustering of individuals by the state in which they were recovered, the PCA was recalculated and plotted again, but no clustering patterns were evident (data not shown).
Values of F ST .-Population differentiation was quantified using pairwise F ST values calculated among both broadscale study regions and among states. Among study regions, the pairwise F ST value after the FDR correction between the GOM and NOR samples was significant (F ST = 0.108, P ≤ 0.001) and values between the GOM and SOU samples was also significant (F ST = 0.107, P ≤ 0.001). The F ST value between the NOR and SOU regions was low but significant (F ST = 0.001, P ≤ 0.001). At the state level, F ST values were highest between the GOM specimens and specimens from all other states (0.109-0.131; Table 1), and all values were significant after the FDR correction. Pairwise comparisons among samples from U.S. East Coast states had low F ST values (0.000-0.003), and the magnitude of estimates did not change with increasing geographic distance. Comparisons that were statistically significant after the FDR correction were between SCGAFL and six other states: RI (0.003, P = 0.013), NY (0.002, P = 0.001), NJ (0.002, P < 0.001), DE (0.003, P = 0.007), VA (0.002, P < 0.001), and NCN (0.002, P = 0.002).
Analysis of molecular variance.-An AMOVA was used to identify the most optimal grouping of specimens by testing three hypotheses ( Table 2). The first hypothesis, AMOVA 1, was the null hypothesis and had no grouping scheme and maintained state-level delimitation. The second hypothesis, AMOVA 2, partitioned samples into three groups to test whether Cape Hatteras is a strong barrier to gene flow: RI-NCN, NCS-SCGAFL, and GOM. The third hypothesis, AMOVA 3, also placed states into three groups, but these groups corresponded to the results from the pairwise F ST values: RI-NCS, SCGAFL, and GOM. The results indicated that the best-supported scenario was AMOVA 3, which maximized variance among groups (F CT = 0.055, P < 0.001) and minimized variance within groups (F SC = 0.002, P < 0.001; F ST = 0.090, P < 0.001).
Cluster delimitation.-The clustering methods in the DAPC and STRUCTURE analyses were used to explore the 6 NELSON ET AL.  most likely number of population clusters and individual sample assignment to each cluster. Attempts to retain principal components for the DAPC using the "xvalDapc" cross-validation function were unsuccessful as the function estimated retaining a low number of principal components, which corresponded with low percentages of the variance explained. Therefore, only visual assessments of the graphs were used to select the number of retained principal components. For DAPC, the lowest BIC value indicated that the most likely number of clusters was K = 2 (BIC = 1,128.1; Supplementary Figure 3), and the plots suggested that one cluster was comprised of the NOR and SOU specimens and the second cluster was comprised of the GOM specimens ( Figure 3A). Cluster delimitation using STRUCTURE also indicated that two clusters were most optimal according to the largest ΔK and mean loglikelihood values (Supplementary Figure 3), and the pattern of assigning the NOR and SOU specimens into one cluster and the GOM specimens into the second cluster was consistent among runs ( Figure 3B). An attempt was made to look for structuring patterns among the NOR and SOU specimens only, but the results identified K = 1 as the most optimal, which is consistent with a single genetic population.

Spatial Autocorrelation
Evidence of spatially autocorrelated allele frequencies was tested using an sPCA. The global test for spatial autocorrelation using the full data set indicated that the probability of global structure was significant (score = 0.028, P ≤ 0.001), but the probability of local structure was not significant (score = 0.718, P = 0.082). The first principal component, PC1, explained the majority of the genetic variance and a high degree of spatial autocorrelation (variance = 149.390, I = 0.917), and PC2 explained a modest amount of genetic variance but a substantial amount of spatial autocorrelation (variance = 18.348, I = 0.422). These two principal components were retained for plotting. The plot of spatial principal component 1 (sPC1) revealed two distinct clusters: one comprised of the GOM   Figure 1 for abbreviations), and color represents the probability of assignment to cluster 1 or 2. In both instances, K = 2, shown here, is the most supported scenario. 8 specimens and the second of the NOR and SOU specimens ( Figure 4A). The plot of sPC2 revealed a weak genetic cline along the geographic range of the NOR and SOU specimens ( Figure 4B).

Genetic Diversity and N e within Identified Populations
Observed heterozygosity (H o ), expected heterozygosity (H e ), and inbreeding coefficient (F IT ) values were calculated for each state individually and for the three putative populations identified by previous analyses (hereafter, referred to as RI-NCS, SCGAFL, and GOM, according to the associated states). The genetic diversity of the GOM population (H e = 0.132) was lower than both the SCGAFL (H e = 0.159) and RI-NCS populations (H e = 0.161; Table 3). Within each population, H o was less than H e . The significance of this observed heterozygote deficiency within each population was tested using a paired ttest, and results were significant for all three populations (RI-NCS: t = −44.84, P < 0.001; SCGAFL: t = −23.146, P ≤ 0.001; GOM: t = −8.910, P < 0.001) and overall Effective population size, N e , was also estimated for each putative population. For the tested allele frequencies (0.01, 0.02, 0.05), the N e estimate and range for the RI-NCS population was 11,859.0 (5,501.6-infinite), 13,018.0 (5,371.4-infinite), and 14,863.6 (6,426.3-infinite), respectively, the N e estimate and range for the SCGAFL population was infinite (infinite-infinite), 18,381.4 (4,415.1infinite), and 10,459.3 (3,653.1-infinite), respectively, and the N e estimate and range for the GOM population was 8,167.6 (1,578.4-infinite), 8,167.6 (1,578.4-infinite), and 2,005.8 (819.6-infinite), respectively (Table 4). In addition, N e was also estimated for a metapopulation that combined both the RI-NCS and SCGAFL populations, which were 16,359.8 (7,485.9-infinite), 19,992.5 (7,735.4-infinite), and 21,157.0 (8,976.2-infinite) for each of the allele frequencies, respectively.

Outlier Data Set
Principal component analysis and F ST .-An outlier loci data set was identified and used to test for genetic differentiation that may be indicative of selection. Using the OutFlank command, 30 outlier SNPs were identified from the full data set and used in subsequent analyses. In the PCA plot of outlier loci, PC1 explained 60.70% of the variation in the data and PC2 explained 9.65% ( Figure  2C). Similar to the PCA plot using the full data set, the NOR and SOU specimens clustered together, while the GOM specimens formed a second cluster, with the two clusters separating along PC1. However, when the GOM specimens were removed and the PCA was recalculated, there appeared to be two clusters ( Figure 2D). One tightly clustered group consisted of mostly NOR and a few SOU specimens and a second, diffusely clustered group with roughly equal proportions of NOR and SOU specimens intermingled.
Similar to the data set using all loci, the F ST values between samples collected from the GOM region and those collected from all other states were the highest (F ST = 0.685-0.802; Table 1) and were significant after FDR correction (P < 0.001). Pairwise comparisons among the U.S. East Coast states had F ST values ranging from −0.039 (RI versus DE, P = 0.834) to 0.195 (NY versus SCGAFL, P < 0.001). Comparisons between specimens from the SCGAFL state and specimens from the following five states were significant after  Analysis of molecular variance.-Three AMOVA groupings identical to the analysis using the full data set were used to identify populations. The results indicated that when samples were grouped as RI-NCS, SCGAFL, and GOM populations, the variance among groups was maximized (F CT = 0.631, P < 0.001) and the variance within groups was minimized (F SC = 0.002, P < 0.001; F ST = 0.001, P < 0.001), consistent with the results from the full data set.
Cluster delimitation.-Using the DAPC method, the most likely number of clusters was either K = 2 (BIC = 5.41) or K = 3 (BIC = −73.3) (Supplementary Figure 4). Using K = 2, the DAPC assigned the GOM specimens to one cluster and the NOR and SOU specimens to the second cluster (data not shown). Using K = 3, the GOM specimens were assigned to one cluster and the NOR and SOU specimens were assigned to one of the remaining clusters, but there was no evident pattern to the assignment of individuals by study region or state. To examine cluster assignment for just the NOR and SOU specimens, the GOM specimens were removed and the DAPC was recalculated. The most likely number of clusters was possibly K = 2, 3, or 4 (BIC = −55.8, −72.6, −83.7, respectively; Supplementary Figure 4). The DAPC bar plots were made for each possible value of K ( Figure 5A), from which there were no clear patterns of assignment to clusters for any value of K.
For assignment tests using STRUCTURE, K = 2 was most optimal (Supplementary Figure 4), where the GOM specimens were assigned into one cluster and the NOR and SOU specimens were assigned into the second cluster (data not shown). Again, the GOM specimens were removed to more closely examine cluster assignment using just the NOR and SOU specimens. According to ΔK values, K = 2 was most optimal but maximum log-likelihood values suggested that K = 3 and K = 4 were possible scenarios and were therefore also examined (Supplementary Figure 4). Cluster assignment did not reveal any discernable pattern in the bar plots of K = 2-4 ( Figure 5B).
Spatial autocorrelation.-Analysis using the outlier data set indicated that the probability of global structure was significant (score = 0.585, P ≤ 0.001), but the probability of local structure was not significant (score = 0.014, P = 1.00). The sPC1 (variance = 15.78, I = 0.877) revealed strong separation between the GOM population and both RI-NCS and SCGAFL populations except for the two southernmost U.S. East Coast specimens ( Figure 4C). When only the RI-NCS and SCGAFL specimens were examined, global structure remained significant (score = 0.40, P = 0.001) and local structure was not significant (score = 0.014, P = 0.962). The sPCA plot revealed a genetic cline among the RI-NCS and SCGAFL specimens that was weakly spatially autocorrelated (variance = 3.473, I = 0.096; Figure 4D).

DISCUSSION
Few studies have examined the genetic structure of skates, and as of this writing, only one other study has utilized genotyping-by-sequencing methods to examine population structure in rajids (Wang et al. 2012). The present study was also the first to examine the rangewide population structure and genetic diversity of the Clearnose Skate, and we posit that they are comprised of three populations based on the genetic data. All analyses supported a significant level of genetic differentiation between samples FIGURE 5. Bar plots of individual Clearnose Skate cluster assignments according to (A) discriminant analysis of principal components and (B) STRUCTURE analyses using the 30 outlier loci data set. State and study region are shown below the bar plots (see Figure 1 for abbreviations), and color represents the probability of assignment to each of the cluster scenarios (K = 2, 3, and 4).

POPULATION STRUCTURE OF CLEARNOSE SKATE
collected from the Gulf of Mexico (GOM) and those collected in other locations. The data also suggest weak genetic separation between samples collected from off Rhode Island to North Carolina (RI-NCS) and samples collected from off South Carolina, Georgia, and Florida (SCGAFL), though this may be confounded by the distribution of specimens throughout the study area.

Evidence of a Distinct Gulf of Mexico Population
The mechanisms driving the observed separation of Gulf of Mexico and Atlantic Clearnose Skate may include both historical geologic patterns and present-day conditions. During periods of glacial minima in the Pleistocene (1.8 million to 12,000 years ago; Holmes 2001), the combined effects of ice retreat, ice melt, and isostatic changes produced eustatic sea level rise and shallow connections across peninsular Florida. However, at the time of the last glacial maximum about 20,000 years ago, eustatic sea level was estimated to be 120 m below current levels, exposing Florida and the gulf-side continental shelf, thereby separating the ocean basins (Briggs 1974;Avise 1992;Holmes 2001). Historical Clearnose Skate populations may have been subject to a series of range contraction and range expansion events caused by fluctuations in sea level. The identification of Florida as a biogeographic break for this species is concordant with results from studies of many taxa, including elasmobranchs, teleosts, marine mammals, and bivalves (Briggs 1974;Avise 1992;Natoli et al. 2003;Hemond and Wilbur 2011;Anderson et al. 2012;Portnoy et al. 2015).
We hypothesize that the present-day isolation of Clearnose Skate in the Gulf of Mexico may be maintained by the geography and hydrodynamics around Florida. First, the continental shelf between Cape Canaveral and the Florida Keys essentially vanishes, thereby reducing habitat availability. This area might be too steep or too deep to pass through or may not host suitable benthic prey items. Second, suitable habitat around peninsular Florida is further reduced by warm water temperatures that are outside the Clearnose Skate's preferred range (9-30°C). This thermal gradient is also hypothesized to prevent Atlantic angel sharks (family Squatinidae) from migrating from their U.S. East Coast range into their Gulf of Mexico range (Driggers et al. 2018). Lastly, the velocity of the Gulf Stream may be strong enough to prevent or deter individuals from passing through, particularly in a countercurrent direction from the U.S. East Coast into the Gulf of Mexico.

Two Populations along the U.S. East Coast
This study found evidence to suggest a low but significant level of genetic structure between two populations along the U.S. East Coast, RI-NCS and SCGAFL, delimited near the North Carolina and South Carolina border. Although the results were not unanimous, this differentiation was evident from the AMOVA and pairwise F ST results and from both the full and outlier data sets. The ability to accurately determine the number of populations, K, and individual sample assignments using DAPC and STRUCTURE with the outlier data set were likely hindered by the low levels of genetic differentiation (Miller et al. 2020). The results from the sPCA using the outlier data set showed increasing genetic differentiation with increasing distance, which suggests that Cape Hatteras is not a barrier to gene flow for the Clearnose Skate and that other factors are responsible for the observed distribution of genetic variation.
Estimated F ST values for marine fishes are typically lower than values for freshwater and anadromous fishes, often owing to the lack of physical barriers that would otherwise separate populations (Ward et al. 1994). Although the reported F ST values using the full data set were low, there was a statistically significant level of differentiation in pairwise comparisons between the SCGAFL samples and those from six other and more northern states, a pattern that was also seen using the outlier data set. Similarly low yet significant pairwise F ST values have been reported in other fish population studies, including Nassau Grouper Epinephelus striatus (F ST = 0.0023, P = 0.0039 and F ST = 0.002, P = 0.0140; Jackson et al. 2014), Leopard Sharks Triakis semifasciata (F ST = 0.003, P = 0.026; Barker et al. 2015), and South Pacific Albacore Thunnus alalunga using both neutral loci (F ST = 0.0027-0.0031, P < 0.0004) and outlier loci (F ST = 0.0157-0.0203, P < 0.0012) . We feel confident that the number of loci used in this study was sufficient to overcome the high signal-to-noise ratio often seen in marine organism F ST estimates but acknowledge that sample size and distribution may bias estimates (Waples 1998;Willing et al. 2012).
The pattern observed between the RI-NCS and SCGAFL populations could be described as having "crinkled connectivity" driven primarily by demographic life history differences (Ovenden 2013) rather than physical barriers, as seen in other fish species (Briggs 1974;Avise et al. 1987b;McCartney et al. 2013;Leidig et al. 2015). Despite continuous distribution along the U.S. East Coast, the dispersal capabilities of the Clearnose Skate may be responsible for the observed heterogeneity (Ovenden 2013). Like other elasmobranchs, skates produce relatively few young and disperse as juveniles or adults, which must then survive, assimilate into, and reproduce with a different population to maintain genetic connectivity. This is in stark contrast to many teleost fishes that reproduce by spawning thousands to millions of eggs, which may be subject to transport through currents or waves and lead to high genetic connectivity. Temperature preferences may also reduce connectivity. Throughout most of their range, Clearnose Skate undergo temperature-driven inshore-12 offshore migrations, except for skates from the SCGAFL area, which exhibit year-round inshore residency, which may act as a barrier to gene flow. Angel sharks also exhibit discontinuous distribution patterns along the U.S. East Coast, which has been attributed to temperature preferences (5.5°C to 26.7°C) and the bottom temperature profiles in that region. (Driggers et al. 2018). The Gulf Stream current carries warm water from Florida northward to Cape Hatteras, North Carolina, thereby making offshore bottom temperatures warmer and more preferable than the cold inshore bottom waters during winter months. Clearnose Skate likely maintain a similar relationship between thermal tolerances (9°C to 30°C) and seasonal water temperature fluctuations, but investigation of environmental variables associated with Clearnose Skate catch data and habitat modeling is needed to confirm this hypothesis. The complex combination of localized site fidelity, broad seasonal inshore-offshore migration, limited dispersal capabilities, and wide geographic range likely explain the low but significant level of genetic differentiation.

Genetic Diversity and N e within Identified Populations
Estimates of allelic diversity and N e suggest that the GOM population is more susceptible to the effects of genetic drift than the other two populations. The GOM population had the lowest gene diversity (H e ) values, the highest number of monomorphic loci (nearly half of all loci in the data set), and smaller estimates of contemporary N e as compared with the RI-NCS and SCGAFL populations. Nonetheless, estimates for the GOM population reported here are larger than contemporary N e estimates based on nuclear markers that have been reported for other elasmobranchs, including the oviparous, seasonally aggregating, and threatened Zebra Shark Stegostoma fasciatum (N e = 377; Dudgeon and Ovenden 2015) and the near-threatened and migratory Spotted Eagle Ray Aetobatus narinari (N e = 1,893; Newby et al. 2014).
When specimens from the RI-NCS and SCGAFL populations were examined as a single metapopulation (RI-SCGAFL), N e estimates were smaller than the combination of each population's individual N e estimates. If the two populations conformed to an island model, where local populations do not exchange migrants, it is expected that the N e of the combined metapopulation would be greater than the sum of each population's independent N e (Waples 2010). Cases when an N e estimate for a metapopulation is similar to the N e estimate of a local population, as indicated by the current study, can be caused by weak genetic differentiation and high migration rates (Waples and Do 2010), the first of which has been demonstrated here. Calculations of N e assume discrete generations, no migration, random mating, no mutation, and no natural selection, and violations of these assumptions, as well as low sample size, can bias estimates of N e and result in "infinite" estimates (Waples and Do 2010;Waples et al. 2016;Marandel et al. 2019). The first and second assumptions are known to be violated; Clearnose Skate may have historic or contemporary migration between populations and have overlapping generations (sexually mature between age 3 and age 4 and live to be at least 6 years old ;Daiber 1960;Fitz and Daiber 1963).

Limitations and Future Work
More research is needed to verify the existence and characteristics of these Clearnose Skate populations. Population structure and biogeography could also be detailed using mitochondrial DNA markers (Avise et al. 1987a;Toews and Brelsford 2012). Population structure can be further assessed using a thorough morphological examination, the results of which could serve as diagnostic tools for identifying an individual skate's population of origin (Byrkjedal et al. 2007;Santos and Capellari 2009). Moreover, investigating whether the genetic patterns presented here correlate with environmental factors, such as bottom temperature, could either confirm or refute the posited hypothesis that Clearnose Skate demographics are partly driven by temperature preferences. In addition, tagging and tracking studies are recommended to clarify our understanding of individual Clearnose Skate migration patterns and dispersal capabilities. This would aid in estimating connectivity between the RI-NCS and SCGAFL populations, as well as inform appropriate sampling strategies (Marandel et al. 2018;Siskey et al. 2018). Finally, rangewide studies of Clearnose Skate life history, such as fecundity, age and growth, and mortality, would identify whether differences in these aspects exist throughout their range as well as ensure consistent methodology.
Specimens for the current study were collected opportunistically by fisheries-independent surveys. Although generally abundant, there may be differences in the relative abundance of Clearnose Skate throughout their range, but this remains unclear. Thus, sample sizes among the study regions and states were unequal and were unevenly distributed within them. Future specimen collection for molecular-based studies should strategically target areas that were underrepresented in this study to allow more robust comparisons. Equalizing the male to female ratio throughout each sample region would aid in determining whether genetic indices differ between sexes, that is, whether gene flow is mediated by one or both sexes, and potentially elucidate reproductive patterns. It is also important to sample the same areas throughout the year as this can address unknowns such as temporal and spatial variation of sex ratios or genetic structure related to seasonal separation or congregations for feeding or migration periods (Wearmouth and Sims 2010).
The results from this investigation of Clearnose Skate population structure can be used to better monitor and POPULATION STRUCTURE OF CLEARNOSE SKATE manage this vulnerable elasmobranch, thereby ensuring the sustainability of U.S. fisheries. Regional management council recognition that Clearnose Skate are not panmictic and that local depletion of one population may or may not be replenished by migrants from another population is imperative. The NEFMC, which maintains the only fishery management plan for Clearnose Skate, could incorporate population structure in developing more accurate estimates of abundance and catch limits. These findings could also be foundational to drafting management plans by the Southeast Fishery Management Council and Gulf States Fishery Management Council. Lastly, all councils should avoid incorporating existing life history literature in which specimens from a different region were examined into their management plans. Instead, councils should urge research into the life history traits of Clearnose Skate specific to their respective population (Siskey et al. 2018).