Speciation with gene flow in marine systems

Over the last century, a large body of literature emerged on mechanisms driving speciation. Most of the research into these questions focussed on terrestrial systems, while research in marine systems lagged behind. Here, we review the population genetic mechanisms and geographic context of 33 potential cases of speciation with gene flow in the marine realm, using six criteria inferred from theoretical models of speciation. Speciation with gene flow occurs in a wide range of marine taxa. Single traits, which induce assortative mating and are subjected to disruptive selection, such as differences in host-associations in invertebrates or colour pattern in tropical fish, are potentially responsible for a decrease in gene flow and may be driving divergence in the majority of cases. However, much remains unknown, and with the current knowledge, the frequency of ecological speciation with gene flow in marine systems remains difficult to estimate. Standardized, generally applicable statistical methods, explicitly testing different hypotheses of speciation, are, going forward, required to confidently infer speciation with gene flow.


Introduction
Since the modern synthesis of biology, speciation without gene flow has been accepted as the default mode of speciation. An extrinsic barrier splits a population and prevents gene flow between subpopulations, while genetic incompatibilities between populations accumulate through genetic drift and gradually result in reproductive isolation (Mayr, 1942(Mayr, , 1963. Whether reproductive isolation can arise without extrinsic barriers has long been a discussion in the field of evolutionary biology. In the absence of, for example, geographic barriers, interbreeding may impede divergence within a population, which was already recognized by Wagner (1868) in his critique on Darwin (1859). However, since the second half of the twentieth century, mathematical models of speciation with gene flow driven by disruptive selection emerged (Maynard Smith, 1966;Dickinson & Antonovics, 1973;Rice, 1984Rice, , 1987Johnson et al., 1996;Kawecki, 1996Kawecki, , 1997Dieckmann & Doebeli, 1999;Kondrashov & Kondrashov, 1999;Doebeli & Dieckmann, 2000;Gavrilets, 2004). Empirical examples of divergence in sympatry, driven by disruptive selection (i.e., divergent selection within a single, panmictic population; Rundle & Nosil, 2005), were reported as well (Bush, 1969;Feder et al., 1988;Rice & Salt, 1988). The paradigm of speciation in isolation as the sole mechanism of speciation slowly shifted towards a view incorporating speciation with gene flow as a plausible alternative (Bolnick & Fitzpatrick, 2007). Currently, it is accepted that speciation can occur in the presence of gene flow (e.g., Hey, 2006;Bolnick & Fitzpatrick, 2007;Nosil, 2008). Literature has mostly focused on the biogeographic patterns, e.g., by comparing the frequency of allopatric vs. sympatric speciation. The focus has recently shifted towards the population genetic process (Fitzpatrick et al., , 2009), where we focus on here.
Most research so far has dealt with speciation in terrestrial systems, while research on evolutionary mechanisms and speciation in the marine realm has been lagging behind (Miglietta et al., 2011;Peijnenburg & Goetze, 2013;Bowen et al., 2016). Based on a study on tropical echinoids, Mayr (1954) concluded that speciation in marine systems did not differ from speciation in terrestrial systems. However, geographic barriers are seemingly less common in marine systems compared to terrestrial systems, and many organisms have a large dispersal potential because of long pelagic larval stages (Palumbi, 1992;Hellberg, 1998;Kinlan & Gaines, 2003;Rocha et al., 2005Rocha et al., , 2007Rocha & Bowen, 2008;Miglietta et al., 2011; but see Chen & Hare, 2011;Peijnenburg & Goetze, 2013;Goetze et al., 2017). Despite the lack of obvious geographic barriers, diversity in marine systems can be very high. The shallow-waters of the tropical latitudes, for example, are amongst the biologically most diverse ecosystems, in some cases rivalling the biodiversity seen in tropical rain forests (Reaka-Kudla, 1997;Plaisance et al., 2011;Fisher et al., 2015). Some authors therefore hypothesised that speciation with gene flow plays a larger role in marine systems compared to terrestrial ecosystems (Miglietta et al., 2011;Bowen et al., 2013). Here, we review potential cases of speciation with gene flow in marine systems. We discuss and compare the potential mechanisms driving divergence with gene flow among different case studies. We further discuss gaps in the current state of knowledge and suggest directions for future research on speciation with gene flow in marine systems.

Mechanisms of speciation with gene flow
Mathematical models of speciation with gene flow mostly follow a similar pattern. Disruptive selection in an initially panmictic population results in non-random mating and reproductive isolation within that population (Maynard Smith, 1966;Bolnick & Fitzpatrick, 2007). With disruptive selection being the driver of divergence, speciation with gene flow can in most cases be considered a special case of ecological speciation, as defined by Schluter (2009). It should be noted though that ecological factors and selection are likely to play a role in other mechanisms of speciation as well (Schilthuizen, 2000;Schilthuizen et al., 2006;Sobel et al., 2010). Disruptive selection on its own is, however, not sufficient to maintain reproductive isolation, as random mating will counteract the genetic diversification and reproductive isolation between subpopulations (Bolnick & Fitzpatrick, 2007). Assortative mating, for example through mate selection or (sub-) habitat differentiation, is required. Recombination between the assortment trait and the selected trait would, however, prevent the evolution of reproductive isolation through disruptive selection (Felsenstein, 1981), which has been coined the 'selection-recombination antagonism' by Rice (1984). Selection should therefore either indirectly, in linkage disequilibrium, or directly affect the assortment trait causing non-random mating for the evolution of reproductive isolation to occur (Kirkpatrick & Ravigné, 2002). Besides evidence for (historical) gene flow after the divergence between populations, evidence for disruptive selection, assortative mating, and a link between the selected trait and the assortment trait should be present in the studied populations to infer speciation with gene flow. We will evaluate cases of potential (incipient) speciation with gene flow against six criteria. First, to identify cases of incipient speciation we ask: 1.
Are populations reproductively isolated? This is the second criterion of Coyne & Orr (2004) for speciation in sympatry. While the fate of cases of incipient speciation is not known, studying systems at different stages of speciation could provide a better insight into the mechanisms playing a role during the process of speciation (Butlin et al., 2008). It is likely that some hybridisation will occur among closely related species. The threshold value of hybridisation, below which reproductive isolation is considered to be 'complete' , is impossible to determine exactly. Determining such a value will always be arbitrary, trying to classify the continuous scale of hybridisation, from panmixia to full reproductive isolation, into discrete categories. Additionally, data on hybridisation are often scarce in less well-studied species, and therefore difficult to quantify and compare among species. This makes assessing this criterion problematic. However, the effect of this lack of data will not affect the conclusions of this study, as the next criteria are more relevant to the question of the frequency of ecological speciation with gene flow.
Secondly, the four requirements for speciation with gene flow are assessed: 2. Is there (potential for) disruptive selection? 3. Do populations mate assortatively in sympatry? 4. Is the selected trait linked to the assortment trait? 5. Is there evidence for gene flow between populations at the time of divergence? These criteria are derived from mathematical models of speciation with gene flow (e.g., Maynard Smith, 1966;Kawecki, 1997;Doebeli & Dieckmann, 2000). The second, third, and fourth criterion were also used by Rundle & Nosil (2005) in their review of ecological speciation. The fifth criterion, which directly tests for the occurrence of gene flow during divergence, was added to exclude cases where ecological speciation occurred without gene flow, in geographically separated populations (Schluter, 2009). However, this criterion is problematic to meet, as it requires the reconstruction of the degree of historical gene flow based on data from modern populations.
Finally, the geographic context of speciation will be discussed: 6. Do geographic ranges of populations overlap? This is the first criterion of Coyne & Orr (2004) for speciation in sympatry. The geographic context of speciation (sympatry versus allopatry) and the degree of gene flow may often be correlated in the field but are not identical and should therefore not be used interchangeably. While the geography of speciation cannot be considered an evolutionary process (Fitzpatrick et al., , 2009, the question how geography influences gene flow is still interesting (Mallet et al., 2009;Maas et al., 2018), especially in marine systems without clear geographic borders, even though recent studies show that in the open ocean geographic as well as bathymetric barriers might be more common than previously assumed (Weiner et al., 2012;Peijnenburg & Goetze, 2013).

Potential cases of speciation with gene flow
We searched the literature for potential cases of ecological speciation with gene flow in marine systems. Studies were included if evidence for gene flow between sister clades (either conspecific, genetically diverged populations or distinct species), was found or if sister clades occur in sympatry. The literature search yielded a list of 33 marine taxa across a wide taxonomic range in which ecological divergence with gene flow potentially may have played a role (table 1,  more extensive data are in supplementary  table S1). Most of these cases involve tropical taxa (21 out of 33 cases) with one additional case involving deep-water species from the tropical latitudes. Almost all cases deal with animal taxa, only one instance of ecological speciation in non-animals was found (a brown alga). Ecological speciation was found scattered across the animal kingdom (figs. 1-2). Breaking the list down further based on taxonomy, fishes are represented with eleven cases, distributed over ten families. The other instances involve six gastropods, two species complexes of shrimp, three stony corals, two octocorals, a barnacle, a sponge, a polychaete worm, an amphipod, a crab, a sea star, a bivalve, and a sea snake. The identified cases differ in the phylogenetic level on which nonallopatric divergence is thought to have occurred, from recently diverged conspecific host races (e.g., Fritts-Penniman, 2016) to older divergence among species within a genus (e.g., Jones et al., 2003).

3.1
Criterion 1: Reproductive isolation The degree of reproductive isolation varies among the taxa. In 16 of the 33 cases of potential ecological speciation with gene flow, diverged populations are described as a single species. Some degree of reproductive isolation has been reported within three of these taxa. Prezygotic (Rolán-Alvarez et al., 1999) and possibly some postzygotic reproductive isolation (Hull et al., 1996;but see Johannesson et al., 2010) has been observed between ecotypes of the intertidal gastropod Littorina saxatilis (Olivi, 1792), as well as a reduced gene flow over hybrid zones between ecotypes (Panova et al., 2006). Limited gene flow was detected between the depth-segregated ecotypes of the scleractinian coral Favia fragum (Esper, 1793)  , as well as in depth-segregated populations of the scleractinian coral Montastraea cavernosa (Linnaeus, 1767) (Serrano et al., 2014). Besides these three taxa, assortative mating between populations, for example as a result of different host-associations in philopatric species, is Downloaded from Brill.com02/25/2020 12:25:48PM via Universiteit of Groningen Table 1 Potential cases of (incipient) speciation with gene flow in marine taxa. Criteria are scored as follows: '+': yes; '-': no; '+/-': partly; '+?': potentially, possibly; '?': unknown, no data. See supplemental table S1 for a more extended assessment of criteria.  Table 1 Potential cases of (incipient) speciation with gene flow in marine taxa. Criteria are scored as follows: '+': yes; '-': no; '+/-': partly; '+?': potentially, possibly; '?': unknown, no data. See supplemental table S1 for a more extended assessment of criteria. (cont.)
cf. nucula. † In Zardus et al. (2006) classified in the genus Deminucula, which has been synonymised with the genus Nucula (Bergmans, 1978). ‡ In Castelin et al. (2012) referred to as Bursa fijiensis, which is a synonym of Bursina fijiensis (Beu et al., 2012). § In Stevens (1990)     seen in many of the taxa discussed here and leads to some of the cases of prezygotic reproductive isolation (supplementary table S1; see also below).
Rates of hybridisation can give some insight into the degree of reproductive isolation. Spawning among Caribbean hamlet species of the genus Hypoplectrus Gill, 1861 is rare, but a low rate of hybridisation has been reported (Ramon et al., 2003;García-Machado et al., 2004;Puebla et al., 2007). Walter et al. (2014) recorded a hybrid of the manta rays Mobula alfredi (Krefft, 1868) and M. birostris (Walbaum, 1792) (both species used to be classified in the genus Manta Bancroft, 1829, see White et al., 2017).

3.2
Criterion 2: Disruptive selection Evidence for disruptive selection was found in eight out of the 33 cases. This evidence mostly comes from F ST outlier tests (Lewontin & Krakauer, 1973;Lotterhos & Whitlock, 2014). Puebla et al. (2014) found one single nucleotide polymorphism (SNP), located in a region of Hox genes, to be under disruptive selection in three species of hamlets: Hypoplectrus puella (Cuvier, 1828), H. nigricans (Poey, 1852), and H. unicolor (Walbaum, 1792). Bernal et al. (2017) identified four independent loci that might have been involved in the divergence between the sympatric sister pair of grunt species Haemulon maculicauda (Gill, 1862) and H. flaviguttatum Gill, 1862, which differ in habitat use. Several studies found F ST outliers in ecotypes of the intertidal gastropod Littorina saxatilis from different localities using a wide range of methods (Westram et al., 2016 and references therein). The L. saxatilis ecotypes are a well-studied model system of speciation, in which selection by crab predation and wave exposure acts on body size and shape (Janson, 1983;Johannesson, 1986;Rolán-Alvarez et al., 1997;Johannesson et al., 2010;Le Pennec et al., 2017). These data suggest that traits on which selection acts are highly polygenic (Hollander et al., 2015;Westram et al., 2016). Finally, Carlon & Budd (2002) found strong divergence between depth-segregated morphotypes of the coral Favia fragum on an allozyme locus encoding the phosphoglucomutase enzyme, indicating potential disruptive selection as well.  found significant heritability in corallite architecture, which varies between morphotypes, confirming a genetic basis for the differences between the morphotypes.
The link between function of these genes potentially under disruptive selection and the (ecological) differences between the diverged species is not obvious. In case of hostassociated divergence, a well-known driver of ecological speciation (e.g., Drès & Mallet, 2002;Matsubayashi et al., 2010), this link may be easier to establish. Indeed, two studies did find evidence for disruptive selection on genes potentially related to survival on new host species. In the corallivorous gastropod Coralliophila violacea (Kiener, 1836), a gene involved in the control of xenobiotic detoxification pathway gene expression, which may be related to the neutralisation of host-specific metabolites, as well as four other genes involved in metal binding ions, were found to be potentially under disruptive selection (Simmonds, 2016). Functionally similar genes are thought to play a role in ecological speciation in insects (Soria-Carrasco et al., 2014;Simon et al., 2015;Simmonds, 2016). Fritts-Penniman (2016) identified a gene encoding a lysosome membrane protein in the coral-associated nudibranch Phestilla minor Rudman, 1981 as potentially being under disruptive selection, which may be involved in the treatment of host coral nematocysts in the nudibranch's digestive tract.
Results from F ST -outlier tests such as these should be interpreted carefully. Several studies reported the prevalence of false positives (Beaumont & Balding, 2004;Narum & Hess, 2011), especially when populations are isolated by distance or have undergone range expansions (Lotterhos & Whitlock, 2014), or in species living in long, linearly shaped habitats such as river, deep-sea, or coastal ecosystems Fourcade et al., 2013). Using a negative control by randomizing the data (Puebla et al., 2014) and searching for consistent outliers across geographically distant populations or across different datasets (Puebla et al., 2014;Fritts-Penniman, 2016) may reduce the rate of false positives. Including a large dataset of putatively neutral loci also increases the performance of F ST -outlier tests (Lotterhos & Whitlock, 2014).
Mechanisms other than disruptive selection may result in F ST -outliers. Loci involved in intrinsic (i.e., independent of ecology) genetic incompatibilities can easily and arbitrarily associate with locally adapted and neutral loci, especially in case of fine-grained environmental heterogeneity (e.g., in case of populations differing in host use), or when ecological selection is weak (Bierne et al., 2011). Even though this association results in a coupling between endogenous and exogenous barriers and may be responsible for the majority of F ST -outliers found in genome scans, this alternative explanation is mostly neglected in genome scans looking for signs of disruptive selection in the context of ecological speciation with gene flow (Bierne et al., 2011). The candidate loci potentially under disruptive selection discussed above (Puebla et al., 2014;Simmonds, 2016;Bernal et al., 2017) might therefore not be driving ecological speciation. Instead, the loci may be involved in or linked to genetic incompatibilites (Bierne et al., 2011;Puebla et al., 2014).
Genomic data may provide further insights in the mechanisms driving divergence between populations. In the case of the Atlantic cod Gadus morhua Linnaeus, 1758, several genes involved in adaption to low salinity, temperature, and oxygen levels have clustered together in discrete genomic regions. These regions, which likely represent large chromosomic inversions, are likely involved in local adaptation and in the divergence of fjord and offshore populations (Berg et al., 2015Kirubakaran et al., 2016;Sodeland et al., 2016;Barth et al., 2017). Genes within these regions are tightly linked and the lack of recombination may enable local adaptation in spite of high population connectivity (Barth et al., 2017).
In the other taxa, no evidence for specific genetic loci under disruptive selection has been found. However, disruptive selection on, for example, habitat choice (Levene, 1953;Maynard Smith, 1966;Gavrilets, 2006;Bolnick & Fitzpatrick, 2007) could potentially play a role in most cases. Transplant experiments on the depth-segregated ecotypes of the Caribbean octocoral Eunicea flexuosa (Lamouroux, 1821) showed for instance a decrease in survival rate of ecotypes in non-native depth ranges, suggesting strong disruptive selection (Prada & Hellberg, 2013). Similar differences in hostor habitat-associations are found in the majority of discussed taxa (supplementary table S1). The genetic mechanisms underlying these differences, as well as the strength of any disruptive selection, remain unknown.

3.3
Criterion 3 and 4: Assortative mating and the link between the selected and assortment trait Disruptive selection is not enough for ecological speciation to occur, both assortative mating and a link between the trait under disruptive selection and the assortment trait are required for reproductive isolation to evolve. In some cases, a single trait is both subjected to disruptive selection and contributing to assortative mating. The selection-recombination antagonism does not apply to these traits, which have been called 'magic traits' (or 'speciation traits') (Gavrilets, 2004;Servedio et al., 2011). While pleiotropy between selected and assortment trait was originally thought to be unlikely (Maynard Smith, 1966), such traits seem to be less rare than previously assumed (Servedio et al., 2011).
In four of the 33 cases, disruptive selection on body size may coincide with assortative mating based on body size. In seahorses, selective constraints on body size, imposed by malepregnancy, combined with size-assortative mating, as observed in seahorse populations (Vincent & Sadler, 1995;Jones et al., 2003), may have led to speciation in two pairs of sympatric species, Hippocampus erectus Perry, 1810 and H. zosterae Jordan & Gilbert, 1882 from the western Atlantic, and H. abdominalis Lesson, 1827 andH. breviceps Peters, 1869 from the Indo-West Pacific, which differ in body size (Jones et al., 2003). Similarly, body-size diverged as a result of diet specialisation in macro-and microcephalic ecotypes in species of the sea snake genus Hydrophis Latreille [in Sonnini & Latreille], 1801, and, being known as a mating cue in some sea snake species (Shine, 2005), could have induced assortative mating as well (Sanders et al., 2013b). In these two cases however, only observational data is available and disruptive selection or assortative mating could act on or be based on a correlated, unknown trait instead of body size or colour pattern (Servedio et al., 2011). Sizeassortative mating has been observed between ecotypes of the gastropod Littorina saxatilis (Johannesson et al., 1995;Hollander et al., 2005;Johannesson et al., 2008Johannesson et al., , 2010, which, as discussed above, are also subject to disruptive selection based on body size . Experimental data on assortative mating and disruptive selection is available in this case (Rolán-Alvarez, 2007;Johannesson et al., 2010), which makes the case for body size as a (single) trait driving divergence between L. saxatilis ecotypes stronger (Servedio et al., 2011). Size-assortative mating, combined with potential disruptive selection on body size driving differences in habitat use, has also been observed between the triplefin fishes Ruanoho decemdigitatus and R. whero (Wellenreuther et al., 2008).
Like body size, colour pattern has potentially been subject to both disruptive selection and assortative mating in at least one of the 33 cases discussed here. Colour-assortative mating has been observed within Caribbean hamlets (Hypoplectrus spp.) (Fischer, 1980;Barreto & McCartney, 2007;Puebla et al., 2007Puebla et al., , 2012, in which aggressive mimicry, where the fish mimic the colour pattern of non-predatory fish, may have resulted in disruptive selection on the colour pattern (Puebla et al., 2007). However, only observational data is available (Servedio et al., 2011), and it is unknown if the genomic region to be found under disruptive selection in hamlets (Puebla et al., 2014), as described above, is linked to the colour pattern of the fish. Colour patterns are a diagnostic trait in snapping shrimps of the Alpheus armatus complex as well (Knowlton & Keller, 1985;Hurt et al., 2013), and may induce assortative mating (Knowlton & Keller, 1983, 1985, as it appears to do in the basslets Gramma loreto Poey, 1868 and G. dejongi Victor & Randall, 2010(Victor & Randall, 2010Lohr et al., 2014). Colour pattern could therefore potentially have played a role in the divergence of these species, but the available data is not enough to draw any conclusions.
Body size and colour pattern can be considered 'classic' magic traits, meaning that disruptive selection is hypothesized to act on a mating cue (Servedio et al., 2011). In case of a classic magic trait, the mating cue is targeted by a preference locus elsewhere in the genome ( fig. 3a). Disruptive selection on a classic magic trait is in some cases sufficient for the evolution of reproductive isolation, such as in species displaying self-referent phenotype matching, using one's own cues as a referent for recognizing kin (Mateo, 2010), while in other cases (specifically, in theoretical, twoallele systems of speciation, see Felsenstein, 1981;Gavrilets, 2004), disruptive selection on the preference locus is required as well (Servedio et al., 2011). Whether speciation in the taxa described above follows the theoretical oneor two-allele system is not known.
In contrast to classic magic traits, 'automatic' magic traits do not require a preference locus ( fig. 3b). Disruptive selection on an automatic magic trait, in theory, automatically leads to assortative mating instead (Servedio et al., 2011). For example, divergence in the mating system of macroalgae Fucus spiralis Linnaeus, 1753, F. guiryi Zardi, Nicastro, Serrão and Pearson, 2011, and F. vesiculosus Linnaeus, 1753, as well as a difference in timing of gamete release, may have induced assortative mating and driven the divergence among these species (Ladah et al., 2008;Billard et al., 2010;Zardi et al., 2011). Species where mating occurs on a host species or within a preferred habitat (i.e., habitat fidelity) are another example where selection directly affects the assortment trait (Kirkpatrick & Ravigné, 2002;Bolnick & Fitzpatrick, 2007). A shift to a novel habitat or host species may pleiotropically bring about prezygotic reproductive isolation and could ultimately result in speciation (Drès & Mallet, 2002;Kirkpatrick & Ravigné, 2002;Bolnick & Fitzpatrick, 2007;Forbes et al., 2017). Host or habitat preference can therefore be considered an automatic magic trait (Servedio et al., 2011). Differences in habitat or host associations are found in 27 out of the cases discussed here (table 1; supplementary table S1).
While differences in habitat-associations are observed in the majority of the taxa discussed here, data on the strength of habitat or host preference is scarce, and the genetic mechanism of habitat preference remains unknown in all of the cases discussed here. Experimental evidence for host preference is only available from the four snapping shrimp species of the Alpheus armatus species complex (Knowlton & Keller, 1985), the amphipod Eogammarus confervicolus (Stimpson, 1856) (Stanhope et al., 1992), and the corallivorous gastropod Coralliophila violacea (Fujioka & Yamazato, 1983;Simmonds et al., 2018). Hostpreference was not perfect in either of these cases, suggesting that while in theory magic traits might induce full assortative mating (and therefore reduce gene flow to zero), in the real world a single trait may not be able to do so.
Habitat-preference is not required to induce (some) assortative mating in some cases. Disruptive selection over gradients of multiple environmental variables (such as temperature, oxygen concentration, and salinity) may have induced assortative mating in the Atlantic cod Gadus morhua (Berg et Neuenfeldt et al., 2013). On much smaller geographic scales, disruptive selection over a depth gradient, resulting in immigrant inviability, could in duce almost complete assortative mating between depth-segregated ecotypes of the octocoral Eunicea flexuosa, as the cumulative decrease in survival from the time of larval settling to the high age of reproductive isolation (Beiring & Lasker, 2000) results in a decrease of 'mismatched' genotypes over time, even when random settling of larvae is assumed (Pfennig, 2013;Prada & Hellberg, 2013). Similar mechanism could potentially play a role in the depth-segregated morphotypes of the octocoral Eunicella singularis (Esper, 1791) studied by Constantini et al. (2016) and the Alpheus armatus species complex (Hurt et al., 2013), but more research is required to confirm this hypothesis.
In three cases, specifically the depthsegregated morphotypes of the broadcast spawning, scleractinian corals Montastraea cavernosa and Seriatopora hystrix Dana, 1846 studied by Serrano et al. (2014) and Bongaerts et al. (2010) respectively, and the Hawaiian limpets of the genus Cellana Adams, 1869 studied by Bird et al. (2011) the presence of a mechanism inducing assortative mating is not clear.

3.4
Criterion 5: Gene flow during divergence The main challenge of inferring speciation with gene flow is the rejection of the alternative hypothesis of speciation in isolation. Specifically, distinguishing speciation with gene flow from gene flow after secondary contact is a challenging task (Gaggiotti, 2011;Strasburg & Rieseberg, 2011;Feder et al., 2013;Strasburg & Rieseberg, 2013).
Coalescence-based 'isolation-migration' (IM) models (Hey & Nielsen, 2004) are used widely to assess gene flow between populations (e.g., Won & Hey, 2005;Niemiller et al., 2008). IM models, or similar Bayesian models for inferring gene flow, were used in six out of the 33 taxa discussed here. Evidence for gene flow after divergence was found among snapping shrimps of the Alpheus armatus species complex (Hurt et al., 2013), between two species of manta rays (Mobula spp.) (Kashiwagi et al., 2012), among four closely related sea snakes of the genus Hydrophis (Sanders et al., 2013b), between depth-segregated ecotypes of the octocoral Eunicea flexuosa (see Prada & Hellberg, 2013), between depth-segregated populations of the deep-sea bivalve Nucula atacellana Schenck, 1939(Jennings et al., 2013, and between depth-segregated populations of the scleractinian coral Montastraea cavernosa (see Serrano et al., 2014). Three of these studies (Kashiwagi et al., 2012;Prada & Hellberg, 2013;Sanders et al., 2013b) used likelihood ratio tests within their IM approach and found higher support for a model including gene flow compared to models without any gene flow. IM models, as well as related methods as used by Serrano et al. (2014), cannot, however, assess heterogeneity in gene flow over time (Niemiller et al., 2010;Gaggiotti, 2011;Sousa et al., 2011;Strasburg & Rieseberg, 2011, 2013. Therefore, secondary contact after speciation without gene flow cannot be precluded based on these results (Hurt et al., 2013).
To test complex evolutionary scenarios, explicitly incorporating secondary contact, other methods are required. Approximate Bayesian computation (ABC; Beaumont et al., 2002;Hickerson et al., 2006;Beaumont, 2010;Csilléry et al., 2010) is a powerful method to test such complex hypotheses of divergence, which could be used in this context (Sousa et al., 009;Bird et al., 2012;Sousa et al., 2012;Hurt et al., 2013). Using ABC, Butlin et al. (2014) found higher support for a model of parallel, repeated divergence with gene flow between the ecotypes of the intertidal gastropod Littorina saxatilis compared to a model assuming gene flow through secondary contact.
The ABC-approach comes however with its own caveats, ABC models have been criticised by multiple authors (e.g., Templeton, 2009Templeton, , 2010Robert et al., 2011). As ABC only compares a subset of all possible models, which introduces subjectivity in the selection of models to be tested, and introduces the risk of not including the true model (Templeton, 2009; but see Beaumont et al., 2010;Csilléry et al., 2010,). Additionally, the use of insufficient summary statistics, used to compare simulation data with actual data, may lead to an unknown loss of information (Robert et al., 2011).
Data on historical gene flow is lacking in the other taxa. Ongoing gene flow, described above, is not necessarily evidence for speciation with gene flow, as, again, hybridisation can be the result of secondary contact. The ecology of species may give some clues as to whether gene flow was likely during divergence. The dispersal potential of species can strongly affect the potential degree of gene flow among populations. The potential for gene flow during divergence may be higher in species that do have a pelagic larval stage and lack genetic structuring over distance Krug, 2011). Biophysical modelling of larval dynamics, for example, showed a potential high connectivity between locally adapted populations of Atlantic cod Gadus morhua (see Barth et al., 2017). Additionally, in panmictic, sympatric sister taxa with a long, pelagic larval stage, panmixia during divergence seems a reasonable hypothesis Krug, 2011). For example, low phylogeographic structuring was found east of the Sunda Shelf in both the corallivorous gastropod Coralliophila violacea (Lin & Liu, 2008;Simmonds, 2016) and the nudibranch Phestilla minor (but see Faucci et al., 2007) (Fritts-Penniman, 2016, both of which have a long, pelagic larval stage (Demond, 1957;Taylor, 1975;Ritson-Williams et al., 2007Lin & Liu, 2008). A lack of phylogeographic structure was found in several other taxa (table 1; supplementary table S1).

3.5
Criterion 6: Geographic context of speciation While the geographic context of speciation is related to the degree of gene flow among populations, both concepts are distinct and should not be confused (Fitzpatrick et al., , 2009). The geographic distribution may give a clue about whether speciation occurred in sympatry, regardless of the evolutionary mechanisms and the degree of gene flow among diverging populations. In all but one (the Atlantic cod Gadus morhua; Barth et al., 2017; but see Neuenfeldt et al., 2013) of the cases discussed here, some overlap in geographic range between species or populations can be found. Nested distributions, where the distribution of one species is encompassed by the distribution of its sister species, are hard to explain in a framework of speciation as a result of geographic isolation. Nested distributions were found in five of the 33 taxa discussed here, namely in the Indo-Pacific coral-associated gobies of the genus Gobiodon Bleeker, 1856 (Munday et al., 2004), the seahorses Hippocampus erectus and H. zosterae (see Jones et al., 2003), the basslets Gramma loreto and G. dejongi (Victor & Randall, 2010;Lohr et al., 2014), sea snakes Hydrophis parviceps Smith, 1935, H. cyanocinctus Daudin, 1803, and H. melanocephalus Gray, 1849(Rasmussen et al., 2012Sanders et al., 2013b), three Hawaiian limpet species form the genus Cellana (Bird et al., 2011), the deep-sea gastropods Bursina fijiensis (Watson, 1881) (= Bursa fijiensis, see Beu et al., 2012) and Bursa quirihorai Beu, 1987(Castelin et al., 2012, and the eastern Pacific gobies Clevelandia ios (Jordan & Gilbert, 1882) and Eucyclogobius newberryi (Girard, 1856) (Dawson et al., 2002).
Of the taxa discussed here, overlap in distribution of pelagic fish that do not differ in habitat, such as the hamlets of the genus Hypoplectrus (see Puebla et al., 2007), is closest to 'pure' sympatry (Mallet et al., 2009). The majority of taxa discussed here differ in either habitat-or host-association, and show 'mosaic' sympatry (table 1; supplementary table  S1), where species do not coexist at a local scale, but have a patchy distribution instead (Mallet, 2008;Mallet et al., 2009). Depending on the dispersal potential of species, mosaic sympatry can be correlated with a decrease in gene flow between populations Mallet et al., 2009).

Importance of ecological speciation with gene flow in marine systems
The strengths of the arguments for ecological speciation with gene flow vary strongly among the taxa discussed here. A strong case for ecological speciation with gene flow can for example be found among the snapping shrimps of the Alpheus armatus species complex (Hurt et al., 2013) and the sea snakes of the genus Hydrophis (Shine, 2005;Sanders et al., 2013b), in which evidence for historical gene flow has been found. However, as discussed above, secondary contact cannot be excluded. Alternatives to ecological speciation with gene flow cannot be ruled out in the majority of the discussed taxa (Howell et al., 2004;Munday et al., 2004;Faucci et al., 2007;Tsang et al., 2009;Bird et al., 2011;Castelin et al., 2012). The likelihood of secondary contact compared with divergence with gene flow has only been tested in the intertidal gastropod Littorina saxatilis, in which Butlin et al. (2014) found a higher support for parallel divergence of ecotypes with gene flow using ABC. With evidence for gene flow during divergence, as well as (experimental) evidence for disruptive selection and assortative mating within ecotypes Servedio et al., 2011;Westram et al., 2016), this is the strongest case for speciation with gene flow discussed here.
Many of the other studies discussed here also provide compelling arguments for speciation with gene flow, such as those studying taxa with nested distributions (Dawson et al., 2002;Jones et al., 2003;Munday et al., 2004;Bird et al., 2011;Castelin et al., 2012;Rasmussen et al., 2012;Sanders et al., 2013b), or in which some evidence for disruptive selection has been found, such as in Caribbean hamlets (Puebla et al., 2014), grunts (Bernal et al., 2017), the corallivorous gastropod Coralliophila violacea (Simmonds, 2016), or the nudibranch Phestilla minor (Fritts-Penniman, 2016). Also, the pattern of pre-and postzygotic isolation in the sympatric Japanese greenlings Hexagrammos otakii and H. agrammus compared with their allopatric congener H. octogrammus, as described above, is strongly indicative for speciation with gene flow (Crow et al., 2010), as postzygotic isolation would be expected to be higher in case of reinforcement after secondary contact (Coyne & Orr, 1989Berlocher, 1998;Schluter, 1998;Turelli et al., 2001;Via, 2001). Speciation with gene flow seems the most parsimonious explanation in all these taxa. However, data on gene flow is not available, and is required to confirm or reject this hypothesis.
Examining the list of potential cases of marine ecological speciation with gene flow as a whole reveals a phylogenetically diverse group of taxa ( fig. 2). However, the taxonomic distribution of potential cases of divergence with gene flow seems to reflect research effort rather than actual occurrence of speciation with gene flow. There are some similarities among the taxa in which ecological speciation with gene flow may have occurred. Many species show differences in either host-or habitat associations or reproductive timing, which may act as a single trait driving divergence when species are philopatric (Servedio et al., 2011). In taxa that do not differ in habitatassociations, some other trait may act as such a 'magic' trait (e.g., body size in seahorses, see Jones et al., 2003; colour pattern in hamlets, see Puebla et al., 2007). As speciation with gene flow is thought to happen most easily in case a single trait is involved in both selection and assortative mating (Maynard Smith, 1966) and such traits are less rare than previously assumed (Servedio et al., 2011), it should be no surprise that such traits are potentially involved in most of the taxa discussed here, even though single traits may not be able to induce perfect assortative mating and decrease gene flow to zero in the real world. The Hawaiian limpet species (Cellana spp., Bird et al., 2011), as well as the depth-segregated morphotypes of the broadcast spawning, scleractinian corals Seriatopora hystrix (see Bongaerts et al., 2010) and Montastraea cavernosa (see Serrano et al., 2014) are the only taxa discussed here where a clear single trait that may potentially drive divergence cannot be indicated.
With the exception of the Atlantic cod Gadus morhua (see Barth et al., 2017), marine pelagic species are absent from the list of case studies. Several studies have reported local adaptation in pelagic species (e.g., the Atlantic herring Clupea harengus Linnaeus, 1758: Limborg et al., 2012; pteropods of the genus Cuvierina Boas, 1886: Burridge et al., 2015; the copepods Acartia tonsa Dana, 1849 andPleuromamma xiphias (Giesbrecht, 1889): Hare, 2011 andGoetze et al., 2017 respectively; the foraminiferan Neogloboquadrina pachyderma (Ehrenberg, 1861): Darling et al., 2007 Norris, 2000;Johannesson & André, 2006;Peijnenburg & Goetze, 2013;Bowen et al., 2016). Evidence for potential (historical) gene flow is, however, absent in these species. Gene flow might, therefore, be restricted for example oceanographic instead of ecological conditions. Local adaptation in pelagic species with high dispersal capabilities, even when populations occur in sympatry, does not necessarily mean that these populations diverged from a sympatric, panmictic population (Foote et al., 2011;Foote & Morin, 2016). More data on historical gene flow among locally adapted pelagic populations is needed. In G. morhua (see also supplementary table S1), some data on connectivity among populations is available. Barth et al. (2017) concluded that connectivity between locally adapted, genetically diverged Atlantic cod populations was high based on biophysical modelling of larval dynamics, these results might potentially be extrapolated to other species in the same region (e.g., the Atlantic herring C. harengus, see Limborg et al., 2012).
Identifying isolated cases of ecological speciation with gene flow, as in the majority of studies discussed here, will not progress towards understanding the frequency of ecological speciation (Coyne, 2007; see also Butlin et al., 2012). Only in a few studies discussed here has speciation been studied in a wider phylogenetic context. The importance of ecological speciation is well-studied in the hamlet genus Hypoplectrus, and is thought to have played a large role in the diversification of the genus (Ramon et al., 2003;Puebla et al., 2007Puebla et al., , 2008Holt et al., 2011;Puebla et al., 2012Puebla et al., , 2014Picq et al., 2016). For example, some potential cases of ecological speciation with gene flow have been identified based on phylogenetic studies of grunts Rocha & Bowen, 2008;Bernal et al., 2017) and seahorses (Jones et al., 2003), where ecological speciation with gene flow may have been the mechanism behind the divergence between some, but not all, pairs of sister species. The patterns from these phylogenetic studies suggest a relatively small role for ecological speciation with gene flow. Similar results were found by Reid et al. (2012) in littorinid gastropods based on phylogeographic patterns. So far, detailed, systematic studies on the importance of marine ecological speciation with gene flow in nonfish taxa are lacking. However, the studies discussed here suggest that ecological speciation with gene flow may be an important driver of divergence in coral reef symbionts (Duffy, 1996b;Munday et al., 2004;Faucci et al., 2007;Tsang et al., 2009;Duchene et al., 2013;Hurt et al., 2013;Fritts-Penniman, 2016;Simmonds, 2016;Simmonds et al., 2018) and depthsegregated populations of corals (Carlon & Budd, 2002;Bongaerts et al., 2010;Prada & Hellberg, 2013;Serrano et al., 2014;Costantini et al., 2016;Serrano et al., 2016;Bongaerts et al., 2017). Therefore, systematic assessments of many more, both species-poor and speciesrich, phylogenetically diverse clades are needed to assess the frequency of speciation with gene flow in marine systems (see also Via, 2001;Butlin et al., 2012).
Finally, an important point to consider when discussing ecological speciation with gene flow is that the degree of gene flow may not be constant during the process of speciation. Speciation is often more complex than a simple scenario of speciation either with or without gene flow (Johannesson, 2010). For example, speciation may start in isolation and be completed in the presence of gene flow, for example following secondary contact (i.e., reinforcement; Liou & Price, 1994;Butlin, 1995;Schluter, 2001;Turelli et al., 2001;Ortiz-Barrientos et al., 2009;Johannesson, 2010). A combination of different models of speciation will most likely also be applicable to most of the cases discussed here, further complicating any general conclusions that may be drawn from these results.

Directions for future research
Arguments for speciation with gene flow are lacking on several points. Taking these points into account in future studies will improve the confidence with which speciation with gene flow could be inferred. For example, the presence of single traits driving divergence is based on observational data only, with the exception of shell size in the intertidal gastropod Littorina saxatilis (Rolán-Alvarez et al., 1997;Rolán-Alvarez, 2007;Servedio et al., 2011). These 'magic' traits, both subjected to disruptive selection and resulting in assortative mating, are hypothesized to drive divergence with gene flow in the majority of cases discussed here. Data from manipulative experiments will increase the strength of evidence for a single trait driving divergence. Data on genetic architecture of adaptive traits may contribute to understanding the genomic mechanism of speciation. Chromosomal rearrangements are, for example, known to play a role in adaptation and potentially speciation (Hoffmann & Rieseberg, 2008;Schwander et al., 2014;Barth et al., 2017). As discussed above, such chromosomal rearrangements are thought to have played a role in the divergence between Atlantic cod Gadus morhua populations and may allow local adaptation despite ongoing gene flow (Berg et al., 2015Kirubakaran et al., 2016;Sodeland et al., 2016;Barth et al., 2017). Similar mechanisms might play a role in the divergence within other species discussed here. Genomic data is however still lacking in many of these non-model species, impeding the detection of such mechanisms.
Distinguishing between the population genetic process and the biogeographic pattern of speciation is important. Terms and objectives should always be defined properly to avoid the confusing of pattern with process (Bird et al., 2012). Mathematical models and statistical, population genomic methods adopt a population genetic approach, testing the degree of gene flow (i.e., divergence with gene flow) instead of the geographic context of speciation (i.e., sympatric speciation) (Gavrilets, 2003;Fitzpatrick et al., 2008). Many studies discussed here explicitly studied the geography of speciation and found evidence for sympatric speciation (e.g., Munday et al., 2004;Crow et al., 2010;Bird et al., 2011). It seems likely that sympatric speciation in most of these cases is equivalent to ecological speciation with gene flow. However, with data on (historical) gene flow lacking, definitive conclusions cannot be drawn.
The main challenge in inferring speciation with gene flow is the rejection of the alternative hypothesis of speciation without gene flow. The timing of gene flow relative to the divergence is not always considered in studies of speciation (Nosil, 2008), which may lead to incorrect assessments of speciation. To distinguish between divergence with gene flow and secondary contact, Crow et al. (2010) studied patterns of pre-and postzygotic reproductive isolation in Japanese greenlings (Hexagrammos spp.) and, as discussed above, found a pattern indicative of divergence with gene flow instead of reinforcement after secondary contact. This approach might not be generally applicable, as it may not be possible to rear certain species in laboratory environments.
Going forward, standardised, generally applicable methods, which formally test contradictory hypotheses and are able to discriminate between divergence with gene flow and secondary contact, should be used to confidently answer questions on speciation, compare results between taxa, and assess the generality of speciation with gene flow (Nosil, 2008;Johannesson, 2010). ABC may serve as a basis for such a method, even though this approach has its limitations as well, and has been criticised by some authors (e.g., Templeton, 2009Templeton, , 2010Robert et al., 2011). Several studies have used ABC to explicitly test phylogeographic models, including divergence with gene flow (e.g., Li et al., 2010;Duvaux et al., 2011;Pettengil & Moeller, 2012;Chu et al., 2013;Rougemont et al., 2016), also in marine systems (e.g., Hickerson & Meyer, 2008;Ilves et al., 2010;Boehm et al., 2013;Roux et al., 2013Roux et al., , 2014. Of the studies discussed here, only Butlin et al. (2014) used ABC, to study divergence between ecotypes of the gastropod Littorina saxatilis. Combining statistical methods such as ABC with genome-wide SNPs from nextgeneration sequencing protocols (Baird et al., 2008;Cornuet et al., 2014;Andrews et al., 2016) provides a powerful set of tools to study speciation in non-model species. Additionally, the development of third generation sequencing, providing fast and easy access to sequencing in the field and generating long reads, will likely provide new methods to study speciation in non-model species and again transform the field of evolutionary biology over the next decade (Bleidorn, 2016).
Based on all studies discussed here, we can state that speciation with gene flow has played a role in the diversification of several marine taxa, and many compelling cases of potential speciation with gene flow are to be found in the marine realm. However, the frequency of ecological speciation with gene flow cannot be assessed with these data. More systematic assessments of a more diverse range of marine taxa are required to understand the factors leading to speciation with gene flow. As speciation occurs at the population-species boundary, focussing analyses at this boundary may be most informative (Losos & Glor, 2003). Formally testing hypotheses of speciation within a wide range of marine clades (Bird et al., 2012;Butlin et al., 2012), using generally applicable, standardized methods, may bring us closer to answering the question on the importance of speciation with gene flow in marine systems.

Acknowledgements
We thank J.A.J. Metz (Leiden University), M. Schilthuizen (Naturalis Biodiversity Center, Leiden University) and an anonymous reviewer for reading an earlier version of the manuscript and providing valuable feedback. We also thank the reviewers K.T.C.A. Peijnenburg (Naturalis Biodiversity Center), K.M. Hultgren (Seattle University) and K. Johannesson (University of Gothenburg) for helpful comments on the manuscript. We also thank Arthur Anker (Federal University of Goiás), P. Larsson and K. Johannesson