Studying the genetic basis of speciation in high gene flow marine invertebrates.

A growing number of genes responsible for reproductive incompatibilities between species (barrier loci) exhibit the signals of positive selection. However, the possibility that genes experiencing positive selection diverge early in speciation and commonly cause reproductive incompatibilities has not been systematically investigated on a genome-wide scale. Here, I outline a research program for studying the genetic basis of speciation in broadcast spawning marine invertebrates that uses a priori genome-wide information on a large, unbiased sample of genes tested for positive selection. A targeted sequence capture approach is proposed that scores single-nucleotide polymorphisms (SNPs) in widely separated species populations at an early stage of allopatric divergence. The targeted capture of both coding and non-coding sequences enables SNPs to be characterized at known locations across the genome and at genes with known selective or neutral histories. The neutral coding and non-coding SNPs provide robust background distributions for identifying FST-outliers within genes that can, in principle, identify specific mutations experiencing diversifying selection. If natural hybridization occurs between species, the neutral coding and non-coding SNPs can provide a neutral admixture model for genomic clines analyses aimed at finding genes exhibiting strong blocks to introgression. Strongylocentrotid sea urchins are used as a model system to outline the approach but it can be used for any group that has a complete reference genome available.

Introduction

Understanding how reproductive barriers evolve between populations represents one of the most fundamental challenges in evolutionary biology (Coyne and Orr 2004; Gavrilets 2004). Recent studies of genome-wide patterns of divergence between closely related species have provided new insights into the identity of genes causing reproductive barriers (Nadeau et al. 2012; Nolte et al. 2013), the semi-permeable nature of species barriers (Rieseberg et al. 1999; The ), and the genetic architecture of species differences (Teeter et al. 2010; Arnegard et al. 2014). Following the discovery of “genomic islands of speciation” in Anopheles mosquitoes by Turner et al, heterogeneous patterns of genome-wide divergence have been described in growing number of species pairs (reviewed by Sousa and Hey 2013; Seehausen et al. 2014). This heterogeneity was initially interpreted as reflecting highly variable levels of introgression across the genome (Via and West 2008; Nosil et al. 2009). However, the conditions producing genomic islands of divergence appear restrictive (Feder et al. 2012b; Flaxman et al. 2013) and recent re-analyses suggest that they might result from recent diversifying or background selection occurring during allopatric divergence (Cruickshank and Hahn 2014).

Many population genomic studies have embraced the concept of “ecological speciation” where diversifying selection represents the initial driver of speciation with, or without, ongoing gene flow (Schluter 2009; Nosil 2012). An important feature of ecological speciation is the reduced fitness of intermediate phenotypes (or migrants) caused by a mismatch between genotype and the environment (i.e., extrinsic postzygostic isolation). The enthusiasm for ecological speciation has shifted attention away from investigating how species originate in allopatry without any marked changes in ecology, physiology, or morphology. The widespread presence of cryptic species in both terrestrial and marine environments (Knowlton 1993; Pfenninger and Schwenk 2007) shows that reproductive barriers commonly evolve between taxa that remain unrecognized until genetic identification. The evolutionary processes and genes responsible for establishing and maintaining ecologically similar species in communities remain largely unknown. In the absence of strong ecological divergence, it is possible that some intrinsic isolation evolves from neutral or non-adaptive processes such as the random loss or movement of duplicated genes (Lynch and Force 2000) or intragenomic conflict (Rice 1998; Crespi and Nosil 2013).

Progress toward understanding the genetic basis of speciation requires identifying the genes responsible for reproductive barriers between species, the selective (or neutral) processes responsible for their divergence and their order of appearance (particularly during the initial stage) (Coyne and Orr 2004; The Marie Curie SPECIATION Network 2012). A growing list of “speciation genes” or “barrier loci” has been discovered that cause reproductive incompatibilities between species (Noor and Feder 2006; Presgraves 2010; Rieseberg and Blackman 2010; Nosil and Schluter 2011). An interesting pattern emerging from the characterization of barrier loci is that they often exhibit signals of positive Darwinian selection driven by genomic conflict (e.g., Ting et al. 1998; Phadnis and Orr 2009; Tang and Presgraves 2009; Nolte et al. 2013). Irrespective of the form(s) of selection involved, the possibility that barrier genes commonly have histories of positive selection has not been systematically investigated on a genome-wide scale. Theory predicts that Dobzhansky–Muller incompatibilities (DMIs; Dobzhansky 1937; Muller 1940) should be more prevalent at genes experiencing accelerated rates of divergence (Navarro and Barton 2003; Coyne and Orr 2004; Payseur and Nachman 2005). However, it is possible that many positively selected genes (PSGs) do not contribute to reproductive barriers and some may even experience adaptive introgression (e.g., Ding et al. 2014; Huerta-Sánchez et al. 2014). A powerful approach for distinguishing between PSGs conferring strictly adaptive benefits and those causing reproductive incompatibilities has been to study their patterns of introgression in hybrid zones (Harrison 1990; Payseur 2010).

