Introduction

The genus Sebastes (Scorpaeniformes:Sebastidae) is a rapidly evolving lineage of marine fishes that includes many commercially and recreationally important species. Although Sebastes can be found in temperate regions of the North Atlantic and southern hemisphere, the greatest diversity is found in the North Pacific (Hyde and Vetter 2007; Love et al. 2002). In the North Pacific, the genus has been characterized as an adaptive radiation with rapid diversification into heterogeneous ecological niches and habitats that often occur along a depth gradient (Heras and Aguilar 2019; Hyde and Vetter 2007; Ingram 2011; Johns and Avise 1998; Love et al. 2002). In such rapidly diverging lineages, the true extent of species diversity can be underestimated as some closely related taxa may be difficult to differentiate morphologically. Indeed, genetic studies within Sebastes have uncovered several cryptic species pairs including S. aleutianus and S. melanostictus (Gharrett et al. 2005; Orr and Hawkins 2008), S. mystinus and S. diaconus (Burford and Bernardi 2008; Frable et al. 2015) and S. miniatus type I and type II (Hyde et al. 2008a, b; Hyde and Vetter 2007). The identification of cryptic species complexes such as these, which clearly represent distinct taxonomic units, is critical to effective fisheries management. Although closely related, sister taxa can often be characterized by dissimilar life history traits, abundances and productivity. For instance, in the rougheye rockfish species complex, S. aleutianus and S. melanostictus exhibit significantly different reproductive parameters (Conrath 2017). These parameters are vital for accurate stock assessment models and the estimation of stock biomass. In cases such as this, treating distinct taxa as a single complex may present significant risk to one or both species and place them at risk of overfishing (Burford et al. 2011; Frable et al. 2015; Hyde and Vetter 2009).

The vermilion rockfish species complex, which includes the vermilion rockfish (S. miniatus) and the yet to be formally described sunset rockfish (S. miniatus type I sensu Hyde et al., 2008a, b), is of particular concern as it is one of the most valuable recreational fisheries on the U.S. West Coast. These sister species, which diverged around 2.3 Ma (Hyde and Vetter 2007), show a high degree of range overlap from central California to northern Baja, Mexico but exhibit dissimilar depth distribution as adults. Vermilion rockfish are more common in shallower waters (< 100 m) in kelp forest habitat while sunset rockfish are typically found deeper (> 100 m) at offshore banks (Hyde et al. 2008a, b; Love and Passarelli 2020). The vermilion complex is the third-most commonly landed recreational species of groundfish on the U.S. West Coast with approximately 625 mt of landings in 2019 (RecFIN extracted 1 Feb 2021) and forms a critical component of the coast’s recreational fishery that contributed over $1.8 billion to the region’s GDP based on the most recent estimate prepared in 2014 (Lovell et al. 2016). The stock status of both sister species are not fully known as previous attempts to assess the assemblage in 2005 and 2013 were not endorsed for setting harvest guidelines by the Pacific Fishery Management Council (PFMC) in part because of uncertainty surrounding the relative contribution of each species within existing data sources.

NOAA’s Northwest Fisheries Science Center (NWFSC) annually conducts a bottom trawl survey and a hook and line survey along the U.S. West Coast that, along with multiple other sources, has acquired thousands of tissue samples nominally identified as vermilion rockfish with associated biological and environmental data (Harms et al. 2008; Keller et al. 2017). Providing a reliable method to disentangle sunset and vermilion rockfish will allow for the development of separate abundance indices, life history profiles, and, potentially, catch histories to support the assessment of this commercially and recreationally important species complex.Footnote 1

