Sex matters: Otolith shape and genomic variation in deacon rockfish (Sebastes diaconus)

Abstract Little is known about intraspecific variation within the deacon rockfish (Sebastes diaconus), a recently described species found in the northeast Pacific Ocean. We investigated population structure among fish sampled from two nearshore reefs (Siletz Reef and Seal Rock) and one offshore site (Stonewall Bank) within a <50‐km2 area off the Oregon coast. Fish from the three sample sites exhibited small but statistically significant differences based on genetic variation at >15,000 neutral loci, whether analyzed independently or classified into nearshore and offshore groups. Male and females were readily distinguished using genetic data and 92 outlier loci were associated with sex, potentially indicating differential selection between males and females. Morphometric results indicated that there was significant secondary sexual dimorphism in otolith shape, but further sampling is required to disentangle potential confounding influence of age. This study is the first step toward understanding intraspecific variation within the deacon rockfish and the potential management implications. Since differentiation among the three sample sites was small, we consider the results to be suggestive of a single stock. However, future studies should evaluate how the stock is affected by differences in sex, age, and gene flow between the nearshore and offshore environments.

These species occur in sympatry from northern California to central Oregon; however, the deacon rockfish has a more northern distribution-extending to Vancouver Island, British Columbia, whereas the blue rockfish is more southern-reaching northern Baja California (Frable et al., 2015). Given the previous coupling of the two species, little is known about demographic, ecological, and genetic variation within deacon rockfish.
Previous studies investigating intraspecific variation in rockfish of the northeast Pacific have focused on the influence of the north-south latitudinal gradient on growth and maturity (Frey, Head, & Keller, 2015;Gertseva, Cope, & Matson, 2010). Clear breaks in population genetic structure have also been reported near oceanographic boundaries such as upwellings (Cope, 2004;Sivasundar & Palumbi, 2010). However, few studies have examined the influence of the east-west longitudinal gradient on intraspecific variation, which is closely related to depth change between the nearshore and offshore environments (Boehlert & Kappenman, 1980). Deacon rockfish inhabit a wide depth range, occurring from the shallow intertidal zone to depths >70 m (Frable et al., 2015 Figure 1). In 2017, fishing at depths ≥55 m opened to a new recreational fixed gear fishery that utilized modified terminal tackle (though deacon retention was prohibited until 2019) (Hannah, Buell, & Blume, 2008). Also in 2017, the markets for fish caught in the offshore midwater trawl fishery re-emerged and the fishery began to operate again. In the nearshore, recreational and commercial hook and line fishing for deacon rockfish never closed and significant catches of deacon rockfish continue to occur.
Unlike the nearshore environment, life history data for deacon rockfish in offshore areas are limited due to the previous fishing restrictions. It is important to determine whether the nearshore and offshore represent distinct fish stocks so that assessment models can account for potential connectivity between the two areas. The need to precisely define stock boundaries in the management of rockfish was demonstrated by a recent population genetic study of three species sampled from Puget Sound and outer coastal waters (Andrews et al., 2018). Based on the genetic results, canary rockfish Sebastes pinniger (Gill, 1864) from Puget Sound and the outer coastal area were concluded to represent a single genetic population (Andrews et al., 2018). This result suggested that the species did not meet the criteria for Endangered Species Act listing, which was previously designated based on evidence from other rockfish species (Andrews et al., 2018).

Deacon rockfish are managed by the Pacific Fishery Management
Council, the National Marine Fisheries Service, and each coastal state through measures such as annual catch limits for each stock or stock complex, harvest guidelines, trip or bag limits, area and gear restrictions, and seasonal closures. These measures are described in the Pacific Coast Groundfish Fishery Management Plan (PFMC, 2016). The most recent stock assessments combined deacon and blue rockfish (e.g., Dick et al., 2017), and it is uncertain how species differences influenced the stock assessment model. In Oregon, the stock boundaries were defined by the California state border to the south and the Washington border to the north.
The definition of a "fish stock" is ultimately a management decision (Carvalho & Hauser, 1994;Hilborn & Walters, 1992). Thus, for the purpose of this paper we follow Cadrin, Karr, and Mariani (2013) and define a fish stock as "an exploited fishery unit" that "may be a single spawning component, a biological population, a metapopulation, or comprise portions of these units. For management purposes, stocks are considered discrete units, and each stock can be exploited independently or catches can be assigned to the stock of origin." Since demography and genetic variation are typically interlinked, fish stocks are often considered genetic populations (Coyle, 1998;Ovenden, Berry, Welch, Buckworth, & Dichmont, 2015;Waldman, 1999). However, variables such as practical limits related to fishing and phenotypic traits (e.g., length and age at maturity) can be of equal importance from a management perspective, and interdisciplinary approaches are therefore necessary to define fish stocks (Abaunza, Murta, & Stransky, 2013;Cadrin & Secor, 2009;Coyle, 1998;MacLean & Evans, 1981;Ovenden et al., 2015). Delineating stocks is important for the sustainable management of fish populations as it allows researchers to investigate the influence and interaction of environmental and anthropogenic factors (Begg, Friedland, & Pearce, 1999;Cadrin & Secor, 2009).
The use of variation in the shape of anatomical structures (e.g., otoliths and scales) to identify fish stocks has been used since Lea's (1929) seminal work on herring. The underlying assumption is that the shape of the structure reflects environmental differences among potential populations. These methods have matured since the advent of image processing methods and the implementation of Fourier transformations to analyze the outline of structures (Bird, Eppler, & Checkley, 1986;Castonguay, Simard, & Gagnon, 1991).
The integration of genetic information into fish stock assessments has been relatively slow, primarily because traditional genetic markers such as microsatellites typically provide limited insight toward recent population genetic change, or local adaptation in marine organisms (Waples, Punt, & Cope, 2008). The advent of high-throughput sequencing methods has significantly increased the amount of data and the resolution of genetic insight for fisheries management in other species (Hauser & Carvalho, 2008;Kumar & Kocour, 2017;Riginos, Crandall, Liggins, Bongaerts, & Treml, 2016;Valenzuela-Quiñonez, 2016). Many studies have attempted to identify neutral and adaptive genetic variation (e.g., Gagnaire et al., 2015;Nielsen, Hemmer-Hansen, Foged Larsen, & Bekkevold, 2009;Ovenden et al., 2015;Valenzuela-Quiñonez, 2016), which has improved the delineation of populations and fish stocks in both migratory species such as Greenland halibut Reinhardtius hippoglossoides (Walbaum, 1792) (Westgaard et al., 2017) and European hake Merluccius merluccius (Linnaeus, 1758) (Milano et al., 2014), and sedentary species such as bluespotted Cornetfish Fistularia commersonii Rüppell, 1838 (Bernardi, Azzurro, Golani, & Miller, 2016). Typically, neutral genetic variation reflects stochastic genetic drift and the degree of gene flow among populations, whereas adaptive variation suggests selective differences among populations (Funk, McKay, Hohenlohe, & Allendorf, 2012). Adaptive genetic variation can reflect differential selection on certain genes among populations, despite the absence of obvious genetic differentiation for other markers. In addition to neutral and adaptive genetic differences, it is important to consider genetic variation associated with sex (Grummer et al., 2019). Such variation and sex biases in sampling can lead to inaccurate interpretations of population genetic differentiation (Benestan et al., 2017), potentially leading to incorrect management decisions.
The aim of this study was to use an interdisciplinary approach (Abaunza et al., 2013) to test for population structure and potential fish stocks among deacon rockfish off the Oregon coast based on otolith shape and genetic variation. We sampled fish from two nearshore reefs (Siletz Reef and Seal Rock) and one offshore area (Stonewall Bank) within a small geographic radius (<50 km 2 ). The three sample sites were analyzed independently, and differences between nearshore and offshore samples were tested as well. In order to assess the influence of sex in our analyses, we tested for otolith shape and genetic differences between males and females.
To disentangle any potential interaction between sample location and sex, we also analyzed genetic variation among the three sample sites using males and females separately.