In broadcast spawning marine invertebrates, the possibility that positive selection at gamete recognition proteins (GRPs) plays an important role in establishing reproductive barriers between species has received considerable attention (see reviews by Swanson and Vacquier 2002; Clark et al. 2006; Palumbi 2009; Lessios 2011; Vacquier and Swanson 2011; Kosman and Levitan 2014). However, the forms of selection responsible for positive selection on GRPs remain largely unknown and their levels of polymorphism and patterns of divergence vary unpredictably among groups (Vacquier and Swanson 2011; Kosman and Levitan 2014). The importance of divergence at GRPs during the early stage of speciation is also not clear because GRPs in allopatric species often fail to exhibit positive selection whereas those in sympatric taxa commonly do (Palumbi 2009; Vacquier and Swanson 2011; however, see Hart et al. 2014). Furthermore, strong asymmetrical patterns of gamete compatibility are common between closely related species (Zigler et al) that are incapable of providing effective barriers to introgression (Lessios 2007). For the sea urchin sperm protein bindin, gamete compatibility between species is correlated with amino acid divergence and Zigler et al. (2005) estimate that sister species can maintain compatibility for up to 5 million years. This is comparable to the mean time to speciation in many marine invertebrate groups (Coyne and Orr 2004), suggesting that more effective barriers must evolve earlier at other unknown genes or genetic elements.

The goal of this paper is to outline a research strategy for studying the genetic basis of speciation in broadcast spawning marine invertebrates that is based on a priori genome-wide knowledge of genes experiencing positive Darwinian selection. Although candidate gene approaches have been used to investigate genome-wide patterns of differentiation (e.g., Andrés et al. 2013; Hebert et al. 2013; Fraïsse et al. 2016), none has been based on a large, unbiased sample of genes with known selective (or neutral) histories. The advantage of this approach is that it provides testable predictions about the identities of genes that might diverge early in the speciation process and be more effectively blocked from introgressing upon secondary contact. This avoids the post hoc characterization of “genomic islands of divergence” that often fail to identify the presumed barrier loci due to unknown linkage associations with the selected gene(s). Strongylocentrotid sea urchins will be used as a representative group to highlight the approach, but it can be applied to any system provided a well-annotated reference genome is available.

The Geography of Speciation in High Gene Flow Marine Invertebrates

The importance of geography in speciation has recently been de-emphasized and replaced with models based on a continuum of gene flow (e.g., Butlin et al. 2008; Fitzpatrick et al. 2009; Harrison 2012). This shift in perspective is necessary for nearshore marine invertebrates whose geographic ranges have undergone dramatic and dynamic changes over the Pleistocene caused by glacial cycles (Jansson and Dynesius 2002; Jacobs et al. 2004; Norris and Hull 2012). In the northeast Pacific region, the fossil record shows that nearshore species respond individualistically to changing ocean temperatures and sea levels, and that community composition is in a state of perpetual flux (reviewed by Valentine and Jablonski 1993; Lindberg and Lipps 1996). The dynamic and complex histories of species past distributions make it difficult to assign a strict geographic mode of speciation, to infer the selective agents causing adaptive divergence, or to assess the role played by reinforcement. The rapidity of environmental change over the latter part of the Pleistocene has apparently not resulted in increased rates of speciation or extinction (Jackson and Johnson 2000).

The high levels of gene flow that occur among populations of many marine invertebrates with long-lived pelagic larvae simplify studying some aspects of their speciation while complicating others. The weak or non-existent population structure commonly observed in these species renders models of “speciation-with-gene-flow” (Feder et al. 2012a; Via 2012) unlikely because extensive gene flow overpowers diversifying selection across most of the genome. If selection does overcome gene flow and cause local adaptation (Sanford and Kelly 2011) it is unclear whether this eventually results in the formation of isolating barriers. The initiation of speciation in high gene flow marine taxa must involve a period of allopatric divergence following long-distance colonization across or between ocean basins (Lindberg 1991; Vermeij 1991) or by major vicariance events such as the rising of the Isthmus of Panama (Lessios 2008). Under these scenarios, complications arise when attempting to identify the genes diverging early in the speciation continuum. Here, the challenge is to identify recently formed allopatric populations separated by large geographic distances (i.e., different ocean basins) that have undergone a small amount of genetic divergence. Although young allopatric populations may not eventually form new species, they do provide snapshots of the initial adaptive divergence following long-distance dispersal events, which is the starting point for how many marine species originate.