Previous efforts to delineate the vermilion rockfish complex to species using mitochondrial DNA (mtDNA) alone yielded inconclusive results due to a historical one-way mitochondrial introgression of vermilion rockfish-type haplotypes into sunset rockfish populations (Hyde et al. 2008a, b). Specifically, a portion of sunset rockfish are characterized by vermilion rockfish mitochondrial haplotypes and sunset rockfish nuclear genotypes, and are referred to here as introgressed sunset rockfish. Although microsatellite markers have proven effective for disentangling this species complex (Hyde and Budrick, unpublished data), this approach remains costly and is not well suited to the scale of samples that must be genotyped (~ 25,000) to achieve a suitable dataset for stock assessment purposes. An attractive alternative is the use of a high-throughput approach, such as restriction site-associated DNA sequencing (RADseq), that can identify thousands of genome-wide markers and potentially diagnostic single nucleotide polymorphisms (SNPs). With proper primer design, diagnostic SNPs can be incorporated into an assay, such as a Genotyping-in-Thousands by sequencing (GTseq) panel (Campbell et al. 2015), that provides a cost-effective assignment method (Larson et al. 2014; McKinney et al. 2020). Such an approach in the vermilion rockfish complex would allow for species-specific demographic and biological analyses to be carried out with the thousands of samples already collected. Additionally, RADseq has proven to be an effective tool in identifying and assessing intraspecific differentiation, especially in cases where traditional markers (i.e., mtDNA and microsatellites) failed to detect signals of structure (Bohling et al. 2019; Gaither et al. 2020; Longo et al. 2020; Morgan et al. 2018; Vaux et al. 2021; Vendrami et al. 2019). Although distinguishing between clear taxonomic units at the level of species is paramount, the identification of significant intraspecific population genetic structure is also important for species-specific management (Bernatchez et al. 2017; Hauser and Seeb 2008; Spies and Punt 2015; Waples et al. 2008). Furthermore, an understanding of population genetic structure allows for insight into the ecological and evolutionary processes shaping differentiation (Bradbury et al. 2008; DeWoody and Avise 2000; Martinez et al. 2018; Romiguier et al. 2014). Previous population genetic work on vermilion rockfish revealed significant population structure and detailed how a shift in ontogenetic migratory behavior in conjunction with oceanographic conditions likely resulted in the observed differentiation (Budrick 2016; Hyde and Vetter 2009). Past population genetic work on sunset rockfish, however, did not detect intraspecific differentiation although the study used a limited number of sampling sites (Budrick 2016). Clearer insight into the population genetic structure of sunset and vermilion rockfishes, respectively, will allow for more informed, species-specific stock assessment and management strategies, and may further the understanding of the mechanisms driving diversification in the Sebastes species flock.

In this study we used RADseq markers to assess the degree of interspecific differentiation between sunset and vermilion rockfishes, and also evaluated and compared population differentiation within each species. With these loci we also identified diagnostic SNPs for the future development of an assay capable of distinguishing between these cryptic species. Our results shed additional light on how life history variation may affect past and ongoing speciation processes in the speciose Sebastes rockfishes.

Methods

Sample selection, RADseq library preparation, data filtering, and genotype calling

We utilized samples collected from the NWFSC’s Southern California Shelf Rockfish Hook and Line Survey and West Coast Groundfish Bottom Trawl Survey, as well as opportunistically collected samples from the Southwest Fisheries Science Center’s tissue archive that were previously identified as vermilion rockfish, sunset rockfish, or as introgressed sunset rockfish based on nuclear and/or mtDNA data (Hyde et al. 2008a, b; Hyde and Vetter 2009; Unpublished data; Table S1). These samples also included “indeterminates” that could not be confidently assigned to one of the aforementioned groups due to inconclusive molecular results but were suspected sunset rockfish, suspected vermilion rockfish, or suspected sunset/vermilion hybrids. Additionally, we assessed samples previously identified only by vermilion rockfish-type mitochondrial haplotypes, which could either be introgressed sunset rockfish or vermilion rockfish due to a past, one-way introgression event (Hyde et al. 2008a, b). Collectively, 388 samples were selected for extraction that were collected from Central Baja, Mexico (~ 29.2° latitude) to the Strait of Juan de Fuca, Washington (~ 48.4° latitude) and from depths ranging from 44 to 314 m. Genomic DNA was extracted from fin clips stored in 95% ethanol or dried fin punches using the Qiagen DNAeasy Blood & Tissue 96 extraction kit (Qiagen, Inc., Valencia, CA) and then quantified using a BioTek FLX800 Microplate Fluorescence Reader. Samples with extractions yielding < 12.5 ng/μl were not sequenced. Individual DNA concentrations for the 288 remaining samples were normalized to 12.5 ng/μl and 125 ng of starting material was used for RADseq library prep. We used the RADseq protocol of Ali et al. (2016) to construct libraries with the following specific details: (1) genomic DNA from each sample was digested with the restriction enzyme SbfI; (2) libraries were sheared to 300–500 bp using a Qsonica sonicator (Newton, CT); and, (3) 100 bp paired-end sequencing was conducted in 3 lanes using a HiSeq 4000 (Illumina, Inc., San Diego, CA) at the University of Oregon Genomics and Cell Characterization Core Facility (GC3F, Eugene, Oregon).

