Skate egg nursery areas support genetic diversity of Alaska and Aleutian skates in the Bering Sea

Skate egg case nursery sites are specific locations on the ocean floor where some species of skates deposit egg cases to incubate for up to 5 yr until hatching. We examined genetic diversity within and among skate egg nursery sites of the Alaska skate Bathyraja parmifera and the Aleutian skate B. aleutica in the eastern Bering Sea to gain a better understanding of how skates utilize these areas. Restriction-site associated DNA (RAD) sequencing libraries were used to obtain single nucleotide polymorphism (SNP) datasets for B. parmifera (5285 SNPs) and B. aleutica (3309 SNPs). We found evidence for significant genetic differentiation among all B. parmifera and B. aleutica nursery areas, with 1 exception. B. parmifera from the spatially proximate Pribilof and Bering Canyons were genetically similar, suggesting that this may represent a large contiguous nursery area. Genetic differences between embryos at distinct developmental stages within nursery areas were not significant. Adult B. parmifera taken from the Bering Sea and Aleutian Islands were genetically distinct from embryo collections, implying that additional genetic types of B. parmifera may exist that were not represented by the nursery areas sampled in this study. Our data also showed evidence for low effective population sizes and low Ne/N ratios. Results indicate that nursery areas support genetically distinct components of each species, which underscores the importance of skate egg nursery areas for conservation of genetic diversity.