The Targeted PSGs Strategy

Here, I propose a research strategy for testing the hypotheses that PSGs (i) diverge early among allopatric populations of a species and (ii) commonly evolve to become barrier loci. This strategy is based on knowing a priori the identity of genes experiencing positive selection across the genome in a group of closely related species. There are 4 requirements for this approach to be successful. First, a well-annotated reference genome must be available for what I call a “focal species”. This is clearly a stringent requirement that excludes virtually all species. However, the advantages provided by a reference genome for studying speciation are substantial (see below) and complete annotated genomes are now available for a growing number of marine invertebrate species (see Table 1). Second, a group of closely related species must be present whose genomes can be sequenced using the focal species’ genome as a reference. The generation of multi-species alignments of single-copy orthologs enables powerful and robust genome-wide tests for positive selection using programs like PAML (Yang 2007) or HyPhy (Kosakovsky Pond et al. 2005). A sufficient number of related taxa must be present (> 6) at an appropriate scale of divergence to provide adequate statistical power (Anisimova et al. 2001). Third, conspecific populations must exist at an early stage of divergence to allow the targeted capture and sequencing of population samples of all single-copy orthologs previously tested for positive selection. This will allow the identification of genes experiencing recent diversifying selection in allopatry (i.e., FST-outliers) and provide insights into the overlap between historical and contemporary selection. Finally, low levels of introgression must occur between some species that enable the identification of barrier loci that are blocked from moving between species. Ideally, multiple hybrid zones or regions experiencing introgression are present to test for the repeatability of introgression or blockage (Payseur 2010). Meeting all of these requirements is obviously extremely challenging. However, strongylocentrotid sea urchins represent one group where all of the above conditions are met and they will be used to illustrate some of the challenges of the approach described in this paper. Other groups with fairly well-resolved phylogenic histories that could be studied in a similar manner include the genus Octopus (Söller et al. 2000; Guzik et al. 2005) and the genus Crassostrea (Yu and Li 2012; Trivedi et al. 2014).

Table 1

Complete genomes of marine invertebrate species with available gene models

Phylum	Species name	Size (Mb)	GC (%)	Number of genes	Reference(s)
Placozoa	Trichoplax adhaerens	105.6	32.7	11,518	Srivastava et al. (2008)
Porifera	Amphimedon queenslandica	166.7	37.5	13,998	Srivastava et al. (2010)
Mollusca	Crassostrea gigas	557.7	35.3	32,261	Zhang et al. (2012)
Echinodermata	Strongylocentrotus pupuratus	814.0	36.9	21,092	Sodergren et al. (2006), Tu et al. (2012)
Cephalochordata	Branchiostoma floridae	521.9	41.2	28,627	Putnam et al. (2008)
Tunicata	Oikopleura dioica	70.5	39.8	17,212	Seo et al. (2001)
Tunicata	Ciona intestinalis	115.2	33.1	15,254	Dehal et al. (2002)
Brachiopoda	Lingula anatina	425.5	36.9	34,000	Luo et al. (2015)
Anthozoa	Nematostella vectensis	356.6	40.6	27,173	Putnam et al. (2007)
Cephalopoda	Octopus bimaculoides	2,282.8	36.1	32,819	Albertin et al. (2015)
Ctenophora	Mnemiopsis leidyi	155.9	38.9	16,548	Ryan et al. (2013)
Annelida	Capitella teleta	333.3	40.4	31,997	Simakov et al. (2013)
Mollusca	Lottia gigantea	359.5	33.3	23,827	Simakov et al. (2013)
Mollusca	Aplysia californica	927.3	42.0	21,312	https://www.broadinstitute.org/scientific-community/ science/projects/mammals-models/vertebrates- invertebrates/aplysia/aplysia-genom
Priapulida	Priapulus caudatus	511.7	45.7	17,096	http://www.ncbi.nlm.nih.gov/genome/?term=Priapulus+caudatus
Hemichordata	Saccoglossus kowalevskii	775.8	38.1	22,073	https://www.hgsc.bcm.edu/acorn-worm-genome-project
Arthropoda	Limulus polyphemus	1,828.3	34.5	22,031	http://genome.wustl.edu/genomes/detail/limulus-polyphemus/