| Sampling
Deacon rockfish were collected from three sites located off the Recreational hook and line gear was used for all collections. At each site, terminal gear included a variety of plastic baits, smallto medium-sized flies, and Sabiki rigs (herring jigs). Prior efforts to collect deacon rockfish off Oregon have shown that Sabiki rigs are capable of capturing a wide size range of individuals (~8-40 cm in this study), which helped offset gear-related bias in size selectivity of typical hook and line fishing gear (Ralston, 1990). At sea, fish were measured to total length and these data were used to ensure a wide range of size classes were sampled. No attempt was made to sex fish at sea. Fin clips for genetics were also taken at sea (see below for methods), and then, whole fish were placed on ice until later dissection of otoliths.
Otoliths were sampled from 676 fish, with 110, 172, and 394 specimens from Siletz Reef, Seal Rock, and Stonewall Bank, respectively (Table 1A). Sampling was conducted between December 2016 and November 2017 during favorable weather periods. At the selected sample areas, a total of 50 individuals were collected every month except January, June, and September (n = 9 per area).
Sampling efforts each month were mostly constrained to a 24-hr period, although low catch rates at Seal Rock meant that fish collected on August 8th and 16th 2017 were combined to achieve adequate sampling. Age was determined for all otolith samples using the break and burn method (Chilton & Beamish, 1982; Figure S1).
Funding allowed a total of 96 fish to be sampled for genetic analysis, with 25 and 23 nearshore specimens from Siletz Reef and Seal Rock, respectively, and 48 offshore individuals from Stonewall Bank (Table 1C). All genetic samples from Siletz Reef and Seal Rock were collected within a single sampling effort on October 4th and TA B L E 1 Deacon rockfish sampled for otolith shape and genetic analyses Note: Sampling is listed for each sample site (Siletz Reef, Seal Rock, and Stonewall Bank), for the tested nearshore and offshore groups, and for males (♂), females (♀), and individuals of unknown sex (U). The dagger sign ( †) indicates where sample sizes were reduced as the sex of a small number of individuals was unknown. A full list of specimens used for each RAD sequencing dataset is provided in Table S1, and a spreadsheet in the online supplement lists all otolith samples. The datasets are as follows: (A) Otolith dataset, used to analyze shape variation among the three sample sites (N = 676), between the nearshore and offshore groups (N = 676), and between males and females (n = 668). (B) Genetic dataset for variation among the three sample sites (n = 73). (C) Genetic dataset for variation between the nearshore and offshore groups (N = 96), and between males and females (n = 94). (D) Genetic dataset for variation among the three sample sites using only females (n = 50). (E) Genetic dataset for variation among the three sample sites using only males (n = 20).
November 6th, respectively. Fish from Stonewall Bank were collected in two even efforts on October 5th and November 6th.

| Otolith shape digitization and analysis
Fish were stored on ice for up to 24 hr, and otoliths were extracted using forceps after cutting the cranium. Otoliths were rinsed with freshwater, air-dried, and stored in binned otolith trays. All otoliths were used to investigate morphometric variation among the three sample sites and between the nearshore and offshore groups (N = 676; Table 1A). The sex of some smaller individuals (n = 8) was indeterminable; therefore, a slightly smaller dataset was used to estimate otolith shape differences between males and females (n = 668;  ; R Core Team, 2018) did not find any significant differences between the right and left otoliths of the sampled deacon rockfish, suggesting that differences observed among groups in this study were unlikely to be influenced by fluctuating asymmetry. A similar ANOVA method was used to estimate fluctuating asymmetry in an otolith shape analysis of lutjanid fishes (Vignon, 2015). Images were scale calibrated using fiji 1.51w (Schindelin et al., 2012).
Otolith shape was analyzed using shaper. The same method was previously used to distinguish two Sebastes species with high accuracy, but intraspecific variation was not investigated (Christensen et al., 2018). Contours were automatically extracted, and shape coefficients were estimated using a wavelet transformation. shaper implements both Fourier and wavelet transformations, and a comparison of the results from either transformation did not result in significantly different results. We decided to use the wavelet transformation. The shape coefficients were then standardized for fish length using the methods of Lleonart, Salat, and Torres (2000) that are implemented in the shaper package, with the aim of controlling for size and potential ontogenetic differences among fish of different ages. Variation among the potential populations was analyzed using a PERMANOVA test, using default settings with 1,000 permutations . No interactions were tested due to the fact that area and sex were confounding. Otolith shape variation among samples was visualized using a canonical analysis of principal coordinates (CAP; Anderson & Willis, 2003) in the vegan 2.5-2 R package (Oksanen et al., 2018).
The effect of sample size on the accuracy of the PERMANOVA test and CAP was investigated. We randomly subsampled, with replacement, the dataset of 676 otoliths 1,000 times and generated datasets increasing in sample size by multiples of 50 up to 500 (i.e., 50, 100, 150, … 500). Each dataset was divided evenly between sampling from the nearshore and offshore areas, and we tested the discrimination of those groups.