We ran Stacks v2.4 (Catchen et al. 2013) components individually in a de novo analysis to discover and identify SNPs from raw sequence data. Raw sequence data were quality filtered, demultiplexed, trimmed to 85 bp, and filtered for PCR clones using the programs process_radtags and clone_filter. We followed methods outlined in Paris et al. (2017) and Rochette and Catchen (2017) to determine the optimal Stacks parameters specific to our dataset to minimize erroneous splitting or lumping of loci while simultaneously yielding a high number of polymorphic loci. Specifically, 12 high coverage individuals that included vermilion, sunset and introgressed sunset rockfish were used for optimization, and three was found to be the optimal value both for the maximum number of bp differences between alleles in a sample (− M) and the maximum number of bp mismatches between sample loci (− n). After initial trimming and filtering for low quality reads and PCR clones, loci were identified in each individual using ustacks with a minimum allele depth (− m) of three. A catalog of consensus loci was constructed with cstacks using 46 high coverage individuals (i.e., ustacks mean coverage ≥ 35x) that included vermilion rockfish, sunset rockfish, and introgressed sunset rockfish. Stacks was then used to match individual sample loci to the catalog. Prior to running populations, we excluded individuals with ustacks mean coverage < 15× resulting in 242 individuals. Within populations, loci were dropped that failed to meet the following criteria: present in ≥ 80% of individuals, minor allele count (MAC) ≥ 3, and maximum observed heterozygosity of 70%. We exported the resulting SNP dataset from Stacks and further filtered using VCFtools v0.1.13 (Danecek et al. 2011). We then dropped all but the first SNP from each RADseq locus (− thin 5000), removed loci in individuals that were below 10 × depth of coverage (− minDP 10), refiltered for loci found in ≥ 80% of individuals (− max-missing 0.8), removed individuals with > 30% missing loci the (− remove), and again filtered for loci with a MAC ≥ 3 (− mac 3) in the final dataset, which was exported for downstream analyses. Finally, we checked for duplicate individuals using a custom R script written by Garrett McKinney (https://github.com/gjmckinney/IDduplicateSamples).

Interspecific analyses

Differentiation analyses

To assess interspecific differentiation and identify hybridization/introgression, we ran genetic structure analyses without any a priori assumptions about species boundaries. We used a model-based Bayesian clustering analysis implemented in STRUCTURE v2.3.4 (Pritchard et al. 2000) where ten replicates were run for each number of genetic clusters tested (K = 1–10), each with a burn-in of 10,000 iterations and 100,000 MCMC replicates with admixture allowed (NOADMIX = 0) and no prior location information (LOCPRIOR = 0). Next, we assessed the most likely number of clusters (K) across replicate runs using the Evanno method (Evanno et al. 2005), which assesses the rate of change in log probability of the data between successive values of K (ΔK) as implemented in Structure Harvester (Earl and vonHoldt 2012). We also used the mean likelihood of the model [L(K)] for evaluating K as the Evanno method cannot detect a scenario of K = 1. Finally, CLUMPP v1.1.2 (Jakobsson and Rosenberg 2007) was used to summarize results across replicate STRUCTURE runs and final plots were created in R with ggplot (Wickham 2016). Average individual membership coefficients (Q values) to each cluster were taken from CLUMPP for species assignment. Additionally, we also used principal component analysis (PCA) to summarize the diversity and variation across RADseq loci using the R package adegenet v2.1.3 (Jombart 2008; Jombart et al. 2010). In the PCA, missing data were replaced by the mean allele frequencies with option NA.method = “mean” in the scaleGen() function and the number of principal components kept was determined using a scree plot of eigenvalues.

Identification of diagnostic loci

In order to confidently identify RADseq loci that unambiguously distinguish between sunset and vermilion rockfishes, apparent hybrids (i.e., Q values > 20% admixture) and clear misidentifications (based on PCA and STRUCTURE results) were removed and the dataset was refiltered in VCFtools for a MAC ≥ 3 prior to downstream analyses. Specifically, PC 1 from the interspecific PCA was used to classify individuals into species groups. Individuals previously identified as vermilion rockfish, indeterminate vermilion rockfish, or as vermilion or introgressed sunset rockfish with a PC 1 score of < − 40 were placed in the vermilion rockfish group. Individuals previously identified as sunset rockfish, introgressed sunset rockfish, indeterminate sunset rockfish or as vermilion or introgressed sunset rockfish with a PC 1 score > 40 were placed in the sunset rockfish group (see Table S1 for sample nomenclature specifics). We then calculated locus-specific F-statistics and overall FST (Weir and Cockerham 1984) between sunset rockfish and vermilion rockfish groups using the R package HIERFSTAT (Goudet 2005). Highly differentiated (i.e., FST ≥ 0.99) loci with SNPs > 20 bp from the 5′ end (to allow for primer design in later amplicon panels) were considered as candidates for diagnostic panel development. We report the number of high-quality candidate markers here but the details of primer design and panel testing will be discussed in subsequent work.

Intraspecific analyses

The interspecific dataset—post removal of the admixed and misidentified individuals, and subsequent MAC ≥ 3 filter—was used as a starting point to evaluate population structure and diversity within species. Individuals from this dataset were split into sunset rockfish-specific and vermilion rockfish-specific datasets based on the criteria described above in the “Identification of diagnostic loci” section. Within the respective vermillion and sunset rockfish datasets, non-polymorphic loci were dropped using a MAC ≥ 3 filter in VCFtools. The following analyses were then carried out on each species-specific dataset independently.

Population structure and genetic analyses

PCA and STRUCTURE analyses were run using the methods described above in the “Interspecific differentiation analyses” section. PCA groupings and STRUCTURE assignment results were then used to designate groups for population genetic analyses. Individuals that showed signs of admixture among populations were excluded from downstream analyses as there was no a priori knowledge of population boundaries. Individual sampling sites were mapped using ggplot2 v3.3.3 (Whickham, 2016) and marmap v1.0.4 (Pante and Simon-Bouhet 2013).

Mean and locus-specific genetic diversity indices, global and locus-specific F-statistics, and population pairwise FST (Weir and Cockerham 1984) were calculated using the R package HIERFSTAT v0.04-30 (Goudet 2005). Also in HIERFSTAT, we used boot.vc to calculate 95% confidence intervals (CI) for the global F-statistics, and used boot.ppfst to test and infer significance for pairwise FST comparisons when confidence intervals of 1000 bootstrap replicates did not overlap with zero. We computed the expected heterozygosity within populations (HS) and tested if differences between populations were significant using 999 permutations with HS and Hs.test, respectively, in adegenet v2.1.3 (Jombart 2008).

Results

Stacks filtering parameters on the entire interspecific dataset resulted in 35,910 total RADseq markers, including monomorphic loci. Using VCFtools, we retained only the first SNP on each polymorphic marker, which reduced the data set to 24,026 SNPs. Subsequent filtering for ≥ 10 × coverage for each locus of every individual and for loci found in ≥ 80% of individuals dropped the number of loci to 10,159. After individuals with more than 30% missing loci were removed, 231 individuals remained. Finally, after filtering for a MAC ≥ 3, 10,003 polymorphic loci were retained for interspecific analyses. No duplicate individuals were detected (see Tables S2, S3, & S4 for filtering steps and SNP dataset depth qualities).

Interspecific results

The scree plot for the interspecific PCA based on 10,003 polymorphic loci suggested the first three eigenvalues best explained the variation (Fig. S1). Principal component 1, which explained 20.64% of the variation, showed a clear distinction between sunset and vermilion rockfishes. Although sunset rockfish showed little intraspecific variation in this interspecific PCA, vermilion rockfish showed clear intraspecific variation with three distinct groupings separated by both PC2, which explained 3.37% of the variation (Fig. 1), and PC3, which explained 1.11% of the variation (Fig. S2). Replicate STRUCTURE runs overwhelmingly suggested K = 2 as the most likely scenario based on (ΔK), while [L(K)] strongly suggested that K = 1 is the least likely scenario (Table S5). The K = 2 plot shows that sunset and vermilion rockfishes are clearly differentiated from each other in most cases (Fig. 2). However, both PCA and STRUCTURE analyses revealed 12 individuals that were either mislabeled or misidentified (e.g., 3 of the 4 indeterminate hybrids clearly grouped with either vermilion or sunset rockfish), as well as three individuals that represent likely hybrids (i.e., one indeterminate hybrid and two individuals previously identified as vermilion; Figs. 1, 2). These individuals were not included in either sunset or vermilion groups for detecting diagnostic SNPs and were also excluded from subsequent intraspecific analyses. Although individuals with low levels of admixture were rare, the direction of nuclear introgression was more commonly from sunset rockfish into vermilion rockfish (Fig. 2).

Fig. 1
figure 1

Principal component analysis on 10,003 interspecific SNPs from 231 samples previously identified (based on nuclear and/or mtDNA data) as vermilion rockfish, indeterminate vermilion rockfish, vermilion or introgressed sunset rockfish, indeterminate hybrid, sunset rockfish, introgressed sunset rockfish, and indeterminate sunset rockfish

Full size image
Fig. 2
figure 2

Bayesian clustering results from STRUCTURE (K = 2) for all 231 samples and 10,003 loci. Individuals are arranged based on previous molecular identification (nuclear and/or mtDNA data) from vermilion rockfish, indeterminate vermilion rockfish, vermilion or introgressed sunset rockfish, indeterminate hybrid, sunset rockfish, introgressed sunset rockfish, to indeterminate sunset rockfish

Full size image

Interspecific F-statistics and diagnostic loci

Sunset and vermilion rockfishes showed strong overall differentiation with an FST = 0.5478 (0.5382–0.5575 95% CI). We detected 203 highly differentiated (i.e., FST ≥ 0.99) SNPs located > 20 bp from the 5′ end of their respective RADseq paired-end contig sequence, 117 of which were fixed (FST = 1). Further testing is required but these loci appear to be promising candidates for developing a diagnostic SNP assay to distinguish sunset and vermilion rockfish.

Sunset rockfish

One hundred fourteen individuals were assigned to sunset rockfish. Removal of non-polymorphic loci using a MAC ≥ 3 filter resulted in 5043 loci.

Population genetic results

The scree plot for the PCA based on the filtered sunset rockfish dataset suggested the first two eigenvalues best explained the variation (Fig. S3). The PCA plot showed two groups separated along PC1 (1.76% of the variation), which we refer to as S-A (PC1 > 5) and S-B (PC1 < 0; Fig. 3). Replicate STRUCTURE runs corroborate the PCA pattern, with both ΔK and [L(K)] finding K = 2 as the most likely scenario (Table S6). Most individuals clearly assigned to one of the two clusters (Q > 0.7), however, one individual showed signs of significant admixture with Q assignment values of ~ 0.52/0.48 to Group S-A and S-B, respectively (Fig. 4). This individual was excluded from population genetic analyses resulting in 113 individuals total (48 and 65 individuals in S-A and S-B, respectively).

Fig. 3
figure 3

Principal component analysis for sunset rockfish based on 5,043 intraspecific SNPs reveal two distinct groupings, S-A and S-B. Individual sampling latitude is shown by color and mitochondrial haplotypes are depicted by shape

Full size image
Fig. 4
figure 4

Bayesian clustering results from STRUCTURE (K = 2) for sunset rockfish based on 5,043 intraspecific SNPs. Individuals are arranged based on distinct PCA groupings, S-A and SB

Full size image

Individuals with introgressed and non-introgressed mitochondrial haplotypes—vermilion and sunset rockfish mtDNA, respectively—were found in near equal numbers in both S-A (25:23) and S-B (34:31; Fig. 3). Individuals in S-A and S-B showed a difference in depth of capture (Fig. S4A) with S-B being more common in intermediate depths while S-A was more common at shallower and deeper depths. When each sunset group was separated by mitochondrial haplotype, introgressed haplotypes were generally more common in more shallow habitat while non-introgressed haplotypes were relatively common throughout sampling depths (Fig. S4B). Notably, S-A and S-B groups appeared to show differences in latitudinal distribution in the PCA (Fig. 3). More specifically, when individual sunset rockfish geographic coordinates were mapped, S-A individuals were found in higher latitudes in close proximity to the coast while S-B individuals dominated the lower latitudes on the offshore banks and islands of the Southern California Shelf (Fig. 5).

Fig. 5
figure 5

Individual sunset rockfish sampling locations withs colors corresponding to associated intraspecific grouping, S-A and S-B

Full size image

Population genetic diversity, as estimated by expected heterozygosity (HS), was similar but slightly lower in S-A compared to S-B (0.2123 and 0.2138, respectively) but not significantly different from each other (p = 0.185). Global F-statistics and mean genetic diversities for sunset rockfish are reported in Table 1. Pairwise FST comparison (between S-A and S-B) yielded a significant value of 0.0122, which was identical to the reported global FST, as expected. Locus-specific estimates of observed heterozygosity, within population gene diversity, overall gene diversity, gene diversity among samples, fixation index, inbreeding coefficient, and Jost’s measure of population differentiation (Jost’s D) are reported in Table S7. Sunset rockfish locus-specific corrected FST (FSTP) values ranged from − 0.0215 to 0.5564.

Table 1 Intraspecific global F-statistics (95% CI) and diversity indices for sunset rockfish and vermilion rockfish based on 5403 and 6406 SNPs, respectively
Full size table

Vermilion rockfish

One hundred five individuals were assigned as vermilion rockfish. Removal of non-polymorphic loci using a MAC ≥ 3 filter resulted in 6406 loci.

Population genetic results

The scree plot for the PCA based on the filtered vermilion rockfish dataset suggested the first two eigenvalues best explained the variation (Fig. S5). The PCA plot showed three distinct groups separated along both PC1 (8.05% of the variation) and PC2 (2.38% of the variation), which we refer to as V-A (PC1 > 45), V-B (PC2 > 20), and V-C (PC1 < − 10; Fig. 6). Replicate STRUCTURE runs show similar patterns of differentiation with ΔK suggesting that K = 2 was the most likely scenario followed by K = 3, while [L(K)] suggested K = 3 was slightly more likely than K = 2 (Table S8). Plots for both K = 2 and K = 3 corroborate the PCA pattern with three distinct cluster assignment groups (Fig. 7). Both PCA and STRUCTURE identify an individual that showed significant admixture between groups V-B and V-C (Figs. 6,7). This individual was excluded from population genetic analyses resulting in 104 individuals total (12, 17, and 75 individuals in V-A, V-B, and V-C, respectively).

Fig. 6
figure 6

Principal component analysis for vermilion rockfish based on 6406 intraspecific SNPs reveal three distinct groupings, V-A, V-B, and V-C. Individual sampling latitude is shown by color

Full size image
Fig. 7
figure 7

Bayesian clustering results from STRUCTURE for vermilion rockfish based on 6,406 intraspecific SNPs for K = 2 and K = 3. Individuals are arranged based on distinct PCA groupings, V-A, V-B, and V-C

Full size image

Depth of capture was not plotted for vermilion rockfish groups as these data were lacking for all V-A individuals and only available for 3 out 17 V-B individuals. Geographic coordinates, however, were available for all samples and, like sunset rockfish groups, vermilion rockfish groups showed a difference in latitudinal distribution in the PCA (Fig. 6). Specifically, when vermilion rockfish individuals were mapped, V-A individuals were restricted to higher latitudes north of Point Reyes, CA, V-B individuals were common in Central California from Cape Mendocino south to Point Conception, and V-C individuals dominated south of Point Conception but were relatively rare to the north (Fig. 8). Notably, V-B was absent from the Southern California Shelf but was sampled at the Mercado Negro fish market in Ensenada, Baja California, Mexico. However, this is the largest fish market in Baja and these Group B individuals may have been caught at a distant location.

Fig. 8
figure 8

Individual vermilion rockfish sampling locations withs colors corresponding to associated intraspecific PCA grouping, V-A, V-B, and V-C (Note: the southernmost V-B was opportunistically sampled from the Mercado Negro fish market in Ensenada, Baja California, Mexico)

Full size image

Population genetic diversities, as estimated by expected heterozygosity (HS), for V-A, V-B, and V-C were 0.1618, 0.1939, and 0.2026, respectively, and were significantly different from each other in all three comparisons (p = 0.001 for V-A/V-B and V-A/V-C, p = 0.024 for V-B/V-C). Notably, HS was highest in V-C, the southernmost group, and progressively decreased to the north in V-B and V-A. Mean genetic diversities and global F-statistics with 95% CI for vermilion rockfish are reported in Table 1. Pairwise FST comparisons between V-A and V-B, V-A and V-C, and V-B and V-C were 0.1099, 0.2326, and 0.1147, respectively; all comparisons were significant. Locus-specific estimates of observed heterozygosity, within population gene diversity, overall gene diversity, gene diversity among samples, fixation index, inbreeding coefficient, and Jost’s D are reported in Table S9. Vermilion rockfish locus-specific corrected FST (FSTP) values ranged from − 0.0695 to 0.8977. Locus-specific metrics were estimated for each pairwise comparison between nominal vermilion rockfish groups ad hoc in order to potentially detect loci with fixed differences (i.e., FST = 1). No fixed differences were detected between V-B and V-C, however five were detected between V-A and V-C, and one was detected between V-A and V-B (data not shown).

Discussion

Here we identified a set of diagnostic markers for species identification in the commercially and recreationally important sunset and vermilion rockfishes. Additionally, many of the RADseq markers generated for characterizing interspecific differentiation also provide powerful insight into the intraspecific population structure of both species and reveal dissimilar levels of differentiation. Specifically, vermilion rockfish exhibit much stronger signals of differentiation compared to the relatively weak but significant differentiation observed in sunset. Notably, introgressed sunset rockfish clearly group with conspecifics (i.e., non-introgressed sunset rockfish) and show little evidence of admixture with vermilion, a clear delineation not possible using solely mitochondrial haplotype data. Within both species, the distinct genetic groups exhibit dissimilar geographic distributions that appear correlated to latitude and/or bathymetry. Notably, both species have a distinct genetic group on the Southern California Shelf.

Interspecific differentiation patterns

Genome-wide RADseq markers clearly distinguish between sunset and vermilion rockfishes in all interspecific analyses. Wright’s fixation index between the sister species is similar to that reported in other closely related Sebastes spp. (Benestan et al. 2020). However, we detected evidence of nuclear introgression, which includes an apparent F1 hybrid (that was previously identified as an indeterminate hybrid), several possible back crosses, as well as individuals with lower levels of introgression (Fig. 2). Interestingly, the direction of introgression appeared disproportionate with more evidence of sunset introgression into vermilion. Previous genetic work in the sister species using microsatellites, however, suggested low levels of bidirectional nuclear admixture (Hyde et al. 2008a, b).

Differences in intraspecific population differentiation

Sunset rockfish exhibit weak but significant population genetic differentiation (Table 1) between two groups with dissimilar geographic distributions which is a novel finding. Specifically, sunset individuals in group S-A are restricted to near shore habitat in the northern portion of the range down to approximately Dana Point and are conspicuously missing from the offshore banks and islands of the Southern California Shelf, except for Catalina Island (Fig. 5). Interestingly, a similar geographic distribution was observed for highly differentiated populations of cowcod (Sebastes levis), although the mechanisms driving the differentiation and distribution were unknown (Hess et al. 2014). The near equal ratio of introgressed to non-introgressed sunset rockfish in both groups demonstrates that mitochondrial haplotypes are not correlated with this pattern of differentiation. Depth, however, may contribute to the observed intraspecific differentiation as the two groups exhibit dissimilar depth distributions (Fig. S1A). Indeed, segregation by depth has played a significant role in the differentiation of closely related Sebastes spp. and the diversity of rockfishes as a whole (Heras and Aguilar 2019; Hyde et al. 2008a, b; Hyde and Vetter 2007; Ingram 2011; Love et al. 2002). However, the possible relationship between depth and the distribution of S-A and S-B genetic groups observed here could be a sampling artifact. Planned work to disentangle the thousands of nominal vermilion rockfish samples, and associated metadata, may yield additional sunset rockfish data that would allow for a more robust investigation into this potential relationship. Undoubtedly, subsequent population genomic work on sunset rockfish should include samples where S-A and S-B overlap (i.e., Point Conception), allowing for more detailed investigations into potential mechanisms behind the apparent differences in depth distribution. Similar distribution patterns have been attributed to factors such as ecological competition and life history differences in other Sebastes spp. (Hyde et al. 2008a, b; Ingram 2011; Larson 1980).

Vermilion rockfish display an order of magnitude stronger signal of population genetic differentiation and are instead characterized by three distinct groups with latitudinally stratified distributions (Table 1, Fig. 8). Inference into whether depth distribution differs amongst the vermilion genetic groups is not possible here as we lack depth data for V-A and V-B individuals. However, as discussed in the previous paragraph, future diagnostic work on the thousands of nominal vermilion rockfish individuals may allow for insight into this question. Group V-C dominates south of Point Conception, which was previously identified as a strong genetic break for this species (Budrick 2016; Hyde and Vetter 2009) and many other taxa, including other Sebastes spp. (Briggs 1974; Buonaccorsi et al. 2004, 2005; Pelc et al. 2009). The observed pattern of progressively decreasing HS moving from the southernmost group (V-C) poleward to the northernmost group (V-A) suggests that V-C, with the highest observed diversity, may represent the ancestral population. It follows that V-B would have differentiated from V-C in an expansion north, with a third poleward expansion of V-A subsequently differentiating from V-B. This pattern, along with the pairwise FST comparisons, corroborates a stepping-stone pattern of gene flow for vermilion rockfish, as previously suggested (Budrick 2016; Hyde and Vetter 2009). However, gene flow amongst these distinct genetic groups appears to be fairly low in general based on the strong pairwise FST comparisons (including fixed differences between V-A/V-C and V-A/V-B) and significant differences in HS between all group comparisons. Indeed, the degree of divergence between vermilion rockfish groups is similar to levels of interspecific differentiation observed between sister Sebastes species, such as the gopher rockfish and black and yellow rockfish (Narum et al. 2004) as well as the blue rockfish and deacon rockfish (Burford and Bernardi 2008). The central coast of California harbors all three distinct vermilion genetic groups, corroborating the observation from Hyde and Vetter (2009) that the highest vermilion diversity exists in this area. Given that adults of all three genetic groups overlap in this region, the potential for inter-group breeding appears highly plausible. Although limited larval dispersal due to parturition timing in conjunction with the seasonal prevailing oceanographic conditions is thought to restrict gene flow amongst regions in vermilion rockfish (Hyde and Vetter 2009), this scenario would not produce a barrier for these apparently sympatric genetic groups. Notably we did detect an apparent hybrid between V-B and V-C groups that suggests interbreeding between groups does occur. However, based on the otherwise strong differentiation patterns, it appears that a barrier to gene flow exists. One potential explanation could be assortative mating, a potentially strong barrier to gene flow. Indeed, Sebastes rockfishes are viviparous with internal fertilization and courtship displays have been characterized in other Sebastes spp. (Helvey 1982; Shinomiya and Ezaki 1991), suggesting mate selection. Further work is required to better understand this strong intra-specific differentiation in vermilion rockfish.

Distinct life history behaviors may contribute to the dissimilar levels of intraspecific differentiation observed in sunset and vermilion rockfish. The weaker intraspecific differentiation characterizing sunset rockfish could be explained by increased connectivity through ontogenetic migration of juveniles to offshore banks and deeper habitat suitable for adults. Canary rockfish, S. pinniger—the sister taxon to the sunset and vermilion rockfish clade—share a similar ontogenetic migration of juveniles to adult habitat (Vetter and Lynn 1997), which is common in Sebastes spp. (Love et al. 1991), and do not exhibit any significant population structure (Andrews et al. 2018; Budrick 2016). Whereas vermilion rockfish, which show strong intraspecific differentiation, may have restricted gene flow due to the lack of this ontogenetic shift, in conjunction with strong biogeographic barriers in its range as suggested by Hyde and Vetter (2009). Specifically, because there is little movement as juveniles or adults, the larval stage is the most likely method of dispersal, which can be greatly affected by prevailing currents, upwelling, and offshore jet (Leis 1991). The absence of an ontogenetic migration is also implicated for precipitating the divergence between sunset and vermilion roughly 2.3 Ma (Hyde et al. 2008a, b; Hyde and Vetter 2007). Interestingly, we find that both sunset and vermilion rockfish show a pattern of intraspecific variation where a single genetic group dominates the offshore banks and islands of the Southern California Shelf. The eddies around the islands and the gyre of prevailing circulation that characterize the Southern California Bight (Bray et al. 1999) may affect larval dispersal patterns but additional samples and work is needed to properly assess this possibility.

Because the underlying purpose of generating these RADseq data was to identify and characterize candidate markers for a diagnostic SNP panel, our intraspecific analyses are limited by relatively small sample sizes. A robust sunset population genetic study that includes more samples, especially from the southern portion of the range is warranted based on the observed structure with this limited dataset. The population genetic structure of vermilion has already been thoroughly evaluated (Budrick 2016; Hyde and Vetter 2009), although further insight into the distribution of the three observed clusters in this study at the periphery of the range could be useful for management purposes. Additionally, focused sampling of vermilion rockfish in Central California from Cape Mendocino to Point Conception, where all three groups co-occur, could shed light on the degree of overlap in depth distributions. This would allow clearer insight into the potential mechanisms restricting gene flow in this region of overlap and further our understanding of the processes driving differentiation in the highly speciose Sebastes rockfishes.

Fisheries management implications

Although the cryptic sunset rockfish was identified over a decade ago (Hyde et al. 2008a, b), this distinct species is still managed collectively with vermilion rockfish as a single species. Managing multiple, discrete stocks as one ostensibly homogeneous population may result in an increased risk of overfishing to one or more of the stocks or the inadvertent (and potentially undiagnosed) loss of genetic diversity if the spatial distribution of fishing effort is misaligned with the population’s underlying stock structure (Hauser et al. 2002; Kerr et al. 2017; Laikre et al. 2005; Okamoto et al. 2019; Spies and Punt 2015; Sterner 2007). Given the current management regime, these and other potential risks to this species pair have not been evaluated which is a significant concern for such an economically and ecologically important assemblage.

This research also furthers our knowledge of the degree of intraspecific genetic structure characterizing each species, which should also be accounted for in an effective management strategy (Reiss et al. 2009; Spies and Punt 2015; Waples et al. 2008). Our findings suggest weak but significant intraspecific differentiation for sunset rockfish, however our sampling was relatively limited and a thorough population genetic study is warranted to better characterize this structure. Although more data are required, our results provide no clear indication of distinct management needs for these genetic groups. Vermilion rockfish, however, show relatively strong signals of differentiation that appear to be correlated with latitude. In California, the vermilion rockfish species complex has historically been managed with boundaries at Cape Mendocino and/or Point Conception, which generally match our findings of where the genetic breaks occur. Specifically, V-A is the only group found north of Cape Mendocino, V-B predominantly occurs between Cape Mendocino and Point Conception, while V-C dominates to the south of Point Conception. Therefore, these management boundaries which have been used for the species complex may be ideal boundaries for vermilion rockfish as suggested by Hyde and Vetter (2009).

This effort’s next step is to use the diagnostic SNPs described here to delineate to species ~ 25,000 nominal vermilion rockfish samples collected during previous NWFSC Southern California Shelf Rockfish Hook and Line and West Coast Groundfish Bottom Trawl Surveys. Separating these existing samples will support finer resolution of the differences in spatial and depth distribution between the two species as well as comparisons of their respective age and growth, reproductive biology, and diet and trophic patterns. It will also provide much needed insight into both interspecific differentiation as well as intraspecific population genetic structure. These analyses will help inform decisions about whether the two species are most appropriately assessed and managed collectively or individually and ensure the PFMC has the best information possible to reduce the risk of overfishing to both components of this valuable fishery.