aCompiled from the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov) in December 2015.

Genome-Wide Screens for Positive Darwinian Selection

Genome-wide tests for genes experiencing positive selection have been conducted on only a limited number of groups including Drosophila (), mammals (Kosiol et al. 2008), birds (Künstner et al. 2010), ants (Roux et al. 2014), and sea urchins (Kober KM, Pogson GH, submitted for publication). Although a number of factors affect the success of genome-wide screens for positive selection, one of the most foundational is that a set of accurate gene models must be available. The well-annotated genome of the purple sea urchin Strongylocentrotus purpuratus provides a good example of this necessity. The S. purpuratus genome contains ∼23,300 genes of which over 9,000 were manually annotated by an international consortium of over 200 members (Sodergren et al. 2006). The initial gene predictions were recently updated by Tu et al. (2012) who performed a deep transcriptome analysis on pooled RNA from 10 different embryonic stages, 6 different larval stages, and 6 adult tissues. These updates incorporated exons missing from the original gene models, excised erroneously included exons, corrected inaccurate exon/intron boundaries, and removed overlapping gene models (Tu et al. 2012). Surprisingly, over half of the original manually curated gene models contained errors. These results highlight the need for direct experimental evidence (i.e., deep transcriptome sequencing) to validate computational gene predictions. The accuracy of the S. purpuratus gene models are on par with those from well-studied model species (humans, mice, Drosophila, and Caenorhabditis elegans) making sea urchins an excellent system for genome studies.

Using the genome of S. purpuratus as a reference, the complete genomes of 8 additional strongylocentrotid sea urchins have recently been obtained (Kober and Bernardi 2013a; Kober KM, Pogson GH, submitted for publication). Illumina paired-end reads (100 bp) were aligned directly to the reference genome, producing high coverage across single-copy protein-coding genes (mean = 58.2X). Although gene tree topologies were heterogeneous, the species tree was clearly resolved from 4-fold degenerate sites at 2,815 genes not experiencing positive selection (see Figure 1; Kober and Bernardi 2013a). By fitting PAML models M7 and M8 (Yang 2007) and applying a conservative false discovery rate of 5%, 1,008 of the 6,520 single-copy orthologs tested (15.5%) showed significant signals of positive selection. A further 824 candidate PSGs were found to exhibit positive selection along the 9 terminal branches of each gene tree (so-called “branch-sites” tests). Recent substitutions along these terminal branches represent prime candidates for barrier loci because of their potential negative epistatic effects (Orr 1995; Figure 2). Consistent with the action of repeated selective sweeps, the candidate PSGs had 68% higher rates of nonsynonymous substitution (dN) and 33% lower levels of heterozygosity than non-PSGs (both groups had similar rates of synonymous substitution, dS). As expected, positive selection was observed at innate immunity genes and reproductive proteins. However, it was also prevalent at membrane receptors, cell adhesion molecules, extracellular matrix proteins, and ion channels that have been implicated as targets of pathogens in other groups (notably mammals). The PSGs represent excellent a priori candidates that might experience diversifying selection in recently formed allopatric populations.

Figure 1.

Species tree of the family Strongylocentrotidae constructed from 4-fold degenerate sites from 2,815 genes not experiencing positive selection (adapted from Kober and Bernardi 2013a). Trees reconstructed using maximum parsimony (MP), maximum likelihood (ML), and Bayesian approaches produced identical topologies. All nodes had MP bootstrap values of 100 and Bayesian poster probabilities of 1. Abbreviations denote geographic locations of species distributions (NEP = Northeastern Pacific, NWP = Northwestern Pacific; CIR = circumpolar).

Figure 2.

The distinction between historical and recent positive selection. (A) Alignment of portion of a hypothetical protein-coding sequence with one codon (red) showing historical positive selection and another (green) showing recent positive selection based on PAML (Yang 2007) sites tests and branch-sites tests, respectively. Historical and recent selections are not mutually exclusive. However, when both operate it is difficult to produce a significant branch-site test because the external branch needs to have a significantly elevated dN/dS ratio compared with the remainder of the tree. (B) Gene tree with locations of amino acid changes marked. Historical selection at the red codon cannot generate incompatibilities between S. droebachiensis and S. pallidus but recent selection at the green codon can (designated by the “X”).

