Quantifying connectivity between endangered or threatened marine populations is critical information for scientifically sound management. Of the 86 species managed by the Pacific Fishery Management Council (PFMC) on the West Coast of the United States, over 75% are rockfishes, and of those species, 27 were once deemed ‘at risk’. Although most stocks have been successfully rebuilt along the Washington Coast, Puget Sound stocks have yet to recover despite fisheries closures. The connectivity of Puget Sound stocks to coastal populations is relatively unknown, despite the potential of recruitment subsidies from the coast and considerable interest in reopening many fisheries for recreational use. The importance of accurate connectivity estimates was demonstrated by recent research on two of three Endangered Species Act (ESA) listed species in Puget Sound, which showed that one species was not sufficiently distinct to warrant listing as a separate Puget Sound distinct population segment (DPS). The common assumption of isolation of Puget Sound populations, which stems primarily from research on three hybridizing species, may therefore be erroneous. This study aimed to provide an analysis of the connectivity of five Puget Sound rockfish populations and identify distinct population segments where appropriate. Samples from five species (Black (Sebastes melanops), Yellowtail (S. flavidus), Redstripe (S. proriger), Greenstriped (S. elongatus), and Puget Sound (S. emphaeus)) were collected in three areas within and one area outside of Puget Sound and analyzed at over 12,000 restriction-site associated DNA sequencing (RADseq) loci. We found unique species-specific patterns of genetic diversity, attributable to multiple extrinsic and intrinsic factors. In particular, Black and Puget Sound Rockfish showed no genetic differentiation; Yellowtail and Greenstriped Rockfish were structured according to known geographic barriers; and Redstripe Rockfish revealed evidence for temporal genetic differentiation, suggesting that irregular recruitment influences population structure. None of the species followed the DPS boundaries generally assumed for rockfish, further emphasizing the importance of species-specific management for the effective recovery of these rockfish populations.

Sampling Procedure 

We used 279 samples from five species of rockfish (Black (S. melanops), Yellowtail (S. flavidus), Redstripe (S. proriger), Greenstriped (S. elongatus), and Puget Sound (S. emphaeus), Figure 1) that were collected in 1999-2021 in previous surveys of WDFW, NOAA-NMFS, and the Department of Fisheries and Oceans (DFO Canada). Individual fin clips were preserved in 95% ethanol or dried on filter paper. Samples were collected from multiple locations (see Figure 3), grouped into four different regions: 1) southern Puget Sound (Puget Sound proper, south of Admiralty Inlet, SPS), 2) British Columbia (Canadian Salish Sea north of the US/Canada border, BC), 3) northern Puget Sound (US Salish Sea north of Admiralty Inlet, NPS and east of the Victoria Sill) and 4) the US west coast (US Pacific Coast west of Victoria Sill, WC). Due to differences in the abundance and distribution of species across this geographic range, we have no Puget Sound Rockfish from WC and one Greenstriped Rockfish from NPS. 

DNA Extraction, Library Preparation, and Sequencing 

Genomic DNA was extracted using the Nexttec DNA isolation kit (Nexttec Incorporated, Middlebury, VT, USA) following the manufacturer’s protocol and quantified using a Qubit Fluorometer (ThermoFisher Scientific, Waltham, MA, USA). DNA concentration was normalized to 125ng in 10 µL of molecular-grade water. Restriction site-associated DNA sequencing (RADseq) libraries were prepared using a version of the Ali et al. (2016) protocol without the targeted bait capture step, referred to in the literature as BestRAD (https://github.com/merlab-uw/Protocols/blob/main/bestRAD). Briefly, genomic DNA was digested using the SbfI enzyme. An adapter (P1) containing a forward amplification primer site, an Illumina sequencing primer site, and an individual 6 bp barcode was ligated to each fragment at the restriction site end. Fragments were then randomly sheared using sonication and size-selected to 300-500 bp in length. Subsequently, P2 adapters were ligated to the reverse end and libraries were amplified by PCR. Each library was assessed for quality on a 1% agarose gel and a Bioanalyzer DNA 1000 kit (Agilent Technologies, Santa Clara, CA). Libraries were pooled in equimolar amounts and sequenced on a NovaSeq (paired-end, 116 bp or 150 bp) at the University of Oregon, either an S4 or SP run type. Ninety-six individuals were randomly included in one of six RADseq libraries to avoid any lane effect (Leigh et al., 2018).

Initial Filtering 

Raw sequence data were quality-checked using multiQC (Ewels et al., 2016). Prior to SNP calling and genome alignment, raw sequences were demultiplexed using process_radtags in the STACKS v2.60 pipeline (Catchen et al., 2011; Rochette et al., 2019). Sequences were trimmed to 104 bases and filtered for quality. Individuals with fewer than 250,000 total reads were excluded from downstream analysis (Krohn et al., 2018). Our paired-end sequences were then aligned to the Honeycomb Rockfish (Sebastes umbrosus) genome from GenBank (NCBI Accession Number: PRJNA562243) with Bowtie2 v2.4. using the ‘very-sensitive’ option (Langmead & Salzberg, 2012). The Honeycomb Rockfish genome is one of only two annotated full genomes and was chosen due to its closer phylogenetic relationship to our species (Hyde & Vetter, 2007a). Following genome alignment, SNP calling and basic population genetics statistics were calculated using the gstacks and populations modules from the STACKS pipeline. SNPs were called if they had a minimum mapping quality of 40. SNPs were filtered following published recommendations (O’Leary et al., 2018) requiring that loci meet the following criteria: minimum genotype depth ≥ 5, mean minimum read depth ≥ 15, genotype call rate ≥ 80%. We subsampled SNPs by choosing one SNP on each RADtag with the highest minor allele frequency. SNPs with genotype frequencies that were significantly different than expectations under Hardy-Weinberg Equilibrium (HWE) were removed using the following procedure: p-values were calculated across samples for each population using the r package pegas v1.1 (Paradis, 2010). P-values were then combined across samples for each locus using Fisher’s combination of probabilities, and adjusted to q-values for the false discovery rate (Benjamini & Hochberg, 1995). Loci with q-values below 0.05 were considered significantly out of HWE and removed from downstream analysis. 

Misidentification Analysis

To identify any cases of 1) misidentification of species during field sampling or 2) interspecific hybridization in our dataset, raw sequences for all individuals were analyzed together for eight species (five from this study, and Brown, Quillback, and Copper Rockfish from Wray et al (in prep)) immediately after genome alignment. SNP calling, and basic population genetics statistics were calculated using the gstacks and populations modules from the STACKS pipeline. SNP filtering followed the same protocol as in the species-specific analyses (O’Leary et al., 2018). Unlike Different from the species-specific analysis, however, we intended to we avoided SNPs with fixed differences between species, because this would likely find differences between only two species. Therefore, we chose the first SNP on each RADtag using the –write-single-snp option in populations. Additionally, we did not filter for HWE because a reduction of heterozygosity due to species-specific (subpopulation) structure would likely influence HWE p-values due to the Wahlund effect (REF). We plotted all eight species together in a principal components analysis (PCA). Any individuals that visually grouped with a species different than its field identification were removed from downstream analysis.