Speciation genomics and the role of depth in the divergence of rockfishes (Sebastes) revealed through Pool‐seq analysis of enriched sequences

Abstract Speciation in the marine environment is challenged by the wide geographic distribution of many taxa and potential for high rates of gene flow through larval dispersal mechanisms. Depth has recently been proposed as a potential driver of ecological divergence in fishes, and yet it is unclear how adaptation along these gradients' shapes genomic divergence. The genus Sebastes contains numerous species pairs that are depth‐segregated and can provide a better understanding of the mode and tempo of genomic diversification. Here, we present exome data on two species pairs of rockfishes that are depth‐segregated and have different degrees of divergence: S. chlorostictus–S. rosenblatti and S. crocotulus–S. miniatus. We were able to reliably identify “islands of divergence” in the species pair with more recent divergence (S. chlorostictus–S. rosenblatti) and discovered a number of genes associated with neurosensory function, suggesting a role for this pathway in the early speciation process. We also reconstructed demographic histories of divergence and found the best supported model was isolation followed by asymmetric secondary contact for both species pairs. These results suggest past ecological/geographic isolation followed by asymmetric secondary contact of deep to shallow species. Our results provide another example of using rockfish as a model for studying speciation and support the role of depth as an important mechanism for diversification in the marine environment.


| INTRODUC TI ON
Understanding mechanisms of speciation in the marine environment remains difficult due to the lack of apparent geographical barriers and high rates of gene flow among populations (Dennenmoser et al., 2017). Most marine species demonstrate high dispersal capabilities and connectivity among populations, which can impede local adaptive processes and differentiation (Bierne et al., 2003;Carreras et al., 2017). However, even with the homogenizing effects of gene flow, there is potential for local adaptation that may be the driving force for genomic differentiation in the marine environment (Whitney et al., 2018).
Early work on speciation in marine fishes was thought to be a consequence of geographic isolating mechanisms. The formation of land barriers (Bermingham et al., 1997;Bernardi et al., 2004), islands (Leray et al., 2010), and physical boundaries generated from oceanographic processes (Gaither & Rocha, 2013;Hubert et al., 2012) were used from a biogeographical perspective to describe speciation patterns in marine fishes. However, the role of pelagic larval duration in contributing to gene flow among populations suggested that allopatric divergence may be rarer in marine fish (reviewed in Bindea et al., 2013). Although large-scale allopatric events can drive marine speciation, there is evidence that other isolating mechanisms occur in the marine environment (Faria et al., 2021;Rocha et al., 2005).
Studies have demonstrated the possibility for closely related species to be sympatrically distributed, indicating ecology may be important (Burford, 2009;Crow et al., 2010;Rocha et al., 2005). Although Rocha et al. (2005) found evidence for sympatrically distributed reef fishes, there remains a lack of knowledge for how ecological speciation applies to temperate marine environments. Furthermore, sympatric or parapatric distribution of species contradicts existing knowledge that gene flow can be a barrier to speciation. This creates an apparent "marine speciation paradox," or how can marine speciation occur in the face of high apparent gene flow (Faria et al., 2021;Johannesson, 2009)?
The model of ecological speciation is important in understanding the speciation process in the marine environment (Bernardi, 2013;Puebla, 2009). Numerous studies have now documented the role that ecological specialization, especially in fishes, plays in speciation. A number of environmental factors have been documented that lead to ecological divergence in the marine environment; these include temperature (Teske et al., 2019) and salinity (Momigliano et al., 2017), which are often correlated with other habitat characteristics (e.g., depth and latitude). Depth has already been identified as a potential factor in the diversification of rockfishes (Behrens et al., 2021;Heras & Aguilar, 2019;Hyde et al., 2008;Ingram, 2010;Sivasundar & Palumbi, 2010), and depth may be important in driving speciation for other marine organisms (Carlon & Budd, 2002;Gaither et al., 2018;Hirase et al., 2021;Prada & Hellberg, 2013). Thus, adaptation to these environmental differences can lead to divergence in life history traits, such as spawning behavior, which subsequently drives reproductive isolation between incipient species. Rockfishes (genus Sebastes) inhabit temperate waters across the Atlantic and Pacific Ocean, with 60 different species found in the North Pacific that have radiated over the past 5 million years (Johns & Avise, 1998). Species are found from rocky intertidal habitats to depths greater than 1500 m (Love et al., 2002). Given the ecological partitioning and habitat similarity between recently diverged forms of rockfish, ecological speciation may have contributed to their divergence (Behrens et al., 2021;Burford, 2009;Pavoine et al., 2009).
Depth has been proposed as an important component in the diversification of rockfishes (Behrens et al., 2021;Heras & Aguilar, 2019;Hyde et al., 2008;Ingram, 2010;Sivasundar & Palumbi, 2010 (Hyde & Vetter, 2007). Our goals are to determine whether depth-segregated speciation is a result of selective sweeps that generates islands of genomic differentiation or "divergence islands" (Via, 2012;Wolf & Ellegren, 2017). We will also examine islands of genomic differentiation to see whether they are shared across species pairs. Sharing of these divergence islands may indicate parallel evolutionary pressures in depth adaptation. Finally, we also investigated the demographic history of speciation in these species pairs, to see whether similar patterns exist and whether these patterns are consistent with ecological speciation.

| Sample collection
Ethanol-preserved fin clips for the following depth-separated species pairs were obtained from S chlorostictus-S. rosenblatti and S. crocotulus-S. miniatus (Table 1). High molecular weight DNA was obtained using standard phenol-chloroform methods followed by ethanol precipitation (Sambrook & Russell, 2006). DNA integrity was checked on a 2% agarose gel and quantified on a Qubit fluorometer before preparing the samples for Illumina library preparation.

| Targeted sequence capture design
We designed a series of oligonucleotide capture baits that could efficiently enrich DNA sequencing libraries across Perciformes, a large order of more than 2200 species that includes Sebastes (Daane et al., 2016(Daane et al., , 2019Nelson et al., 2016). We targeted protein-coding exons, conserved non-coding elements (CNEs), miRNAs and ultra-conservative elements (UCNEs) for enrichment.

| Exome sequencing
Exome sequencing of 20 individuals from four different species was done using a pool-seq approach ( Table 1) (Daane et al., 2016;Schlötterer et al., 2014). DNA from 20 individuals was pooled in equimolar amounts, and libraries were constructed using the Kapa HyperPlus kit (Table S1). Enrichment of exome sequences was done following the approach of Daane et al. (2019) pri) with bwa-mem (Li & Durbin, 2010). Resulting BAM files were converted to mpileup format with samtools (Table S1) (Daane et al., 2019), and regions surrounding indels were masked with the identify-indelregions.pl script for subsequent analysis. Allele frequencies were estimated with Popoolation2 (Kofler et al., 2011), and loci with a minor allele frequency (MAF) less than 0.05 were removed for subsequent analyses. We use this MAF cutoff to avoid any bias from fixed or nearly fixed variants between species.
We identified divergence islands across species pairs using outlier approaches. We estimated F ST for species pairs using the parameters in Popoolation2: -suppress-noninformative -min-count 6 -min-coverage 80 -max-coverage 500 -min-covered-fraction 1 -window-size 1 -step-size 1 -pool-size 20 (Kofler et al., 2011;Li et al., 2009). Results from the F ST analysis were then plotted utilizing qqman to generate Manhattan Plots in R (Turner, 2018). Our comparisons include two species pairs that span a range of divergence across the speciation continuum (Behrens et al., 2021;Hyde et al., 2008). Identification of clear regions of divergence in the more diverged species pair possesses additional challenges, as drift and recombination may erode any signals associated with speciationrelated divergence (Quilodrán et al., 2020).

| Sliding window analysis
To identify divergence islands, we applied an approach similar to Holliday et al. (2016) and Renaut et al. (2013). We performed a sliding window analysis of F ST across each chromosome, using a window of 10 adjacent SNPs and sliding the window every two SNPs, to identify the number of regions that contained SNPs greater than the top one percentile for F ST of the genome-wide analysis. To assess significance, we randomly sampled 10 SNPs from across the genome with replacement for their F ST values 100,000 times. For each of these subsamples, we estimated the proportion of top one percentile F ST SNPs present. The proportion of top one percentile SNPs in each window, over the resampled dataset, was used to determine significance for the original dataset (p < .001) to reduce the signal from false positives.

| Gene ontology enrichment analysis
We identified genes found in outlier windows from our sliding window analysis. These served as our candidate list of genes that was then compared with the background list of all annotated genes found from our exome dataset. To test for enrichment of molecular pathways and/or function between the background and candidate lists and to allow for the visualization of the gene ontology networks, we used Cytoscape-CLUEGO and its plug-in CLUEPedia (Bindea et al., 2009(Bindea et al., , 2013. Cytoscape-CLUEGO utilizes hypergeometric testing followed by Bonferroni multiple testing corrections between a candidate gene list and a custom background list (Bindea et al., 2009). For Cytoscape-CLUEGO, all annotations were made using the D. rerio genome from the Gene Ontology Consortium (Ashburner et al., 2000), provided as it was the closest related species to genus Sebastes in this analysis package. Finally, Cytoscape-CLUEGO groups genes by GOterm to avoid redundancy in the results.

| Demography of speciation
We utilized folded site frequency spectra in δaδi (Gutenkunst et al., 2009) to determine the demographic history of speciation in each of the species pairs. To generate datasets for this analysis, we used all identified SNPs from Popoolation2 (before filtering for MAF as above). To assure independence (linkage equilibrium) of each SNP, we randomly sampled one SNP every 1,000,000 bp across the genome. Raw SNP frequencies were converted into δaδi SNP format using genomalicious (https://rdrr.io/githu b/j-a-thia/ genom alici ous/). We explored seven simple two-population models in δaδi: no migration, symmetric migration, asymmetric migration, symmetric migration followed by isolation, asymmetric migration followed by no migration, secondary contact with symmetric migration, and secondary contact with asymmetric migration. We hypothesize that if ecological speciation has occurred in these species pairs, due to adaptation to different depths, we would expect to observe a demographic history with at least some gene flow. We utilized the δaδi optimization procedure from Portik et al. (2017) (https://github.com/dport ik/dadi_pipeline). We ran four iterations of optimizations for each model with 10, 20, 30, and 40 replications, respectively. For each species pair, we compared models using AIC and Δ AIC. We did not convert parameters from the best fit model into biologically meaningful values, as our goal was simply to reconstruct a reasonable demographic scenario for each species pair.

| Genome assembly and coverage
We obtained sequences for pooled sequenced from each species that contained 20 individuals. The average read depth was 10.36 across species and >99% of the reads mapped back to the S. umbrinus genome (Table S1).  (Figure 1; Tables S2 and S3). There were two windows that were shared in both comparisons, one on chromosome 6 and the other on chromosome 12 (Figure 1, Tables S2-S4).

| Demography of speciation
In order to assess the demographic history of speciation within our two species pairs, we tested seven models of population divergence using δaδi. We used pruned datasets that consisted of 5368 and 5310 SNPS for the S. rosenbaltti-S. chlorostictus and S. miniatus-S. crocotulus species pairs, respectively. For both species pairs, we found the secondary contact with asymmetric gene flow from the deep to shallow species to be the best model (Tables 4   and 5).

| DISCUSS ION
The presumed lack of geographic barriers and propensity for high levels of gene flow poses challenges for studying speciation in the marine environment. Our results build on a growing body of work   .01 .04 [4,5,6,7,8,9,10,11,12,13] Group3 6.12 3.00 [arid1ab, arid1b, tfpt] rosenblatti comparison including cation homeostasis and the neuropeptide signaling pathway. The identification of genes in the neuropeptide signaling pathways in the S. chlorostictus-S. rosenblatti species pair is more directly related to ecological divergence and speciation in this group. Wang et al. (2021) point out the interplay between chemosensory divergence and ecological speciation. They focus on the importance of chemosensory drive in this process with a particular focus on diet adaptations. While dietary differences may exist in the S. chlorostictus-S. rosenblatti pair, it is likely adaptation to depth-related features is more important. Hyde and Vetter (2007) proposed a mechanism by which depth segregation could lead to divergence in Sebastes. In Sebastes, juveniles are attracted to speciesspecific habitat types during settlement, and homing to a different depth-related habitat could contribute to sensory drive that would eventually lead to reproductive isolation via habitat differences (Heras et al., 2015).
We can only speculate on the relative importance of these pathways to the ecological divergence of this species pair. Genes involved in homeostatic functioning are likely crucial in maintaining cellular stability in the face of environmental differences. While there has not been adequate characterization of depth-related habitat differences for any of the species we studied or the physiological demands of these environments, we can hypothesize that differences in temperature, pH, and salinity exist along the depth gradient that may drive local adaptation for these and other species.
Overall, the enrichment of candidate genes near significant F ST windows of genomic islands of differentiation for the S. chlorostictus-S.
rosenblatti pair provides functional descriptions of genes and gene networks related to behavior, development, homeostasis, and immunity, supporting previous work on rockfish that has found similar evolutionary evidence corresponding to depth driving their speciation (Behrens et al., 2021;Hyde & Vetter, 2007;Ingram, 2010).
We found fewer enriched gene ontology groups for the S. crocotulus-S. miniatus species pair, concordant with the finding of fewer significant outlier windows. This finding is likely due to the increased divergence of this species pair compared with S. chlorostictus and S. rosenblatti. It is possible that our approach is not applicable to more diverged species pairs and would benefit from methods that could account for levels of intraspecific variation (Cruickshank & Hahn, 2014). The amount of divergence between S. crocotulus and S. miniatus has likely eroded most of the signal associated with speciation due to the increased effects of drift and recombination (Quilodrán et al., 2020). The significant gene ontology terms that we did identify for this pair were associated with basic housekeeping functions (actin regulation of polymerization, the SWI/SNF pathways, histone lysine methylation pathways, tight junction, and regulation of actin filament polymerization) and cannot be associated with differences in depth for this species pair. The SWI/SNF pathways act as ATP-dependent chromatin remodelers that repress and activate genes and are associated with cardiovascular development (  (Mihola et al., 2009). Histone lysine methylation could indicate a postzygotic barrier in this species pair (Sha et al., 2020). It may be that the enriched genes found in this species pair are indicative of postzygotic barriers as chromatin remodeling genes are known to cause hybrid sterility and these molecular functions may be more present in longer diverged species of rockfish; however, additional work needs to be done to assess the validity of these findings.
We found a limited number of shared outlier windows between the two species pairs, a window on chromosome 6 and one on chromosome 12. These overlapping windows contained only four annotated genes (Table S4), some of which are involved in the basic metabolism. The lack of clear shared divergence islands across species pairs is likely due to the difference in estimated divergence times. This pattern was also observed by Behrens et al. (2021) in a comparison of genome-wide divergence in three Sebastes species pairs and demonstrates the limitations of using F ST in these types of studies (Quilodrán et al., 2020).

| Comparison to previous work
A recent study of the genomic architecture of speciation in Sebastes found evidence for two "islands" across two different chromosomes (Behrens et al., 2021). As in our study, Behrens et al. (2021)

| Demography of speciation
Examination of a limited set of demographic models indicated that the same model, secondary contact with asymmetric gene flow, was most likely for both species pairs. This suggests that the invasion of a novel habitat (deeper water in this case) is followed by a period of isolation in these species pairs. Models of ecological speciation predict that gene flow should persist upon invasion of the new ecological space. However, a recent study on depth-segregated ecomorphs in S. mentella found support for demographic models that were similar to those found in this study (Benestan et al., 2021). Benestan et al. (2021) suggested that divergence in S. mentella was relatively recent (0.5 MYA) and driven by changes in sea level during the Pleistocene. Support for the secondary contact model found in this study is also supported by other studies of speciation history in marine organisms (Fairweather et al., 2018;Filatov et al., 2021;Leder et al., 2021).
Another aspect of our analysis is that the directionality of asymmetric gene flow, following isolation, went consistently from the deeper species to the shallower species. It is unclear whether this pattern will hold across other depth-segregated species pairs in Sebastes, but could be indicative of climatic shifts impacting the depth distribution of these species. Future work on Northeastern Pacific Sebastes will determine whether the overall pattern of isolation followed by secondary contact holds across depthsegregated species pairs.

| Limitations
Our work is limited in that we utilized a pool-seq approach and only focused on enriched exome, CNE and UCE sequences. The main advantage of the pool-seq approach is that it reduces the overall cost of sequencing (Schlötterer et al., 2014). It does have the disadvantage that allele frequency estimates can be biased, but this is overcome with increased sequencing coverage (Schlötterer et al., 2015).
We intentionally utilized high coverage regions in our SNP discovery steps (40-500× coverage) to reduce any error. On top of this, we utilized enriched sequences in our analysis, which has the advantage of reducing sequencing efforts to protein-coding regions of the genome but would potentially be missing signals from extragenic regions. We were also limited in any inference from the comparison between the more divergent species pair (S. crocotulus-S. minatus), and future work in this area should focus on more recently derived species.

| CON CLUS IONS
Our exome scan of two Sebastes species pairs revealed a handful of genes and pathways associated with depth-related divergence.
There were a small number of shared islands of divergence between the pairs, but islands of divergence were more readily detected in the pair with more recent divergence. In the S. chlorostictus-S. rosenblatti pair, we found enrichment for the neuropeptide synthesis pathway in outlier loci, which suggests that chemosensory drive may be involved in depth-related speciation for this pair. Our analysis of demography of speciation revealed support for a similar model of divergence for the two pairs (isolation followed by secondary contact), which has been observed in other marine taxa. These results build on the growing knowledge of speciation history in the genus Sebastes and suggest Sebastes will continue to be a valuable model in understanding mechanisms of speciation in temperate marine fishes.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
All Illumina sequence data have been uploaded to the NCBI-SRA (project # PRJNA839756). All code for pool-seq analysis is located on: https://github.com/aagui l67/aacod e/blob/main/Sebas tes_Oliva res_1/Filte ring_Mappi ng_Popoo lation.