Caution must be exercised when interpreting genome-wide tests for positive selection because errors in sequencing, alignment, and annotation are known to generate significant numbers of false positives (Schneider et al. 2009; Fletcher and Yang 2010; Markova-Raina and Petrov 2011; Jordan and Goldman 2012). Although sequencing errors can be minimized by applying strict base quality metrics and high coverage cut-offs, mistakes in gene annotations and alignments are more difficult to address. To minimize annotation and alignment errors, a conservative reference-based alignment method should be used that excludes all codons not present in the reference genome (caused by insertions) and removes all missing codons (caused by deletions) prior to testing for positive selection. This should largely avoid alignment errors caused by gains or losses of amino acids. In Drosophila, Markova-Raina and Petrov (2011) observed that alignment and annotation errors were common near exon borders and regions with insertions or deletions of amino acids that resulted in significant spurious signals of positive selection. Given these results, it is essential that all genome-wide screens for positive selection test for the possible enrichment of false positives in these regions. In sea urchins, Kober and Pogson (submitted for publication) subjected 3, 5, and 7 amino acid windows adjacent to exon borders and missing data regions for enrichment tests and obtained negative results. These results are important because they eliminate the possibility that the signals of positive selection are largely caused by annotation or alignment errors and provide further testament to the accuracy of the gene models.

Targeted Genome Scans for FST-Outliers

In many high gene flow marine invertebrate species inhabiting a single geographic region, it is problematic to test for FST-outliers because of the limited population genetic structure present. However, FST-outlier tests can be performed on marine invertebrate taxa that have broad geographic distributions ranges that may span different ocean basins. Here, populations often exhibit moderate levels of genetic differentiation when sampled at large spatial scales. For example, 2 strongylocentrotid species (Strongylocentrotus droebachiensis and Strongylocentrotus pallidus) have broad arctic–boreal distributions across the north Atlantic and north Pacific oceans (Mortensen 1943; Jensen 1974). No population structure has been described among populations of S. droebachiensis within the same geographic region but moderate genetic differentiation occurs between ocean basins (Palumbi and Wilson 1990; Addison and Hart 2004, 2005). Populations of S. droebachiensis from the northeast Pacific, northeast Atlantic, and northwest Atlantic exhibit significant pair-wise FST-values at the mitochondrial COI locus ranging from 0.213 to 0.325 (Addison and Hart 2005). For microsatellites, these same populations exhibit a mean overall FST-value of 0.087 (Addison and Hart 2004). Comparable studies have not been conducted on S. pallidus but its similar geographic distribution and pelagic larval duration to S. droebachiensis has likely resulted in similar levels of population divergence. Recently, genome-wide single-nucleotide polymorphisms (SNPs) were compared between mussel populations (Mytilus trossulus and Mytilus edulis) sampled from either side of the north Atlantic and were found to exhibit pair-wise FST-values ranging between 0.131 and 0.362 (Fraïsse et al. 2016). These results suggest that other marine invertebrates with similar broad distributions might be amenable for FST-outlier studies aimed at identifying recently diverged loci.

Excluding humans, the majority of genome-wide scans for FST-outliers have used randomly located SNPs, implicitly assuming that patterns of polymorphism across the genome reflect neutral and demographic processes (reviewed by Oleksyk et al. 2010; Haasl and Payseur 2016). One of the long-standing problems faced by FST-outlier studies is the difficulty of accurately quantifying the neutral FST distribution. Some have questioned whether, in some species, a neutral background actually exists (Comeron 2014; Burri et al. 2015). The large effective population sizes of many marine invertebrates may have resulted in the pervasive influence of linked selection across the genome, similar to that observed in Drosophila (Begun et al. 2007; Langley et al. 2012; Comeron 2014). Targeting protein-coding genes for genome scans confounds this problem because these regions have the potential to directly experience diversifying selection, background selection, balancing selection, and even weak selection on patterns of synonymous codon usage (see Kober and Pogson 2013b).

With a reference genome and prior information on genes known to have experienced positive selection there are several ways to estimate a neutral FST distribution. One is to score SNPs at single-copy genes that do not exhibit positive selection (non-PSGs) or exhibit extremely low rates of nonsynonymous substitution. Another approach would be to score thousands of SNPs in non-coding, non-repetitive intergenic regions at large distances (e.g., >50 kb) from any known structural gene. Here, it would be possible to avoid highly conserved non-coding elements, transcription factor-binding sites, and to exclude all regions containing any PSGs. It would be worthwhile evaluating whether null distributions provided by non-PSG SNPs or those based on non-coding SNPs provide the most robust identification of FST-outliers. Because of differences in GC content, recombination rates may differ between non-coding and coding regions of the genome (Eyre-Walker 1993). The neutral FST distribution based on non-PSG SNPs might thus result in a more accurate detection of genes experiencing recent diversifying selection. Irrespective of this outcome, it will also be necessary to confirm the outlier status of SNPs using an absolute measure of divergence (such as d) to avoid the spurious inflation of FST-values caused by reduced levels of polymorphism (see Noor and Bennett 2009; Cruickshank and Hahn 2014).