INTRODUCTION
Egg case nursery sites are localized areas where oviparous skates in the Bering Sea and Aleutian Islands lay egg cases in high densities (> 50 000 eggs km −2 ) in contact with the ocean bottom or other permanent structure to incubate for up to 5 yr until hatching without parental care (Hoff 2008). Skate egg and juvenile nursery areas fall within the broader definition of elasmobranch nursery areas, i.e. regions of the ocean where elasmobranchs are encountered in greater numbers than in other regions, may remain or return over extended periods, and use the habitat repeatedly across years (Heupel et al. 2007, Martins et al. 2018. Nursery areas are critical to the survival of many elasmobranch species, although the mechanisms and extent to which each species depends on nursery areas vary (Beck et al. 2001, Hoff 2016, Heupel et al. 2019. For example, nursery areas may be more important for larger and slower-growing elasmobranchs that produce fewer offspring (Heupel et al. 2019).
The reproductive behavior of skates associated with their nursery areas is not fully understood. Whether nursery sites are reproductively isolated or whether there is genetic connectivity among nursery sites is not known. It has also been suggested that some proportion of skates may deposit egg cases outside of nursery areas (Matta 2015). Understanding the genetic importance of nursery sites and the genetic stock structure of skate species is essential for best management practices and for determining how to protect essential habitat for skates (Waples 1998, Spies & Punt 2015. Like other elasmobranchs, skates grow slowly, mature late, and have low fecundity. The Alaska skate Bathyraja parmifera is found throughout the Bering Sea to southeastern Alaska and comprises > 90% of the skate biomass, estimated at approximately 600 000 t in 2020 (Ormseth 2020) on the eastern Bering Sea shelf. The Aleutian skate B. aleutica is found from northern Japan to the Aleutian Islands and southeast Alaska, and is the most abundant skate on the eastern Bering Sea slope, representing approximately onethird of the skate biomass in that region (Stevenson et al. 2008, Hoff 2016, Ormseth 2020. Fertilization in all skates found in the Bering Sea is internal, and embryos are protected by tough proteinaceous egg cases laid by females in pairs on the seafloor (Matta 2006, Hoff 2009. The age at 50% maturity for B. parmifera has been estimated at 9 yr and fecundity is estimated between 21 and 37 eggs yr −1 (Gunderson 1997, Matta 2006. The age-structured model used for the stock assessment of B. parmifera extends through age 25, as the maximum age observed is at least 20 yr old (Ormseth 2020). Ageing is difficult for this species, and a higher maximum age is expected based on ob served population dynamics data.
In the eastern Bering Sea, 26 nursery sites have been identified for 6 species of skates, and many of these sites are used by multiple skate species (Hoff 2008, Rooper et al. 2019. These nursery sites are primarily found along the interface between the continental shelf and slope, where environmental conditions such as depth, current, speed, temperature, and aspect of the slope relative to current direction are favorable (Rooper et al. 2019). For these species of skates, egg case nurseries are directly tied to survival in their early life history, and represent specific habitat needed for developing embryos within egg cases during a protracted development time of at least 3 or 4 yr for some species (Hoff 2008, Rooper et al. 2019. The high density of egg cases present at nursery sites may be important as a 'predator swamping' survival strategy (Hoff 2009).
Multiple genetically diverse nursery areas may help preserve the biodiversity of skates if climate change or other factors reduce the population size (Anderson et al. 2015). Biological diversity stabilizes ecosystem processes and maintains resilience across populations. According to the 'portfolio' theory, when sharks inhabit a range of habitats at early life stages, it may help to reduce variation in adult population sizes (Yates et al. 2012). Similarly, diversity across hundreds of individual populations of sockeye salmon Oncorhynchus nerka reduces variability in salmon runs to Bristol Bay, Alaska (Schindler et al. 2010).
The US North Pacific Fishery Management Council has designated 6 skate nursery sites as Habitat Areas of Particular Concern (HAPC; NOAA 2015), but this listing confers no particular protection to skate egg case nursery areas in the Bering Sea. There is no directed fishery for skates, but adults are caught incidentally in substantial numbers, particularly in the Pacific cod Gadus macrocephalus hook and line fishery (Ormseth 2020). Trawl fisheries pose risks to skate nursery areas because of their potential to disrupt habitat on the ocean floor and physically remove or damage skate eggs. Climate change is also a potential risk factor; growth and development times for skate embryos are directly related to temperature (Hoff 2008). Nursery sites occur within a narrow range of temperature with low variability, but projections indicate that temperature and temperature variability will increase for the Bering Sea shelf in the near future (Hermann et al. 2016).
Our objectives in this study were to evaluate genetic diversity in the 2 most common skates in the eastern Bering Sea, B. parmifera and B. aleutica, within and among egg case nursery sites. We also examined a small collection of adults to examine how they were related to egg case nursery areas. Four of the 5 nursery sites we examined are currently listed as HAPC sites (Fig. 1), and a general description of the substrate and temperature observed is presented to provide context for habitat use. We generated 2 single nucleotide polymorphism (SNP) datasets based on low-coverage sequencing of restriction site-associated DNA sequencing (RADseq) libraries (Andrews et al. 2016). We used these SNP datasets to test the following 2 null hypotheses: (1) skate embryos collected from different nursery sites represent a single panmictic population, and (2) relatedness is the same within and among nursery areas. The first hypothesis was designed to test whether genetically distinct nurseries exist to maintain genetic diversity in the species, while the second hypothesis provided information on site fidelity and population structure.

Sample collection and habitat
Embryos (within skate egg cases) were collected from each nursery site in 2006, 2007, and 2016 using an 83−112 eastern otter trawl net similar to that used for Alaska Fisheries Science Center (AFSC) standard eastern Bering Sea bottom trawl surveys (Lauth & Acuna 2007). Embryos of Bathy raja parmifera were collected from the following sites: Bering Canyon, Pervenets Canyon, Pribilof Canyon, and between Pribilof and Zhemchug Can yons (IPZ1 and IPZ2) ( Fig. 1, Table 1). Bering Can yon was sampled twice, in 2007 and 2016, to exa mine the degree of site fidelity (Table 1). Embryos of B. aleutica were collected from Bering and Pervenets Canyons in 2016 (Table 1). Trawls were as short as possible (10 min) to minimize im pact to the nursery area. All embryos were preserved in 95% ethanol, and embryo total length was recorded after preservation and converted to developmental stage (2−5) following the methodology of Hoff (2009). An equal number of male and female embryos were collected for each cohort of stages 3−5 (stage 2 embryos were not differentiable by sex). Information on substrate type and temperature were collected at each nursery area. Adults were sampled opportunistically during AFSC bottom trawl surveys conducted from June through August in 2002 and 2006 in the Aleutian Islands and Bering Sea: 8 adults were obtained from Pervenets Canyon, 3 from IPZ2, 3 from the eastern Aleutian Islands, and 1 from the central Aleutian Islands ( Fig. 1; Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m669 p121 _ supp. xlsx). These samples consisted of a piece of tissue preserved in 95% ethanol for genetic analysis, supported by vouchered whole specimens or a photograph for specimens discarded at sea. Embryos from unique combinations of nursery site and year are referred to as collections (Table 1).

Data processing and quality checking
Data were downloaded from the GC3F website and processed using Stacks v1.46 (Catchen et al. 2013). Raw data were archived at the Sequence Read Archive (SRA), which is accessible from the National Center for Biotechnology Information (NCBI) at www. ncbi.nlm.nih.gov/bioproject/4306953, BioProject ID PRJNA646172. Data from embryos and adults of B. parmifera were processed separately from data from B. aleutica to maximize the number of SNPs for each species. Raw data were quality filtered and demultiplexed using the 'process_radtags' subroutine. All sequences were trimmed to 95 bp to accommodate the difference in sequence lengths (100 and 150 bp), and a population map ('popmap') was used to ensure that loci were equally represented in data sequenced on the HiSeq 2500 and the HiSeq 4000.
The 2 most closely related elasmobranch genomes available for alignment at the time of analysis were the whale shark Rhincodon typus (Read et al. 2017, Hara et al. 2018) and the white shark Carcharodon carcharias (Marra et al. 2019). The white shark ge -nome was selected for alignment because it contained fewer missing metazoan genes than the whale shark and fewer scaffolds (19 000); fewer scaffolds indicate a better understanding of chromosome structure (Marra et al. 2019). Alignment success of the B. aleutica data to the white shark was too low (< 0.5%) to acquire an aligned dataset.
A 2-step process was performed to establish a de novo baseline set of SNPs without duplicates, and effectively combine data run on different platforms. A de novo analysis pipeline was selected that increases the number of loci and reduces the SNP calling error rates (Mastretta-Yanes et al. 2015). The 10 individuals with the most sequence data from each collection were used in catalog construction during the first and second steps to reduce the detection of false polymorphisms. First, sequence alignment, SNP identification, and construction of a catalog containing all SNPs were performed using the subroutines 'ustacks,' 'cstacks,' 'sstacks,' and populations with default values except −m 3 (minimum depth of coverage required to create a stack). Default parameters options included −M 2, which was se lected to avoid large, repetitive sequences and paralogs (Mastretta-Yanes et al. 2015), where −M is the maximum distance in nucleotides allowed between stacks.
The 'fasta' file of all selected loci was converted to a database using Bowtie software v. 2.3.4.1 (Langmead & Salzberg 2012). The Bowtie database was aligned to itself to check for duplicate loci or loci that were inadvertently split (−v 3, −k 2), where v is the number of mismatches regardless of quality, and k is the number of valid alignments to report. Identical loci were culled so that all loci in the database were unique using the function 'RemoveDups.R' (https:// github.com/31ingrid/RemoveDups), which was run in R version 4.0.3 (R Core Team 2020). In the second round of analysis, data from all individuals were aligned to the new database using 'pstacks' with flags (−m 3, −bound_high 0.05, −max_locus_stacks 3), where bound_high is the upper bound for the error parameter epsilon, and max_locus_stacks is the maximum number of stacks allowed per single locus. These were followed by subroutines 'cstacks' (−g) and 'sstacks.' Finally, the subroutine 'populations' was run with no restrictions on the number of SNPs per locus (−m 10, min_maf 0.05), where min_maf is the minimum minor allele frequency. The r = 0.5 para meter was also incorporated, which required that 50% of individuals in a sample collection were genotyped to process a locus for that population.
The SNP with the highest minor allele frequency was then selected from each sequence read. Further filtering took place in PLINK software v1.90b5.3 format (Purcell et al. 2007) to identify loci with missing call rates, individuals with missing data, individual heterozygosity, and individual observed and expected homozygous genotype counts (−mind, −geno, −het). For datasets of both species, loci genotyped in less than 80% of individuals (−geno 0.20) were removed, be cause this provided sufficient loci genotyped in a high proportion of individuals (Table  S2). Individuals with less than 60% genotyping rate were re moved (−mind 0.4), which also balanced sequencing success while retaining as many samples as possible (Table S2). Individuals with anomalously high hetero zygosity were removed if they appeared to be low F IS outliers (<first quantile − 1.5× the interquartile range), as this was indicative of contamination (mixture of more than a single individual). Read depth and heterozygosity were plotted for all samples of B. parmifera to check for inconsistencies among sequencing platforms ( Fig. A1 in the Appendix), and heterozygosity by collection is shown as a series of boxplots for comparison (Fig. A2).
Prior to examining whether loci conformed to expectations of Hardy-Weinberg equilibrium (HWE) in B. parmifera, pairwise F ST and Jost's D were calculated and used to establish a prior grouping of collections. This step was important for the B. parmifera data because information on population structure was unknown, and uninformed removal of SNP loci due to HWE could bias marker selection. Prior clustering was not performed for the B. aleutica data be cause only 2 geographically distant nursery areas were sampled, and the expectation of Hardy-Weinberg proportions in each collection was biologically reasonable. The R package 'StAMPP' (Pembleton et al. 2013) was used to perform pairwise comparisons of F ST among all collections, calculated according to Weir & Cockerham (1984) with 1000 bootstraps across loci (Table S3). Pairwise Jost's D was also calculated among collections to provide a second prior estimate of population structure using 'MMOD' in R (Table S3; Winter 2012). Jost's D uses a multiplicative partitioning of diversity based on the effective number of alleles rather than expected heterozygosity (Meirmans & Hedrick 2011) and is considered a useful measure for examining allelic differentiation, particularly when mutation rates are low (Jost 2008). We tested for hierarchical genetic clustering using a method specifically tailored to groups with low levels of differentiation such as marine species (D'Aloia et al. 2020). Clusters were evaluated using the 'pvclust' package with the pairwise matrices of F ST and Jost's D, with Euclidean distance, average cluster method, and support for each cluster was based on 10 000 bootstrap replicates (Suzuki & Shimodaira 2006). Probabilities were presented as approximately unbiased values (AU) and bootstrap probabilities (BP); AU is considered more accurate than BP (Suzuki & Shimodaira 2006).
Once clusters were established, conformation of genotype proportions to HWE expectations was examined for each cluster of collections using 'hw.test' in the R package 'pegas' (Paradis 2010). A chi-squared test compared expected genotype frequencies calculated from allelic frequencies to actual genotypes. The false discovery rate (FDR) method (Benjamini & Hochberg 1995) was used to correct for multiple comparisons at an alpha level of 0.05. Loci were discarded if they were out of HWE in more than 1 out of 6 clusters identified for B. parmifera (following cluster analysis) and if they were out of HWE in both nursery site sample collections of B. aleutica.

Candidate loci for diversifying selection
Candidate loci under natural selection were examined using 2 methods, OutFLANK and BayeScan, with an FDR of 0.05. OutFLANK (Whitlock & Lotterhos 2015; https://github.com/whitlock/OutFLANK) estimates the distribution of F ST estimates based on a trimmed set of F ST estimates, and assumes that most loci in the center of the distribution are neutral; loci that fall outside the distribution are considered candidates for selection. The lower and upper 0.05 ends of the distribution were trimmed. BayeScan 2.1 was also used to assess the presence of loci under selection (Foll & Gaggiotti 2008). BayeScan utilizes differences in allele frequencies between populations in a Bayesian framework and provides F ST estimates with a q-value, the FDR analog of the p-value. This used the following collections: Bering 2016, Pribilof 2016, IPZ1, IPZ2, and Adults. Loci were considered as candidates for selection if they were identified as significant using both methods. A second test served as a check for scoring data compatibility in B. parmifera between the HiSeq 2500 and HiSeq 4000 platforms; all collections run on the HiSeq2500 were compared with collections run on the HiSeq4000. Loci were removed that were considered under selection among samples analyzed on different platforms.
Candidate loci under selection ( Fig. 3) were aligned to the Amblyraja radiata annotated genome, which consisted of 49 chromosomes and was published subsequent to analysis. This species is more similar to the genus Bathyraja than sharks, as both are classified in the order Rajiformes. Alignment to the A. radiata genome was performed using BLASTn (https:// blast.ncbi.nlm.nih.gov).

Inbreeding coefficient, basic statistics
The inbreeding coefficient, F IS , was calculated for each skate collection using the function 'genedivFis' in the R package 'genepop' (Rousset 2008). F IS was estimated for each collection overall and separately for each embryonic developmental stage. Significance of F IS was determined using χ 2 = n × F IS 2 , which follows the chi-squared distribution with a single degree of freedom, and n represents the number of individuals (Nei 1987). Rarefied allelic richness was calculated over all loci using the R package 'hierfstat' and reported as the mean over loci by population (Goudet 2005).

Tests for genetic differentiation
Pairwise F ST and Jost's D were performed on B. parmifera (for a second time based on the final dataset) and B. aleutica. Analysis of molecular variance (AMOVA) between nursery collections (Table 1), between stages within collections, between individuals within stages, and within collections was calculated using the R packages 'poppr' and 'ade4' (Kamvar et al. 2014, Dray & Dufour 2007. AMOVA testing was performed on embryo collections of B. parmifera performed in 2 ways. First, 5 groups were tested (N = 5), retaining Pribilof and Bering Canyons as distinct groups: IPZ1, IPZ2, Pervenets Canyon, Pribilof Canyon, and Bering Canyon (years combined). Second, 4 groups were tested (N = 4), by combining Pribilof and Bering Canyon collections: IPZ1, IPZ2, Pervenets Canyon, and Pribilof/Bering Canyons (years combined). All stages and collections of embryos of B. parmi fera were included, and collections from Bering Canyon were combined. Collections of B. aleu tica consisted of a single stage, so they were excluded from this analysis. The observed variance component (σ) at all hierarchical levels, and the pvalue to determine significance were calculated.
Isolation-by-distance was evaluated among 2016 embryo collections from Bering Canyon, Pribilof Canyon, IPZ1, and IPZ2. Other embryo collections taken from Bering Canyon in 2006 and Pervenets Can yon in 2007 were excluded due to small sample sizes. Distances among nursery areas were calculated as the shortest geographic distance, estimated using Google Earth (earth.google.com/web). A Mantel test with 1000 permutations was applied to the linearized matrix of pairwise F ST and geographical distance matrix.
The R package 'adegenet' was also used to conduct discriminant analysis of principal components (DAPC) to visualize the genetic relationship among collections for B. parmifera, with and without the adult collection (Jombart 2008, Jombart & Ahmed 2011. Cross validation ('xvalDAPC') was used to identify a reduced number of principal components that discriminated observations reliably and consistently.

Kinship
Two methods were used to determine relatedness, that of the software KING 2.1.3 (Manichaikul et al. 2010) and COLONY (Jones & Wang 2010). COLONY v.2.0.6.5 considers sibships among all sampled individuals jointly, while KING considers the relationship for a single pair of individuals in isolation. If 2 or more siblings (sharing the same father, or mother, or both) were present in a collection, COLONY used information (about their inferred common parent or parents) to find probability of full or half siblings (Jones & Wang 2010). The number of SNPs used in COLONY was reduced to 231 to increase computation time by selecting a minor allele frequency of 40% or greater (-maf 0.4) and a genotyping rate of at least 95% (-geno 0.05). COLONY was run with polygamy for males and females, medium precision for full likelihood, and no inbreeding. KING estimated 2 statistics: the proportion of SNPs with zero identity-bystate (IBS0) between a pair of individuals (e.g. AA and GG) and the kinship coefficient. Higher kinship coefficients indicate greater relatedness, and a negative kinship coefficient indicates no relatedness. All statistics were estimated within and among collections of B. parmifera and B. aleutica using full SNP datasets. Kinship coefficients estimated by KING and sibship probabilities from COLONY among nursery sites were used as a threshold probability for sibship within nursery sites, as the level of genetic differentiation among nursery areas renders the possibility of sibship among nursery areas ex treme ly unlikely, except between Bering and Pribilof Canyons.

Assignment testing
A mixed stock analysis in RUBIAS was performed to assign adult B. parmifera to nursery area of origin (Moran & Anderson 2019). RUBIAS is a Bayesian hierarchical genetic stock assignment approach that accounts for population structure among baseline populations using a mixture set of unknown origin and a reference set of known origin (Anderson et al. 2008). Posterior density curves and 95% credible intervals for mixture proportions were created from Markov chain Monte Carlo output. The accuracy of the individual assignment analysis was evaluated by testing self-assignment of simulated individuals of known origin. Known simulated proportions for each reporting group were compared with the numbers estimated by RUBIAS to test the accuracy of the individual assignment analysis. A Z-score was computed from each individual's log-likelihood and compared with a simulated normal density to test whether adults originated from nursery areas outside the reference set, and a Kolmogorov-Smirnov test was used to compare the Z-score with the normal density. RUBIAS required unscored loci to be removed from individual samples in mixture sets of unknown origin, which required the removal of 12 loci.

Effective population size
We used the program LDNe (Waples & Do 2008) to estimate effective population sizes, with a random mating model and a jackknife method for estimating confidence intervals. This method uses the amount of linkage disequilibrium within a population to estimate effective population size and corrects for bias due to sample size. We selected a minor allele frequency of 0.05 to ensure that the results were not influenced by the presence of rare alleles. Collections of B. parmifera from Bering and Pribilof Canyons taken in 2016 were combined for this analysis.

Sample collection and habitat
The substrate at Bering Canyon was soft sand, while the substrate at nursery site IPZ1 was boulders and sand. The substrate at the other nursery sites (Pribilof Canyon, IPZ2, and Pervenets Canyon) was sandy with bits of small to medium gravel. Observed temperatures at nursery areas fell within a range of <1°C (3.8−4.5°C), with the warmest temperatures found at the lowest latitude within Bering Canyon ( Table 1). Embryos of Bathyraja parmifera ranged in size from 23 to 241 mm total length, and all embryo collections contained at least 2 developmental stages, ranging from Stage 2 through Stage 5 (

Data quality filtering in B. parmifera
A total of 34 048 SNPs were identified following the initial stacks analysis of B. parmifera. Two individuals were removed in the stacks pipeline due to insufficient data or missing barcodes. For the other 405 individuals, 33 342 SNPs remained following selection of the single SNP with the highest minimum allele frequency where multiple SNPs were identified within the same sequence read. After filtering for missing genotypes, there were 5289 SNPs (Table S2). During filtering for individuals with < 60% sequencing success, 73 individuals were re moved (Table S2). Two additional samples were re moved: one was a duplicated sample (iden-  Table S1). The read depth and heterozygosity for individuals run on the HiSeq2500 fell in the middle of the range of values seen in individuals sequenced on the HiSeq4000 (Fig. A1), indicating read depth and heterozygosity were not skewed or deficient for the Pervenets 2006 and Bering 2007 collections. Samples sequenced by the HiSeq2500 had a read depth of 64.9 on average and mean heterozygosity of 0.34, while samples sequenced on the HiSeq4000 had an average read depth of 90.8 and mean heterozygosity of 0.31. OutFLANK identified 3 candidate loci under positive selection in the dataset for B. parmifera due to differences in scoring among platforms (Table S4). These loci were removed, leaving 5285 loci in the dataset used to test for clustering.
Estimates of pairwise F ST (prior to filtering for HWE) indicated that embryos of B. parmifera from IPZ1 and IPZ2 were significantly differentiated (the latter highly so) from all other nursery areas except Bering 2007 and Pervenets 2006 (Table S3). F ST between IPZ2 vs. Bering, Pribilof, and IPZ1 collections was 0.0027, 0.0027, and 0.0020 respectively, and all were significant at alpha = 0.05 following FDR pvalue correction (Table S3). F ST between IPZ1 vs. Bering, Pribilof, and IPZ2 were 0.0003, 0.0002, and 0.002, respectively, and were all significant as well (Table S3) (Table S3).
Dendrograms provided > 95% support for adult and IPZ2 collections as distinct from other groups, while Bering 2016 and Pribilof 2016 collections clustered together with > 95% support (Fig. 2). IPZ1 clustered separately from other nursery areas with 100% support based on . AU is considered more accurate than BP, as the latter is a simple statistic computed by bootstrap resampling.

Data quality filtering in B. aleutica
Stacks identified 41 009 SNPs in the dataset for B. aleutica at a stack depth of 10 or greater following data filtering and specified constraints, and 25 227 re mained after retaining the locus with the highest minimum allele frequency per read. There were 23 903 loci remaining after HWE filtering, and 3309 SNPs remained in the final dataset after filtering for missing genotypes. PLINK identified 10 individuals with > 60% missing data, which were pruned from the dataset (Table S2). Three individuals were removed due to identical genotypes identified by KING. This produced a final dataset for B. aleutica with 84 individuals, 3309 SNP loci, 89.6% genotyping success, and a mean stack depth of 77 (Table 1).

Candidate loci for diversifying selection
After removing loci that were out of HWE, 5 loci were identified as candidates for positive selection by OutFLANK (Table S4). BayeScan identified 3 candidate loci under positive selection, 1 of which was also identified by OutFLANK (69495_52, Fig. 3; Tables S4 & S5). The most extreme allele frequencies of the 3 SNPs identified by BayeScan (including the one also identified by OutFLANK) occurred in the Adult and IPZ2 collections (Table S5). OutFLANK identified 26 loci as candidates for positive selection in the dataset for B. aleutica. One of these was also the only locus identified by BayeScan (Table S4). Up to 10 matches to the Amblyraja radiata genome were presented for each short read, with percent identity > 80%, as well as with E-values <10 −14 , and gene name and function (Table S6). Two short reads (25 904) in B. parmifera and sequence 162 908 in B. aleutica produced many matches, while 69 495 had no matches and 34 342 produced only one (Table S6).

Inbreeding coefficient, basic statistics
None of the F IS estimates was considered significant following FDR correction for multiple tests (Benjamini & Hochberg 1995). For B. parmifera, the inbreeding coefficient, F IS , was positive in collections from all nursery areas except Pervenets (Table 2). This trend did not change when individual stages were considered separately. The highest F IS was ob served in the adult samples (0.3657, Table 2), indicating more admixture than relatedness among samples, although these values were not significantly different from zero. The same trends in F IS were reflected by heterozygosity; adult B. parmifera and IPZ2 had notably lower levels of heterozygosity than other collections of B. parmifera (Table 2; Fig. A2). Similarly, the lowest rarefied allelic richness, A R , was ob served in the adult collection (A R =1. 6747), followed by the IPZ2 collection (A R = 1.6817) ( Table 2). F IS was also positive (although not significantly different from zero) in B. aleutica collections, and in the range observed among B. parmifera embryo collections, 0.0874 for Bering Canyon and 0.0566 for Pervenets Canyon (Table 2). Allelic richness was slightly higher for these 2 collections than B. parmifera, both greater than 1.99, but not comparable with B. parmi fera because the datasets were analyzed with separate pipelines.

Tests for genetic differentiation
Pairwise F ST and Jost's D were calculated a second time with the final dataset containing 4802 SNP loci after filtering for HWE (Table 3). Pairwise tests also included 2016 Pribilof and Bering collections combined, which showed significant differences with all other collections except Bering 2007 and Pervenets 2006. Given their small sample sizes, these collections may not be representative of the nursery sites from which they were taken. We also calculated F ST among cohorts with 1000 bootstrap replicates from the same nursery area (Table S7), and none was significant among stages 2, 4, and 5 from Bering Can yon (2016), among stages 2, 3, and 5 from Pribilof Can yon (2016), or among stages 2, 3, 4, and 5 from IPZ1 (2016).
AMOVA results consistently showed significant variation among embryo collections of B. parmifera from different nursery areas. The pvalues decreased from 0.023 to 0.007 when Bering and Pribilof samples were combined into a single group, and the percentage of total variance attributed to differences among collections increased from 0.08 to 0.12, providing support for grouping samples from Bering and Pribilof Canyons. In neither case was significant variation observed among stages within collections (Table 4). Higher-level comparisons, among individuals within stages and within collections were not significant in either case (p = 1.000).
Adult B. parmifera were distinct from embryos in the DAPC, in which 200 principal components were selected using cross-validation (Fig. 4a). A single adult (Bparm209) taken from Pervenets Canyon clustered with embryos from the Pervenets, Bering, and Pribilof Canyon nursery areas (Fig. 4a). The remaining adults, including 7 from Pervenets Canyon, 3 from IPZ2, 3 from the eastern Aleutian Islands, and 1 adult from the central Aleutian Islands, were located in a separate half of the plot, indicating relatively  Table 4. Results of AMOVA testing on embryo collections of Bathyraja parmifera performed 2 ways: (a) 5 groups (N = 5): Inter-Pribilof-Zhemchug 1 and 2 (IPZ1, IPZ2), Pervenets Canyon, Pribilof Canyon, and Bering Canyon (years combined), and (b) 4 groups (N = 4): IPZ1, IPZ2, Pervenets Canyon, and Pribilof combined with Bering Canyon (years combined). The results include the observed variance component (σ) at several hierarchical levels, the percent of the total variance (%), and the p-value to determine significance, based on 1000 permutations of the data large differentiation between adult and embryo samples. Removing adults from the DAPC plot ( Fig. 4b;  100 principal components) highlighted differentiation between IPZ1 and IPZ2, as well as a bifurcation in the samples from IPZ2, with some clustering with embryos from other sites and others situated distally. Embryo collections from the other canyons (Bering, Pervenets, and Pribilof) clustered together, including temporally distinct samples from Bering Canyon. The test for isolation-by-distance was not significant (p = 0.309). This was in part due to high levels of differentiation between the IPZ2 nursery area collection and all other collections, including embryos from the closest nursery area, IPZ1, 134 km away (Table S8).

Kinship
Kinship coefficients for B. parmifera were negatively correlated with IBS0, meaning that higher kinship was associated with higher probability of identity-by-state (Fig. 5). The range and mean of kinship coefficients differed among B. parmifera and B. aleutica collections (Fig. 5). The highest mean kinship coefficient for B. parmifera em bryos was among Pervenets Can yon samples (0.004), followed by Bering 2007 (−0.013) and then Bering/ Pribilof 2016 combined (−0.032) (Fig. 5a). The lowest mean kinship coefficient for B. parmi fera was among adult individuals (−0.786), consistent with random sampling of that collection. The IPZ1 and IPZ2 collections had the lowest kinship coefficients of all nursery areas (−0.182 and −0.319). Mean relatedness among all B. parmi fera nursery areas (−0.286) was lower than mean relatedness within all nursery areas except IPZ2, indicating some level of admixture in that sample.
The mean kinship coefficients for embryos of B. aleutica were also negatively correlated with IBS0 and were within the range of values observed for B. parmifera: −0.270 for Bering Can yon and −0.065 for Pervenets Canyon (Fig. 5b). Relatedness between B. aleutica from Bering and Pervenets Can yons was lower than relatedness within either nursery area (−0.355) (Fig. 5b).
The highest kinship coefficient calculated with KING 2.1.3 among nursery collections was as high as the highest kinship coefficient within nursery areas (0.15). A similar result was found using COLONY; the highest probability of full siblings among nursery areas was 0.703, and the probability of full siblings within nursery areas was 0.703 or lower. Given the large population size of the B.  Table S1 in the Supplement). IPZ: Inter-Pribilof-Zhemchug parmi fera population in the Bering Sea, we ruled out the possibility that re lated individuals were sampled from different nursery areas. Therefore, sibship was considered un likely within any of the B. parmifera nursery areas because no pairwise comparisons were higher within nursery areas than among nursery areas.

Assignment testing
The Kolmogorov-Smirnov test was used to examine a mixture dataset including all adults (except Bparm209) against a reference dataset that included 3 groups: IPZ1, IPZ2, and a third group containing embryos from Bering, Pervenets, and Pribilof Can - yons. The Kolmogorov-Smirnov test was significant (p < 0.0001), providing evidence that the adults in the mixture set likely originated from nursery area(s) outside the reference set. However, simulated re-assignment success was low for the adult samples and for IPZ2 (Fig. 6). A second test assigned adult Bparm209 to a reference dataset that included IPZ1, IPZ2, embryos from Bering, Pribilof, and Pervenets combined, and all adults except itself. Bparm209 was assigned with highest probability to the Bering-Pribilof-Pervenets embryo group, 79.1%, with a credible interval of 0.187−1.000 (Fig. 7, Table 5). Low assignment success of simulated adult and IPZ2 samples implies uncertainty in this result.

Effective population size
Effective population size was estimable among col-  ture of skate egg nursery areas. Egg case nursery areas are essential habitat for eggs to develop for up to 5 yr before hatching, but whether females deposit eggs randomly or in specific nursery areas is unknown. These questions are important for understanding the conservation value of skate egg case nursery areas and for determining appropriate management units. This represents the first study to examine the genetic connectivity among skate egg nursery areas in the Bering Sea. We examined thousands of SNPs to understand whether genetic diversity differed significantly among the nursery areas and to compare relatedness within and among nursery areas. Our most important finding was that skate egg nursery sites in the eastern Bering Sea do not represent a single panmictic population. Embryos of B. parmifera from IPZ1 and particularly IPZ2 were genetically distinct from other collections, and differences were not associated with geographic distance or isolation-bydistance (Table S8). In contrast, embryos of B. parmifera from the Bering and Pribilof Canyon nursery sites sampled 229 km apart appeared to represent a large interconnected region for egg case deposition by B. parmifera. We also observed significant genetic differences among embryos of B. aleutica from Pervenets and Bering Canyons (840 km apart). Embryos of different developmental stages within nursery sites were genetically similar, in contrast to significant differences in F ST observed among several egg case nursery sites (Table 3, Fig. 4; Table S7). The F ST values among nursery sites were relatively low in magnitude (0.000−0.003) at distances from 130− 530 km, but low F ST values may still be statistically significant and biologically meaningfu l (e.g. Waples 1998, Knutsen et al. 2011. Genetic differentiation among adults of other species of skates and rays was higher than observed here. Chevolot et al. (2006) observed genetic differentiation among adult thornback rays Raja c lavata (global estimate of F ST = 0.042) in the eastern North Atlantic between the Mediterranean, the Azores, and continental shelf of Europe at distances of 500−1300 km. O'Connell et al. (2019) found evidence for fine-scale population structure (F ST = 0.00−0.12) in winter skates Leucoraja ocellata but not in little skates L. erinacea between Georges Bank and the Mid-Atlantic Ocean at distances on the order of 300 km.
Another key finding was that relatedness within the Bering, Pribilof, and Pervenets Canyon B. parmifera collections and within the B. aleutica collections was higher than among embryo collections, indicating that in general, skate embryos taken from the same nursery area were more related than skate embryos taken from different nursery areas (Fig. 5). Higher relatedness within nursery areas indicates some level of site fidelity, at least by females, which could be further tested through analysis of mitochondrial DNA. While relatedness was typically higher within than among nursery collections, an exception was found at IPZ2, which had the lowest mean relatedness. The IPZ2 nursery area also displayed the highest F IS and largest spread of heterozygosity of all embryo nursery collections, and a diffuse pattern in the DAPC plot (Fig. 4), suggesting mixed usage by unrelated groups of B. p armifera (Table 2; Fig. A2).
While we predicted that adult samples would be closely rel ated to the sampled embry os, we had the unexpected result that only 1 of the 15 adults of B. parmifera taken from the Bering Sea and Aleutian Islands was genetically similar to embryos from sampled egg case nursery sites (Fig. 4). This adult male, which was caught in a single haul with 4 other males to which it was not related, had high assignment probability to the Bering and Pribilof Canyon nursery areas (Figs. 6 & 7; Table S1). It should be noted that although the probability of reassignment was high (80%) for embryos from Bering and Pribilof Canyons, it was much lower for adult samples and embryos from IPZ1 and IPZ2 (Fig. 7), indicating that assignment test results may be unreliable and larger sample sizes may be needed. The genetic distinctiveness of the adults (F ST : 0.005−0.009) in our collection indicated that the nursery areas we sampled are not a complete representation of the genetic diversity in this species. Adult samples were less related to each other than were embryos within collections from different nursery areas (Fig. 5), which is consistent with random sampling over a large area. Adult samples were also associated with positive F IS (Table 2), which can be the result of low heterozygosity, potentially resulting from the Wahlund effect or subpopulation structure. Whether the genetic distinctiveness of the adult skates examined could be due to longdistance migration from natal habitat is not known. Limited dispersal appears to be the rule in L. erinacea, L. ocellata, and big skates Beringraja bino culata, with a few individuals performing long-distance migrations that may be related with mating or egg deposition (King & McFarlane 2010, Farrugia et al. 2016, O'Connell et al. 2019. A tagging study of 1238 R. binoculata recovered 75% of individuals within 21 km of their tagging location, and only a small proportion of skates (1.5%) performed long-range migrations up to 2350 km (King & McFarlane 2010). These long-range migrations were primarily per-formed by females that were maturing or had recently matured (King & McFarlane 2010).
Finally, results indicated low effective population sizes and ratios of effective population size to cen sus size (N e /N ): for B. parmifera, N e = 1257 at IPZ1 and 890 at Bering/Pribilof nursery areas, and for B. aleutica, N e = 3722 in Pervenets Canyon and 236 in Bering Canyon. These estimates of N e strongly suggest low N e /N ratios in B. parmifera given the estimate of over 200 million individuals of this species in the Bering Sea, although precise estimates of N e /N ratios are not directly calculable (Ormseth 2020). The N e /N ratio is typically small in marine fish populations (10 −3 or lower), whereas sharks exhibit high N e /N ratios, approaching 1 (Portnoy et al. 2009, Dudgeon & Ovenden 2015. Our results are supported by estimates of low N e /N ratio in thornback rays, between 9 × 10 −5 and 6 × 10 −4 (Chevolot et al. 2008). It is not understood whether the causes of the low N e /N ratios in thornback rays and species of Bathyraja exa mined here result from skewed sex ratios, fluctuations in population size, or variance in reproductive success (Chevolot et al. 2008). Complex reproductive patterns may be a likely explanation, resulting in unequal sex ratios or high variance in reproductive success; polyandry has been observed in manta rays Manta birostris and southern stingrays Dasyatis americana (Yano et al. 1999, Chapman et al. 2003, Chevolot et al. 2008. The presence of multiple male participants appears to be a conserved behavior during courtship of bato ids (Mc-Callister et al. 2020). High juvenile mortality due to protracted egg incubation time could also provide an explanation for low N e /N ratios. The effective population size (N e ) is an important parameter for conservation because it represents the size of an idealized population with the same genetic diversity as the population of interest. Higher N e /N ratios are characteristic of more vulnerable populations, as removal of fewer individuals reduces the overall effective population size. Effective population sizes greater than 500 will prevent loss of genetic variation, and N e values lower than 1000 may result in accumulation of deleterious alleles (Franklin 1980, Lande 1994, Dudgeon & Ovenden 2015. Although B. parmifera appear to have a low N e /N ratio, the effective population sizes by nursery area for B. parmifera and B. aleutica are also on the low end of the threshold for conservation. Similarly, low effective population sizes were observed in thorn back rays in the Irish Sea (N e = 283) (Chevolot et al. 2008). Samples consisting of multiple cohorts are not expected to have a large effect on the effective population size because they are genetically similar; however, mixed groups at IPZ2 may down-wardly bias estimates of effective population size (Waples et al. 2014). It should be noted that small sample sizes are a potential reason for failure to estimate effective population size, and may provide a negative bias in effective population size; therefore, future efforts should use larger sample sizes (Marandel et al. 2019).
Future work is needed to resolve new and unanswered questions resulting from this study. Although we identified candidate markers under selection for both species, further work is needed to understand whether there is a genetic component to local adaptation in these species. Also, sequence read depth was over the 40× coverage threshold considered sufficient to minimize errors in SNP calling (Catchen et al. 2011), but the lack of resolution in the relationship between the Bering 2007 and Pervenets 2006 collections was apparent in clustering dendrograms (Fig. 2) and could be a result of low sample sizes (20 and 18, respectively). While this study answered several questions about the genetic population structure of B. parmifera and B. aleutica, it raised other questions. The presence of comingling genetically distinct adult skates remains a mystery, as does their nursery area of origin (Fig. 4a). Tagging studies of these species would greatly improve our understanding of skate behavior, migration, and reproduction. This study also highlighted the need to understand causes for high variability in reproductive success in these species leading to low effective population size estimates and remarkably low N e /N ratios. Finally, additional sampling from all B. parmifera nursery sites will help us estimate the sizes of populations served by these nursery sites.
All egg case nursery sites that we sampled are currently listed as HAPCs, except IPZ1 (NOAA 2015). Given the genetic distinctiveness of embryos from IPZ1 and its low effective population size, we suggest that this site also merits consideration as an HAPC. Further, our research indicates that strengthening this designation to include protection from mobile bottom fishing gear and its known impacts on seafloor habitat and ecosystems is important for all egg case nursery areas (Kaiser et al. 2002, Thrush & Dayton 2002, Jennings et al. 2005. Commercial trawling within skate egg case nursery areas could potentially destroy large numbers of eggs that are crushed either during retention in the net or by the trawl footrope. For those egg cases that survive the mechanical aspects of being captured as bycatch and discarded during trawl operations, they will likely be displaced from their nursery area and may have increased mortality outside those specific habitats.

CONCLUSIONS
This study revealed new information on the population structure and use of nursery habitats by Alaska skates Bathyraja parmifera and Aleutian skates B. aleutica. Data showed that skate egg nursery sites in the eastern Bering Sea do not represent a single panmictic population, and that egg case deposition by females is non-random. Adult samples indicated that more nursery areas may exist, although our data did not provide evidence for the number or nature of those sites. Genetic differentiation did not follow an isolation-by-distance pattern. Nursery sites were genetically distinct, with the exception of the adjacent Bering and Pribilof Canyons, and may support relatively low effective population sizes, on the order of 1000. Results of this study also indicated very low N e /N ratios in skates of the genus Bathyraja, similar to a study of thornback rays, which may be attributable to complex mating behavior. Multiple genetically diverse nursery areas may work in tandem to maintain genetic diversity in these species, similar to a meta population complex or the portfolio effect. Further complexity in samples from IPZ2 suggests that several distinct groups may utilize this nursery area. The results of this work justify the importance of the conservation of nursery areas as essential habitat for maintaining genetic diversity in B. parmifera and B. aleutica, and may serve to inform future studies of the population dynamics of these and other skates in the Bering Sea.