| DNA extraction and sequencing
A piece of caudal fin was taken from each fish and stored in a 5ml vial filled with 95% ethanol. Whole genomic DNA was extracted following the protocol and buffer solutions described by Ivanova, Dewaard, and Hebert (2006). DNA was quantified using Qubit highsensitivity dsDNA fluorometric quantitation (Life Technologies, Thermo Fisher Scientific Inc.).
All available DNA samples were used to analyze variation between the potential nearshore and offshore populations (N = 96; Table 1C), and males and females (n = 94; Table 1C). The three sample sites were analyzed independently with approximately the same number of samples per site (n = 73; Table 1B). Two further subsampled datasets were used to analyze variation among the three sample sites using only males (n = 20; Table 1D) and only females (n = 50;

| Processing of RAD sequencing data
A total of 342,702,104 DNA read pairs were sequenced (Table S2). to process reads, identify loci, and estimate genotypes. Forward and reverse reads from each index were demultiplexed into separate inline barcodes using the process_radtags component of the stacks pipeline. Simultaneously, the process_radtags step was used to remove reads with low-quality read data and ambiguous barcodes and RAD tags, resulting in a total of 291,066,365 read pairs being retained (Table S2). This step included the rescue barcode and RADtag parameter (-r) to retrieve additional reads. Only single-end (forward, R1) reads containing the SbfI restriction site were analyzed downstream.
Reads were assembled into stacks of similar DNA sequences and then into catalogs of reads for each investigated dataset using ustacks and cstacks. Following the recommendation of Mastretta-Yanes et al. (2015), many of the parameters in ustacks, cstacks, and populations were modified to examine the RAD sequencing data comprehensively (Table S3). In ustacks and cstacks, the default parameters (-m 2 -M 2 -N 4 -n 1) used in our final analyses were judged to provide a high number of stacks, consistent with other parametrizations (Table S3), with a low risk to introducing a high rate of erroneous reads. In ustacks, this meant that the minimum depth of coverage used to create a stack was two (-m 2), the maximum distance (in nucleotides) allowed between stacks was two (-M 2), and the maximum distance allowed to align secondary reads to primary stacks was four (-N 4). A bounded SNP model was applied with the error rate not being allowed to exceed 5% (--bound_high 0.05). In cstacks, the number of mismatches allowed between sample loci when building a catalog was one (-n 1). Locus coverage depth per individual was similar across the tested datasets, although some individuals yielded more loci with adequate coverage than others. As an example, Figure S2 presents mean (with standard deviation) and maximum coverage depth per individual for the dataset comparing nearshore and offshore fish (N = 96), where the overall mean coverage depth was 18.7 reads (SD 12.2, mean maximum 95.1).
Population genetic variation was estimated using the populations component of the stacks pipeline. In populations, the minimum stack depth required for individuals at a locus was set at five (-m 5). Samples were organized into multiple, independent datasets, which differed in the number of individuals and designated populations used to construct a loci catalog ( Table 2, Table S1).
The datasets were three independent sample sites, nearshore versus offshore, male versus female, female-only three sample sites, and male-only three sample sites. The minimum number of populations that a locus needed to be present in (-p) was set to the same number of populations set for each dataset (e.g., nearshore vs. offshore -p 2; three sample sites -p 3; sex -p 2). The minimum percentage of individuals in a population required to process a locus for a given population was set at 60% (-r 0.6). A minimum allele frequency of 5% was enforced for loci (--min_maf 0.05). Only the first SNP of each locus was included (--write_single_snp). All variant SNPs were biallelic.
Putative paralogous sequence variants (PSV) were identified using the python and R scripts for paralog-finder 1.0 (Mortiz, 2018), which is based on hdplot (McKinney, Waples, Seeb, & Seeb, 2017) and accounts for varying degrees of missing data per locus. Loci estimated to be in linkage disequilibrium (LD) were identified using plink 1.9 (Purcell et al., 2007). Putative PSVs and one locus for each loci pair estimated to be in LD were organized into a blacklist (-B; Catchen et al., 2013), and the populations component of stacks was rerun (same settings as above) so that these sites were removed from subsequent analyses ( Table 2). We tested for conformance to Hardy-Weinberg equilibrium (HWE) using vcftools 0.1.16 (Danecek et al., 2011; Table 2). This HWE estimation used an exact test (Wigginton, Cutler, & Abecasis, 2005), and we corrected for multiple testing by using a false discovery rate (FDR) adjustment for p-values with a critical threshold of <5% (Allendorf, Luikart, & Aitken, 2013;Bouaziz, Jeanmougin, & Geudj, 2012;Storey, 2002;Waples, 2015).

| Genetic variation
After removing loci estimated to be putative PSVs or in LD, we estimated observed (H O ) and expected heterozygosity (H E ) for each group tested in the separate loci datasets, using the adegenet 2.1.1 R package (Jombart, 2008;Jombart & Ahmed, 2011). Using the same R package, we also estimated allelic richness and the inbreeding coefficient (F IS ) for each tested group. The whoa 0.01 R package was used to investigate genotype frequencies, as well as the relationship between read depth per locus and heterozygote miscall rates (Anderson, 2018 Note: Loci estimated to be paralogous sequence variants (PSVs) and in linkage disequilibrium (LD) were removed. No loci were identified to be out of Hardy-Weinberg equilibrium (HWE). For an exact list of individuals included in each dataset, see Table S1.
Genetic variation among sampled individuals and groups was explored using principal components analysis (PCA), again implemented in adegenet (Jombart, 2008;Jombart & Ahmed, 2011). For PCA, outlier loci were not removed but rather all loci within each dataset were analyzed together. We determined the number of "meaningful" principal components (PCs) to retain for interpretation and downstream analyses by using the broken-stick test on PC Eigen values in the vegan R package. Retained PCs are axes that explain more variance among samples than expected by chance alone (Cangelosi & Goriely, 2007;Jackson, 1993).

| Outlier loci and genetic differentiation
Outlier loci were estimated independently using four genome scan programs: fsthet 1.01 (Flanagan & Jones, 2017a), outflank 0.2 (Whitlock & Lotterhos, 2015), bayescan 2.1 (Foll & Gaggiotti, 2008), and pcadapt 4.0.3 (Luu, Bazin, & Blum, 2017). Four separate programs were used because the stringency of outlier classification, false discovery rate, and the fit of applied models to particular patterns of genetic variation are known to vary among methods (see discussion Ahrens et al., 2018;Flanagan & Jones, 2017a;Luu et al., 2017;Narum & Hess, 2011). Outlier loci identified by fsthet, outflank, and bayescan are F ST outliers, which are sites that exhibit higher genetic differentiation among groups than expected by a neutral model (Ahrens et al., 2018;Foll & Gaggiotti, 2008). However, these programs use different statistical methods. fsthet analyzes the empirical relationship between F ST and observed heterozygosity (Flanagan & Jones, 2017a), whereas outflank analyzes the distribution of a special form of F ST that does not correct for sample size (Whitlock & Lotterhos, 2015). bayescan uses Bayesian maximum-likelihood to analyze differences in allelic frequencies among groups (Foll & Gaggiotti, 2008). In contrast, pcadapt does not consider F ST , and instead, loci are identified as outliers with respect to population structure among sampled individuals, using PCA (Luu et al., 2017).
Default settings were used for fsthet and outflank (Flanagan & Jones, 2017a;Whitlock & Lotterhos, 2015). In bayescan, we used default parameters and a prior of 100, with a q-value threshold of 0.05 (analogous to an FDR of 5%; Foll & Gaggiotti, 2008), and output data were investigated using the boa 1.1-8-2 R package (Smith, 2007 (Foll & Gaggiotti, 2008;Whitlock & Lotterhos, 2015). As a conservative precaution against Type I error, only loci identified as outliers by all four programs were organized into separate outlier loci datasets. Without information from annotated genome or selection studies, we interpret all loci as putatively adaptive and presumed neutral (Shafer et al., 2015), and it is not yet possible to exclude the potential influence of genetic incompatibilities or genetic surfing upon observed allelic frequencies (Bierne, Welch, Loire, Bonhomme, & David, 2011;Excoffier & Ray, 2008).
We examined genetic differentiation among the groups tested in all five loci datasets, including the comparison of the three sample sites, the nearshore and offshore groups, and males and females.
Neutral and outlier loci were analyzed independently using Weir and Cockerham's (1984) pairwise fixation index (F ST ), implemented in the stampp 1.5.1 R package using 5,000 bootstraps (Pembleton, Cogan, & Forster, 2013). We included an FDR adjustment for p-values with a critical threshold of <5% for the F ST p-values, using the same method as used for HWE estimation.
Genetic population structure was investigated among groups using Bayesian genotypic clustering in structure 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). We tested for up to five potential genotypic clusters among individuals (K = 1-5). For each value of K, five replications of the admixture model with independent allele frequencies were applied, with an MCMC length of 50,000 generations and a 10% burn-in. The optimal number of clusters was determined by examining estimates of mean K probability for a given value of K (Pritchard et al., 2000), and deltaK, the rate of change in logarithmic probability of the data (Evanno, Regnaut, & Goudet, 2005) implemented in structure harvester 0.6.94 (Earl & vonHoldt, 2012).

| Identity of outlier loci
The genomic identity of outlier loci was investigated by exporting the 143-bp FASTA format consensus sequence of each outlier locus from the stacks catalog and using NCBI blastn (Altschul, Gish, Miller, Myers, & Lipman, 1990)

| Otolith shape analysis
The PERMANOVA test found statistically significant differences in otolith shape between groups (nearshore and offshore) and among sites (Siletz, Seal Rock, and Stonewall Bank; and Stonewall Bank were similar when analyzed independently, or organized into nearshore and offshore groups. Using CAP to visualize shape variation among individuals, it appears that all sites could be distinguished, despite substantial overlap (Figure 2). In the CAP, otolith shape variation among the three sample sites was differentiated using two axes, whereas variation between the nearshore and offshore groups used only one axis ( Figure 2).
There was a significant otolith shape difference between males and females, with a larger pseudo F-value compared to the differences among sample sites (11.2; Table 3C). Graphs produced by CAP indicated that males and females exhibited distinctive distributions, but there was considerable overlap between them (Figure 2). Variation between males and females was restricted to a single axis (Figure 2).
In the investigation of potential sample size effects, probability density plots revealed that the central tendency was relatively consistent as sample size increased for both the nearshore and offshore groups ( Figure S3). However, CAP plots for the discrimination of the (c) Comparing males and females: males (n = 231) and females (n = 437). One CAP axis represented 100% of variation among individuals two groups showed that the distinction of the groups increased significantly as sample size increased from 50 to 350 individuals ( Figure   S4). After 350 individuals, however, increased sample size did not appear to have a significant effect upon discrimination ( Figure S4), and the mean value of the first CAP axis experienced only minor changes ( Figure S5). Results from the PERMANOVA test also indicated that the average F-statistic changed little after sampling ≥350 individuals ( Figure S5).

| Genetic variation
After controlling for multiple testing, no loci showed significant deviation from Hardy-Weinberg equilibrium ( Table 2). We removed loci likely to be PSVs or in LD from our datasets ( Table 2). The number of loci in each dataset was similar (~14,000-16,000), despite the varying number of individuals and groups (Table 2).
Across all five datasets, values for observed and expected heterozygosity were similar, with slightly fewer heterozygotes observed than expected (average difference across datasets of 0.02; Table 4).
This slight deficiency in heterozygotes was reflected with positive F IS values (Table 4). This result could be due to genotyping error or the sampling of a higher number of related individuals than expected by chance. Based on analyses in whoa (Anderson, 2018), genotypic frequencies did not appear to be significantly biased and an increased heterozygote miscall rate for loci with low read coverage is unlikely to have significantly affected the results (Figures S6-S8). According to the Wang relatedness estimator, however, we found no evidence for high relatedness (≥0.25) within and among the three sample sites using 73 individuals (Table S4). Allelic richness was also similar among groups tested in each dataset (Table 4).
No obvious population structure was revealed by PCA in comparisons of the three sample sites or the nearshore and offshore groups, and all PCs failed the broken-stick test-meaning that any patterns observed were likely to be a product of chance. However, in the comparisons of males and females, it was obvious that the first PC in each dataset reflected sex (Figure 3a,b). This suspicion was confirmed in the dataset comparing males and females (Figure 3c), where the separation across PC1 was clearly associated with sex. When genetic variation was examined among the three sample sites using only males or only females, PC1 (and all further PCs) for each dataset no longer exhibited any obvious structure ( Figure S9).

| Outlier loci and population genetic differentiation
No outlier loci were identified by all four genome scan programs (fsthet, outflank, bayescan, or pcadapt) when comparing the three sample sites or the nearshore and offshore groups (Table 5). In contrast, 92 outlier loci were identified by all four genome scan programs for the dataset comparing males and females (  Figure S10. Using presumed neutral loci, all pairwise F ST estimates were statistically significant for comparisons of the three sample sites and between the nearshore and offshore groups (  (Table S5).
For the female-only dataset comprising the three sample sites (n = 50), we found statistically significant differences for two out of three pairwise comparisons based on variation at the neutral loci (Table 6D) Table 6E). Although significance varied in these smaller datasets, F ST values indicated that the degree of genetic differentiation among males or females from the three sample sites was similar to the results in the three sample site dataset that included both sexes (Table 6).
Genotypic clustering results estimated by structure were similar to the PCA and pairwise F ST results for each dataset (Figure 4). The optimal number of clusters for the three sample site and nearshore versus offshore dataset was one or two (Figures S11 and S12), and the two clusters identified for both datasets separated most males and females (Figure 4). In the sex comparison dataset, the optimal number of clusters for the 92 outlier loci was two and the optimal number for the remaining neutral loci was one or two (Figures S11 and S12). The optimal number of clusters for the female-only and male-only three sample site datasets was one ( Figures S13 and S14).
Overall, the structure results indicated that males and females could be distinguished using both neutral and outlier loci, but that none of the sample sites across the remaining datasets could be distinguished, even when organized into nearshore and offshore groups.
F I G U R E 3 Scatterplots presenting genetic variation among deacon rockfish as estimated by principal components analysis (PCA), based on the SNP genotype data for all loci in each dataset. Plots show results for each dataset: (a) three sites, (b) two groups, and (c) sex. Sample sizes and the coloration of individuals are explained in the key, although it should be noted that PCA only visualizes the variance among samples. The first two PCs presented represent 2.1% and 1.7% of variation among individuals in (a), 1.7% and 1.3% of variation in (b), and 1.7% and 1.4% of variation in (c)

| Identity of outlier loci
We investigated the genomic identity of the 92 outlier loci detected for the sex comparison dataset. Some of the 92 outlier loci had an allele unique to males, but most loci instead exhibited an allele that was present in all 26 males but rare among females (e.g., 3/68 females were previously estimated as quantitative trait loci associated with sex determination and body growth (Loukovitis et al., 2015). Using BLASTN, our outlier loci did not appear to match any of the 26 malespecific, sex-linked loci identified in black-and-yellow and gopher rockfish by Fowler and Buonaccorsi (2016). Using the BWA-MEM algorithm, 91 of our 92 outliers successfully mapped to 86 unannotated scaffold sequences of the flag rockfish. Using BLAT, these scaffold sequences were mapped against the annotated genomes of three-spined stickleback, fugu, and Japanese medaka, with 62 aligning to chromosomes 2 and 3, 13 and 22, and 3 and 17 in each species, respectively.

| Concordance between otolith shape and genetic variation
The three sample sites (Siletz Reef, Seal Rock, and Stonewall Bank) could be distinguished using both otolith shape and genetic data when analyzed independently or organized into the nearshore and offshore groups. Although there was substantial overlap, statistically significant differences in otolith shape were found among the sample sites and the tested groups (Table 3, Figure 2).
Estimated pseudo F-values for the PERMANOVA test (Table 3) were similar to results reported in previous shaper otolith shape analyses of fish populations sampled over larger geographic distances (Berg et al., 2018;Lee, Brewin, Brickle, & Randhawa, 2018;Libungan, Slotte, Husebø, Godiksen, & Pálsson, 2015;Soeth et al., 2018). Additionally, a similar pattern of otolith shape variation was reported for two rockfish species sampled across the North Atlantic Ocean (Stransky, 2005). This comparability suggests that the differences in otolith shape observed for deacon rockfish are similar to those observed for populations spanning the entire North Atlantic Ocean, suggesting that they may be substantial. Resampling analyses also indicated that sample size was unlikely to have influenced the otolith shape results ( Figures   S3-S5). Previous otolith shape studies have sampled a large number (≥350) of individuals (e.g., Cañas, Stransky, Schlickeisen, Sampedro, & Fariña, 2012;Stransky, Murta, Schlickeisen, & Zimmermann, 2008b;Stransky, Naumann, et al., 2008a), but many others have sampled ≤100 individuals (e.g., Duncan, Brophy, & Arrizabalaga, 2018;Zhuang et al., 2015), which therefore may have had results affected by sample size effects. The analysis of potential sample size effects indicated that our total otolith sampling (676 specimens) was likely to represent accurate biological variation among the tested groups without significant sample size bias ( Figures S3-S5).
We investigated genetic differentiation using pairwise F ST and found that the three sample sites and the nearshore and offshore groups were significantly different based on neutral loci (Table 6).
However, the estimated genetic difference among potential groups was very low (pairwise F ST range of 0.0004-0.013; Table 6). Low  Signif. Signif.
* Note: Estimated 95% confidence intervals (in parentheses) and p-values are listed below each F ST estimate (using 5,000 bootstraps). Statistically significant results are marked with an asterisk (*). and the results of the power analysis conducted by Martinez et al. (2017), it is likely that our >15,000 loci (Table 5) provided sufficient statistical power to detect fine-scale population genetic structure, and that our genetic results are accurate.
Despite the statistical significance of our low F ST results, the broken-stick test did not identify any significant axes of variation in the geographic group datasets. There was no obvious geographical structure among samples when these nonsignificant PCs were visualized (Figure 3a,b). No outlier loci were identified in the datasets comparing the three samples sites or nearshore and offshore groups, despite analyzing >15,000 loci (Table 5). Similarly, genotypic clusters estimated by structure analysis did not align with geographic sampling (Figure 4). These results suggest that there are no substantial adaptive genetic differences among these sample sites or between the nearshore and offshore groups. The observed genetic variation based on neutral loci may reflect genetic drift and stochastic demographic changes (i.e., population size and migration rates).

| Influence of sex
Sex had an observable effect in our genetic datasets (Figures 3 and 4).
This result is consistent with previous reports that have stressed the importance of accounting for sex in reduced representation sequencing studies, where sex-linked variation can cause erroneous estimates of population genetic differentiation (Benestan et al., 2017;Catchen, 2017). We identified 92 outlier loci associated with sex (Table 5C).
Males and females were significantly different based on both the presumed neutral and outlier loci (Table 6C; Figure 4). Although significant, the F ST value based on neutral loci (0.0036) was low and comparable to results for the sample sites, whereas the F ST estimate based on the outlier loci was much higher (F ST 0.45; Table 6C). Since three of the applied genome scan programs identified F ST outliers, the higher F ST value estimated for outliers in the sex comparison dataset are expected. Removing 92 outlier loci associated with sex from the datasets comparing the three samples sites and the nearshore and offshore groups did not have a significant effect on F I G U R E 4 structure bar graphs estimated for genetic variation among deacon rockfish. Each vertical bar represents a separate individual in each dataset and the genotypic clusters estimated among individuals are shown in different colors (light gray and charcoal). The height of a cluster within each vertical bar indicates the confidence that a particular individual is assigned to a given genotypic cluster (referred to as the membership coefficient). The bar graphs show results for each dataset: (a) three sites, (b) two groups, and (c) sex. For sex, putatively adaptive and presumed neutral loci are separated. A colored bar above the graphs denotes the sex of individuals, with males (blue), females (purple), and unknown sex (orange), and bars below the graphs indicate group membership pairwise F ST results (Table S5). This finding indicates that genetic differentiation among the sample sites was therefore not solely driven by variation for these 92 outlier loci associated with sex.
Pairwise F ST estimates and PCA indicated that males and females exhibited a similar pattern of genetic variation among the three sample sites (Table 6D,E; Figure S9). The F ST estimate was statistically significant for the Siletz Reef and Stonewall Bank comparison in both the female-only and male-only datasets, and for Siletz Reef and Seal Rock comparison in the female-only dataset (Table 6D,E).
However, all other pairwise comparisons were not statistically significant. The lack of genetic differentiation between Seal Rock and Stonewall Bank in both the male and female-only datasets could be attributed to the shorter geographic distance between these sample sites (24 km) compared to the distance between each site and Siletz Reef (>40 km). On the other hand, the lack of differentiation between Siletz Reef and Seal Rock in the male-only dataset may be attributed to the small sample size (n = 20). In concordance, genotypic cluster analysis in structure also did not identify any population structure among the three sample sites using either sex-specific dataset. Altogether, these results suggest that genetic differences between male and females are unlikely to have influenced comparisons of the three sample sites.
Potential genetic differences between males and females were not examined in the previous RAD sequencing studies of other rockfish species (Andrews et al., 2018;Martinez et al., 2017). In the study of grass rockfish Sebastes rastrelliger (Jordan & Gilbert, 1880), sex was not recorded and PCA plots were not presented (Martinez et al., 2017). However, sex was recorded for canary rockfish and bocaccio Sebastes paucipinis (Ayres, 1854), and no obvious geographic structure was observed for either dataset using PCA (Andrews et al., 2018). In yelloweye rockfish, sex was again recorded, but the first two PCs reflected obvious geographic structure among sampling regions (Andrews et al., 2018). Comparison of these findings suggests that the genetic difference estimated between male and female deacon rockfish is relatively strong within the Sebastes genus, but the potential cause is unknown until further biological information is available for the species. Overall, our results indicate that variation associated with sex should be explicitly investigated in future population genetic studies of rockfish to avoid potential misinterpretation of data.

| Identity of outlier loci
The 92 outlier loci identified between males and females are associated with sex, but do not appear to be strictly sex-linked, as the outlier loci are variant positions shared between males and females. Some of the outlier loci had an allele exclusive to males, but conversely many loci instead exhibited an allele in all males as well as a small number of females. This bias in genetic variation suggests that there may be differential selection (intralocus sexual conflict) between males and females for autosomal or pseudoautosomal regions (segments of sex chromosomes that recombine) in the deacon rockfish genome.
Adaptive genetic differences can occur between the sexes if males and females are placed under different selection pressures for traits such as reproduction and behavior (Bonduriansky & Chenoweth, 2009;Cox & Calsbeek, 2009;Kasimatis, Nelson, & Phillips, 2017).
The genetic basis for such traits can be attributed to variation at a single locus (Bonduriansky & Chenoweth, 2009;Bonduriansky, Maklakov, Zajitschek, & Brooks, 2008;Cox & Calsbeek, 2009;Kasimatis et al., 2017;Mank, 2017;Parker & Partridge, 1998;Rowe, Chenoweth, & Agrawal, 2018), which can lead to high estimates of F ST (Flanagan & Jones, 2017b;Lucotte, Laurent, Heyer, Ségurel, & Toupance, 2016), as observed in this study (Table 6C). A recent RAD sequencing study of gulf pipefish Syngnathus scovelli (Evermann & Kendall, 1896) found that males typically possessed the minor allele, whereas females had the major (Flanagan & Jones, 2017b), which is the same pattern observed for most of the 92 deacon rockfish outlier loci in this study. Since the deacon rockfish was only recently discovered, there is currently insufficient biological information to hypothesize potential drives for selection between males and females, and this subject warrants future investigation. The sex determination system of deacon rockfish is currently unknown. Previous karyotype and genetic research on other rockfish species have indicated both XY and ZW heterogametic systems within the genus (Anderson, 1979;Fowler & Buonaccorsi, 2016;Ida, Iwasawa, & Kamitori, 1982). An alternative hypothesis for these 92 outlier loci is that they occur on an X chromosome and that male deacon rockfish are the heterogametic sex with XY sex chromosomes.
Under this scenario, males would have a single X chromosome copy of these loci that may have been misinterpreted as homozygous sites compared to the same loci in females where there are two X chromosome copies and potential for heterozygosity. If true, we would expect individual males to be consistently homozygous for either allele observed in females. Instead, most loci were heterozygous for members of both sexes, and some loci exhibited an allele exclusive to males or an allele that were frequent for males but rare among females. Given this pattern, it seems more likely that intralocus sexual conflict is the cause. Future investigations should investigate the influence of heterogametic loci when identifying loci and biallelic SNPs in RAD sequencing programs such as stacks.
If the 92 outlier loci occur within a pseudoautosomal region of potential sex chromosomes in deacon rockfish, we could expect most of the outliers to occur within the same genomic region. We tested this hypothesis by modifying the approach of Fowler and Buonaccorsi (2016). All but one of our 92 outlier loci successfully mapped to 86 unannotated scaffold sequences of the flag rockfish, and in turn, 62 of these flag rockfish scaffold sequences aligned to chromosomes 2 and 3, 13 and 22, and 3 and 7 in three-spined stickleback, fugu, and medaka, respectively. In contrast, most scaffolds associated with the 26 male-specific, sex-linked loci for black-andyellow and gopher rockfish, identified by Fowler and Buonaccorsi (2016), mapped instead to chromosomes 14, 12, and 6 in each respective BLAT species. These results suggest that our outlier loci occur on two chromosomes, and not on the equivalent chromosome to the Y chromosome identified in black-and-yellow and gopher rockfish (Fowler & Buonaccorsi, 2016).
It is therefore uncertain whether the outliers represent two autosomal regions of the deacon rockfish genome strongly associated with differential selection between the sexes, or if the outliers are pseudoautosomal and deacon rockfish exhibit different sex chromosomes to black-and-yellow and gopher rockfish. In addition, the small but significant F ST value for the presumed neutral loci between males and females suggests that the sexes may differ for many other positions throughout the genome. Since the karyotypes of other rockfish species have indicated both XY and ZW sex determination systems (Anderson, 1979;Ida et al., 1982), and the last common ancestor of the black-and-yellow, gopher, and the blue rockfish is estimated to have occurred ~6.2 Mya (Hyde & Vetter, 2007), it is possible that deacon rockfish have evolved a different set of sex chromosomes or use an alternative sex determination system altogether (e.g., temperature). Intriguingly, a high intensity of intralocus sexual conflict within loci may drive gene duplication and the evolution of new sex determination systems (Gallach & Betrán, 2011;van Doorn, 2009), so the observed 92 outlier loci in deacon rockfish may reflect the ongoing evolution of sex chromosomes in Sebastes rockfish (Fowler & Buonaccorsi, 2016).
If true, a variable number of outlier loci associated with sex should be observed among rockfish species with different sex determination systems. Differential selection between males and females and change in sex determination systems could also be related to the high level of speciation in Sebastes rockfish (Gavrilets, 2014;Parker & Partridge, 1998

| Secondary sexual dimorphism in otolith shape
There was a significant difference in male and female otolith shape when all 660 available samples were analyzed together, which was an order of magnitude higher than the results comparing the three sample sites, or the nearshore and offshore groups (  (Bird et al., 1986;Castonguay et al., 1991). In contrast, although otolith shape variation in Atlantic cod Gadus morhua Linneaus, 1758, was mostly influenced by differences in growth rates among populations, a small but significant shape difference was observed between males and females (Campana & Casselman, 1993). A recent study by Parmentier, Boistel, Bahri, Plenevaux, and Schwarnzhans (2018) reported substantial secondary sexual dimorphism in the hearing apparatus and otolith shape of the ophidiid Neobythites gilli (Goode and Bean, 1885). Since males and females of this species demonstrated similar hearing ability, it was hypothesized that differences in habitat preference (and associated environmental variables) were responsible for the observed dimorphism (Parmentier et al., 2018).
A potential cause for secondary sexual dimorphism in the otolith shape of deacon Rockfish is uncertain, although otolith shape variation in other rockfish species has been associated with differences in growth rates, habitat usage, and hormone levels (Tuset et al., 2015(Tuset et al., , 2016. No previous studies have compared otolith morphology and population genetics in rockfish, but Tuset et al. (2016) noted that otolith morphology is more strongly influenced by ecological and biogeographical factors rather than phylogeny. As in other fish lineages (Gauldie & Nelson, 1990;Lombarte & Lleonart, 1993), otolith length and width increases with age in Sebastes rockfish, but once fish stop growing in body size, otoliths barely grow in length and instead increase in thickness (Love et al., 2002). Like other rockfish, female deacon rockfish appear to reach larger body sizes than males (Hannah et al., 2015;Love et al., 2002). This difference could cause male otoliths to grow thicker for a greater proportion of their lifespan, generating a relative difference in the shape of male and female otoliths.
A future study with adequate sampling should be able to determine whether age has an influence on the otolith shape difference between males and females. Unfortunately, there is an age bias in our current otolith sampling, with most representatives of the oldest age classes originating from the offshore site of Stonewall Bank ( Figure S1). This bias means that we cannot disentangle the sample site variation from the potential influence of age. However, since the sex ratio for the nearshore and offshore groups was similar (1.87 and 1.96, respectively), it is unlikely that the otolith shape difference between these groups was influenced by sex. Ultimately, the linkage of otolith shape to biogeography and ecological variation suggests otolith shape could be a useful tool for stock discrimination in the genus.

| Significance of results and implications for fisheries management
Our otolith shape and genetic results indicate a small difference between two potential stocks of deacon rockfish in the nearshore and offshore, which corresponds with the current de facto management for the species. Regardless of whether deacon rockfish were organized into nearshore and offshore groups, morphological and genetic differences were statistically significant but small among the sample sites. Although our morphometric and genetic results were comparable to findings from other marine fishes sampled over larger geographic distances (Benestan et al., 2017;Berg et al., 2018;Lee et al., 2018;Soeth et al., 2018), previous stock assessment using similar methods has relied upon stronger patterns in data to delineate a stock boundary (Siegle, Taylor, Miller, Withler, & Yamanaka, 2013;Ward, 2000).
Although differentiation was low, the fact that we detected statistically significant otolith shape and genotypic differences over such a small geographic scale (<50 km 2 ) seems remarkable. This statistical significance may reflect the large amount of information (>15,000 loci) provided by the RAD sequencing method. Further genetic sampling across the range of deacon rockfish may help to interpret the scale and significance of variation observed in this study, and improve the distinction of potential stocks. Future genetic studies of rockfish should record the sex of samples and take into account the sex-ratio and age distribution of sample groups. Inattention to such factors may have already biased the results of previous genetic studies on rockfish (Waples et al., 2008). Future studies of deacon rockfish could investigate the effect of environmental variation on otolith shape across sample sites, as well as differences in habitat usage and feeding preferences. The differences observed for both otolith shape and genetic loci identified in this study may reflect demographic variation within deacon rockfish, which also warrants further investigation. It should be noted that the most recent stock assessment of deacon rockfish did not consider differences in life history parameters between fish from the nearshore and offshore (Dick et al., 2017). However, it should also be remembered that genetic and demographic variation are not always concordant (see discussion by Lowe & Allendorf, 2010;Waples et al., 2008).
Nonetheless, given that we ultimately regard a "fish stock" as a practical management decision (Carvalho & Hauser, 1994;Hilborn & Walters, 1992), it may be prudent to treat fish from the nearshore and offshore as distinct fish stocks until the future genetic, environmental, and demographic research is conducted. Deacon rockfish occurring at depths deeper than 55 m have been effectively not been harvested since 2003, and therefore, fish from the nearshore and offshore areas are essentially managed as separate stocks. Managing the nearshore and offshore fish separately is therefore unlikely to change the current status of the biological populations involved, and this seems to be the most conservative approach until further information is available. We suggest that future assessments and management decisions consider the ramifications of managing this species as one or two stocks. Future research sampling across the entire range of the deacon rockfish, and sequencing DNA from a larger number of individuals, may provide statistical power to differentiate potential nearshore and offshore stocks.

| Conclusions
We found small but statistically significant otolith shape and genetic differences among deacon rockfish sampled off the Oregon coast, regardless of whether the three sample sites were analyzed independently or organized into nearshore and offshore groups. We suggest that deacon rockfish from the nearshore and offshore are managed separately until further genetic, environmental, and demographic data are available, requiring no practical change from current management practices.
Sex mattered in our otolith shape and genetic analyses. We found evidence for secondary sexual dimorphism in otolith shape, which may reflect differences in the growth, age, and lifespan of males and females. Males and females were readily distinguished in our genetic data, although this is unlikely to have affected comparisons of the sample sites. Our results concur with previous studies that sex should be considered in population genetic research (Benestan et al., 2017;Catchen, 2017), particularly for Sebastes species.
We identified 92 outlier loci that are associated with sex in deacon rockfish. These sites likely reflect differential selection between males and females, which should be investigated in other rockfish species with potentially different sex determination systems. A possible biological cause for this selective difference is uncertain in deacon rockfish, due to the recent discovery of the species. This subject warrants further investigation, as the genetic variation may reflect differences between males and females for habitat usage, which in turn could result in different management requirements.
The data generated in this study can contribute to future, more extensive studies of Sebastes rockfish diversity. The sequence data are compatible with reads from previous RADseq studies of other rockfish species that also used the SbfI restriction enzyme (Andrews et al., 2018;Martinez et al., 2017). The shaper otolith digitization method  easily allows morphometric datasets from other regions or species to be combined as well. In particular, further research is needed to investigate biological and management requirement differences between deacon and blue rockfish.
This study provides a first step toward the investigation of intraspecific variation in a recently described species, and the results emphasize the potential of RAD sequencing to provide substantial population and sexual genetic information for species that have not been previously studied.

ACK N OWLED G M ENTS
The We are grateful to Dr. Andres Aguilar and Dr. John Fields for providing information, additional specimens, and for genetically testing the species identity of all small individuals. We thank Oregon State University's Center for Genome Research and Biocomputing for conducting the DNA sequencing, and for providing assistance with their computational infrastructure. Lastly, we are grateful for the feedback on this manuscript provided by the Editor-in-Chief Gareth Jenkins and two anonymous reviewers.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
Felix Vaux conducted the genetic analysis and wrote much of the manuscript. Leif K. Rasmuson designed the study, collected samples, conducted the shape analysis, and wrote much of the manuscript. Lisa A. Kautzi processed and photographed the otoliths