As the cost of next-generation sequencing continues to fall, it is tempting to collect complete genome data for SNP identification and calling. However, because protein-coding genes are the focus of this research strategy, a more cost-effective approach would be to use high-throughput hybridization-based targeted sequence capture of exons (reviewed by Mamanova et al. 2010). Exon capture has proven to be a cost-effective genotype-by-sequencing method that, although widely used in human genomics, has not been widely applied to non-model species (see Jones and Good 2016). Many studies have shown that sequence capture using RNA “baits” generates longer contigs, produces higher quality SNP calls, and provides more even coverage than other genotype-by-sequencing methods like Rad-Seq (Baird et al. 2008) or variants thereof (see Tewhey et al. 2009; Harvey et al. 2013). Although some of the SNPs scored at non-PSGs or in non-coding regions may not be truly neutral, they are expected to provide a conservative null distribution because linked selection will act to inflate their variance and make FST-outlier calls more conservative.

FST-outlier studies are susceptible to false positive signals of diversifying selection due to a variety of factors including population bottlenecks (Foll and Gaggiotti 2008), hierarchical population structure (Excoffier et al. 2009; Fourcade et al. 2013), cryptic hybrid zones (Bierne et al. 2011), and range expansions (Hofer et al. 2009). The proposed targeted capture of SNPs in both coding and non-coding regions of the genome offers a number of advantages over standard genome scans. First, the ability to detect FST-outliers should be significantly improved over the practice of scoring random SNPs in unknown locations. Rather than relying on linkage disequilibrium between a randomly located SNP and a selected gene, the protein-coding regions captured are potentially the direct targets of selection. Second, characterization of a large number of putatively neutral SNPs enables application of the empirical P-value approach described by Lotterhos and Whitlock (2014) where the significance of an FST-outlier is based on its quantile in the empirical distribution of the neutral SNPs (Figure 3). Lotterhos and Whitlock (2014) found that this empirical P-value approach significantly reduced false positives, even under non-equilibrium scenarios such as recent range expansions. Third, the FST-outliers in coding regions can be more easily traced to specific mutations at genes with known selective histories and the impact of linked selection in these regions can be examined in detail (reviewed by Sousa and Hey 2013; Vatsiou et al. 2016). Fourth, alleles that have introgressed between species can be identified by reconstructing gene genealogies and by applying powerful new tests for detecting introgression (e.g., Geneva et al. 2015; Rosenzweig et al. 2016). Discovering that genes with histories of positive selection have successfully moved between species might provide evidence favoring adaptive introgression (Hedrick 2013). Finally, the neutral SNPs could provide robust insights into the species demographic histories (Excoffier et al. 2013; Liu and Fu 2015), which can provide important insights when interpreting the patterns of diversifying selection.

Figure 3.

Identification of FST-outliers by the targeted sequence capture of non-coding and coding SNPs. (A) Two gene models are shown on a scaffold that are targeted for solution-based exon capture. Two non-coding regions >10 kb from either gene are also captured (but could be much further away). (B) Hypothetical distributions of pairwise FST-values for genome-wide coding and non-coding SNPs are presented. The arrow designates the 5% false discovery rate cut-off estimated directly from the non-coding SNPs. Any genes or gene regions with FST-values beyond this threshold are candidates for recent diversifying selection. Null distributions for FST-outlier tests could also be generated from the subset of protein-coding genes that do not exhibit histories of positive selection.

The FST-outliers detected among widely separated populations of marine invertebrates could provide important insights into the genes that first diverge following the colonization of new habitats by long-distance dispersal. A broad overlap between the FST-outlier genes and those known a priori to have experienced positive selection would confirm that diversifying selection is ongoing and that there is a tight coupling of historical and contemporary selection. This might be expected for genes experiencing sexual selection or conflict where relentless cycles of adaptation and counter-adaptation are expected. However, the link between selection at different time scales may be weak. If so, this might suggest that novel forms of selection are being imposed on the new allopatric population, such as that mediated by new suites of pathogens. However, the detection of diversifying selection at a locus does not provide any insight into whether it causes a genetic barrier between diverging lineages. This can only be assessed by performing controlled crosses or by examining the patterns of naturally occurring introgression.

Genome Scans for Barrier Loci

Hybridization is common in the marine environment, occurring at similar levels to that described in terrestrial systems (reviewed by Gardner 1997; Arnold and Fogarty 2009). In his review, Gardner (1997) observed that only 34 of the 108 documented cases of marine hybridization involved hybrid zones or hybrid swarms. The majority of hybridization events in the sea thus occur in the absence of distinct hybrid zones. Many hybridizing broadcast-spawners have broadly overlapping geographic ranges and heterospecific matings must occur randomly when breeding coincides in space and time. Such appears to be the case for S. droebachiensis and S. pallidus. The reproductive barriers between the 2 taxa appear fairly well established and F1 hybrids in nature are rare. For example, using morphological characters Vasseur (1952) identified 3 hybrids in a sample of 562 individuals (0.54%) from northern Norway. Although these observations suggest that gene flow between the two species is limited, even small amounts of backcrossing can result in significant introgression (Barton 2001). Low levels of asymmetric introgression has been recently observed from S. pallidus into S. droebachiensis at a small number of nuclear genes in populations from the northeastern Pacific region (Addison and Pogson 2009; Pujolar and Pogson 2011). The prevalence and identity of genes experiencing introgression across the genomes of this species pair is unknown.

The genomic clines method of Gompert and Buerkle (2009, 2010) can be used to study the patterns of introgression between diverged lineages to identify putative barrier loci. This approach uses multinomial regression to estimate the probability of a specific genotype at a marker given an observed level of genome-wide admixture (i.e., a hybrid index). The genomic clines method relies on estimating a model of neutral introgression. This is done by randomly permuting genotypes among loci within individuals, or by using genetic data from populations outside of the hybrid zone to estimate the probability of genotypes in the admixed population (see Gompert and Buerkle 2009). Similar to the problem faced by estimating the neutral FST distribution, using SNPs from protein-coding genes to estimate the neutral admixture model is inappropriate. However, the random permutation of non-PSG or non-coding SNP genotypes among markers within individuals should provide an accurate characterization of the neutral pattern of admixture. This in turn should increase the power for identifying barrier loci that exhibit significantly reduced levels of introgression between species.

A growing number of hybrid zone studies have reported highly variable degrees of hybridization and patterns of introgression in different geographic regions (e.g., Nolte et al. 2009; Teeter et al. 2010; Mandeville et al. 2015). Interpreting heterogeneous patterns of introgression is particularly challenging because it can result from geographic differences in a number of pre- or post-zygotic isolating mechanisms, changes in relative abundances, or from the presence of population-specific barriers. Since most broadcast spawning marine invertebrates lack strong mate choice, one of the most basic ways for patterns of introgression to vary among geographic regions are through differences in the degree of spatial or temporal isolation. If this information is available, it can be used to predict how introgression may differ between geographic regions. For example, in the northeast Pacific, the reproductive cycles of S. droebachiensis and S. pallidus overlap substantially (Levitan 1998) and F1 individuals produced from these populations in the laboratory are viable and fertile (Strathmann 1981). In the northeastern Atlantic, there is greater spawning asynchrony between S. droebachiensis and S. pallidus, with the former spawning earlier in the year (Vasseur 1952; Hagström and Lönning 1967). Therefore, in this system one might predict that, all else being equal, higher levels of introgression might occur in the northeast Pacific than the northeast Atlantic.

Another factor complicating patterns of introgression among populations is that the distributions of closely related species may vary geographically. Recent studies on Lake Victoria cichlids (Keller et al. 2013), Colorado River catostomid fishes (Mandeville et al. 2015), and blue mussels belonging to the genus Mytilus (Fraïsse et al. 2016) have also illustrated how hybridization and introgression often needs to be assessed from a multi-species perspective. It is becoming increasingly clear that a geographical context is essential for properly interpreting patterns of introgression. In strongylocentrotid sea urchins from the northeast and northwest Atlantic, heterospecific matings are only possible between S. droebachiensis and S. pallidus. However, in the north Pacific different assemblages of species broadly co-occur along the Asian and North American coasts (see Figure 1). If reproductive barriers are incomplete among other species, the patterns of introgression might vary in different and unpredictable ways across the north Pacific.

Some Limitations

The research strategy described here for evaluating the role of positive selection in producing barrier loci faces several important limitations. First, tests for positive selection based on dN/dS ratios are known to be conservative (Anisimova et al. 2001; Gharib and Robinson-Rechavi 2013; Lu and Guindon 2013). It is thus likely that many genes undergoing adaptive diversification fail to exhibit statistically significant signals of positive selection and their omission would result in an unknown number of false negatives. Less stringent criteria could be used to classify candidate barrier genes such as those exhibiting the highest rates of nonsynonymous substitution, irrespective of positive selection per se. Second, recently formed paralogous genes create challenges for genome assemblies and alignments and thus will be underrepresented in the tests for barrier loci that use conserved single-copy orthologs. The importance of new genes in forming reproductive barriers between species is largely unknown, but they appear to diverge quickly and often acquire new functions (Chen et al. 2013; Long et al. 2013). Finally, although candidate barrier genes may be successfully identified, the forms of selection responsible for their divergence may be difficult to discern and additional functional tests (and genetic crosses) are needed to confirm their direct involvement in establishing species barriers. In some cases, the identities of barrier loci might provide clear insights into a specific selective agent. For example, observing that innate immunity genes often experience positive selection and are strongly blocked from introgressing would implicate pathogens as the drivers of divergence. Distinguishing between pre- or post-zygotic barriers will also be challenging, but in some cases can also be discerned from gene identity. For example, a strong barrier at a GRP would clearly implicate a post-mating prezygotic barrier, but the form of selection responsible remains unknown.

Broadcast spawning marine invertebrates often possess extensive levels of genetic polymorphism that may complicate some of the analyses proposed here. Multi-species alignments of single-copy genes will contain large numbers of heterozygous mutations that must be filtered prior to testing for positive selection. This is because many mutations represent low frequency nonsynonymous changes (likely deleterious) that, if included, would lead to overestimation of dN and dN/dS ratios and generate false signals of positive selection. High levels of heterozygosity can result in the loss of considerable data. For example, at the 6,520 single-copy genes studied by Kober and Pogson (submitted for publication) in sea urchins, codons containing 1.72 million high-quality SNPs were filtered prior to testing for positive selection. High levels of standing variation may also lead to the confounding of divergence with polymorphism if only one individual is sequenced for each species. Sequencing the transcriptomes of several individuals per species might improve the identification of fixed differences and allow McDonald–Kreitman tests (McDonald and Kreitman 1991) to be performed that would complement the tests for selection based on d/d ratios.

Long-distance colonization events involving marine invertebrate populations with high levels of standing variation also leads to the prediction that soft sweeps (Hermisson and Pennings 2005; Pritchard et al. 2010) could be more common than hard sweeps. Since soft sweeps generate diminished signals of genetic hitchhiking, the ability to detect FST-outliers may be challenging especially in regions with high rates of recombination. Interpreting soft sweeps may be greatly improved if the biogeographical history of the group is known and if the beneficial alleles are still segregating in known ancestral populations. In the case of S. droebachiensis and S. pallidus, populations in across the north Atlantic are known to have originated from the north Pacific (Durham and MacNeil 1967; Addison and Hart 2005). It is thus possible that haplotypes experiencing sweeps in the north Atlantic may be identifiable in populations from the north Pacific.

Concluding Remarks

Confirming that a diverse array of genes experiencing positive selection diverge early in the speciation process and commonly result in reproductive incompatibilities would represent an important step forward in understanding of the genetic basis of speciation. The power of the approach outlined here results from the targeted capture of both con-coding and coding SNPs in known locations and at genes with known selective (or neutral) histories. The neutral coding and non-coding SNPs provide more accurate estimations of the neutral FST distribution for FST-outlier tests and neutral admixture rates for the genomic clines analyses identifying barrier loci. It is also worth noting that negative results would be equally important to speciation research. Discovering that positive selection does not typically produce barrier loci would direct attention to genetic elements in non-coding regions responsible for establishing and maintaining species boundaries (e.g., Jones et al. 2012; Nadeau et al. 2012; Hebert et al. 2013). Moving beyond the post hoc characterization of genome-wide patterns of divergence to a hypothesis testing framework targeting known genes offers many advantages. It can provide direct and detailed insights into the identity of genes responsible for establishing species boundaries and the importance of gene flow in preventing or facilitating the adaptive divergence of lineages (Abbott et al. 2013). Although the ability to simultaneously study early divergence and patterns of introgression between young marine invertebrate species is clearly challenging, the approach outlined here might provide an exciting way forward.