Speciation results from the progressive accumulation of mutations that decrease the probability of mating between parental populations or reduce the fitness of hybrids—the so-called species barriers. The speciation genomic literature, however, is mainly a collection of case studies, each with its own approach and specificities, such that a global view of the gradual process of evolution from one to two species is currently lacking. Of primary importance is the prevalence of gene flow between diverging entities, which is central in most species concepts and has been widely discussed in recent years. Here, we explore the continuum of speciation thanks to a comparative analysis of genomic data from 61 pairs of populations/species of animals with variable levels of divergence. Gene flow between diverging gene pools is assessed under an approximate Bayesian computation (ABC) framework. We show that the intermediate "grey zone" of speciation, in which taxonomy is often controversial, spans from 0.5% to 2% of net synonymous divergence, irrespective of species life history traits or ecology. Thanks to appropriate modeling of among-locus variation in genetic drift and introgression rate, we clarify the status of the majority of ambiguous cases and uncover a number of cryptic species. Our analysis also reveals the high incidence in animals of semi-isolated species (when some but not all loci are affected by barriers to gene flow) and highlights the intrinsic difficulty, both statistical and conceptual, of delineating species in the grey zone of speciation.

Isolated populations accumulate genetic differences across their genomes as they diverge, whereas gene flow between populations counteracts divergence and tends to restore genetic homogeneity. Speciation proceeds by the accumulation at specific loci of mutations that reduce the fitness of hybrids, therefore preventing gene flow—the so-called species barriers. Importantly, species barriers are expected to act locally within the genome, leading to the prediction of a mosaic pattern of genetic differentiation between populations at intermediate levels of divergence—the genic view of speciation. At the same time, linked selection also contributes to speed up differentiation in low-recombining and gene-dense regions. We used a modelling approach that accounts for both sources of genomic heterogeneity and explored a wide continuum of genomic divergence made by 61 pairs of species/populations in animals. Our analysis provides a unifying picture of the relationship between molecular divergence and ability to exchange genes. We show that the "grey zone" of speciation—the intermediate state in which species definition is controversial—spans from 0.5% to 2% of molecular divergence, with these thresholds being independent of species life history traits and ecology. Semi-isolated species, between which alleles can be exchanged at some but not all loci, are numerous, with the earliest species barriers being detected at divergences as low as 0.075%. These results have important implications regarding taxonomy, conservation biology, and the management of biodiversity.

Multilocus analyses of the process of population divergence have been achieved in various groups of animals [ 26 , 27 ] and plants [ 28 – 30 ] for which genome-wide data are available, revealing a diversity of patterns. These case studies, however, are limited in number and have taken different approaches, such that we still lack a unifying picture of the prevalence of gene flow during early divergence between gene pools. Here, we gathered a dataset of 61 pairs of populations/species of animals occupying a wide continuum of divergence level. Species were selected in order to sample the phylogenetic and ecological diversity of animals [ 31 ], irrespective of any aspect related to population structure or speciation. We investigated the effects of genomic divergence between populations on patterns of gene flow, paying attention to the ability of ABC methods to distinguish between competing models and the influence of model assumptions.

Importantly, introgression rates alone do not govern local patterns of genetic differentiation [ 19 ]. Linked selective processes, such as hitchhiking effects [ 20 ] or background selection [ 21 ], are expected to affect the landscape of population differentiation by lowering polymorphism levels at particular loci, especially in low-recombining or gene-dense genomic regions. Neglecting this confounding effect tends to inflate the proportion of false positives in statistical tests of ongoing gene flow [ 19 ] and to mislead inferences [ 22 , 23 ]. Linked directional selection is expected to locally increase the stochasticity of allele frequency evolution, a process sometimes coined genetic draft [ 24 ]. Its effect can therefore be modeled by assuming that the effective population size, N e , which determines the strength of genetic drift, varies among loci [ 25 ].

Migration tends to homogenize allele content and frequency between diverging populations. This homogenizing effect, however, is often expected to only affect a fraction of the genome. This is because the effective migration rate is impeded in regions containing loci involved in assortative mating, hybrid fitness depression, or other mechanisms of isolation—the so-called genetic barriers [ 12 ]. Consequently, gene flow is best identified by models explicitly accounting for among-locus heterogeneity in introgression rates, as demonstrated by a number of recent studies [ 13 – 16 ]. When homogeneous introgression rate across the genome is assumed, distant lineages that have accumulated a large number of genetic barriers can be inferred as currently isolated, whereas they actually exchange alleles at a minority of loci unlinked to barriers [ 14 ]. On the other hand, neglecting heterogeneity in introgression rates between closely related lineages can result in a failure to identify some regions of the genome that are already evolving independently [ 16 , 17 ]. Heterogeneous introgression models therefore appear necessary according to the genic view of speciation [ 18 ].

As genomic data have become easier and less expensive to obtain, sophisticated computational approaches have been developed to perform historical inferences in speciation genomics (i.e., estimate the time of ancestral separation in two gene pools, changes in effective population size over evolutionary time, and the history of gene flow between the considered lineages [ 8 – 10 ]). Simulation-based approximate Bayesian computation (ABC) methods are particularly flexible and have recently attracted an increased attention in speciation genomics. One strength of ABC approaches is their ability to deal with complex, hopefully realistic models of speciation and test for the presence or absence of ongoing introgression between sister lineages. This is achieved by simulating molecular data under alternative models of speciation with or without current introgression and choosing among models based on their relative posterior probabilities [ 11 ].

Besides taxonomic aspects, the grey zone has raised an intense controversy regarding the genetic mechanisms involved in the formation of species [ 5 – 7 ]. Of particular importance is the question of gene flow between diverging lineages. How isolated must two gene pools be for speciation to begin? How long does gene flow persist as lineages diverge? Is speciation a gradual process of gene flow interruption or a succession of periods of isolation and periods of contact? These questions are not only central in the speciation literature but also relevant to the debate about species delineation, with the ability of individuals to exchange genes being at the heart of the biological concept of species.

An important issue in evolutionary biology is understanding how the continuous-time process of speciation can lead to discrete entities—species. There is usually no ambiguity about species delineation when distant lineages are compared. The continuous nature of the divergence process, however, causes endless debates about the species status of closely related lineages [ 1 ]. A number of definitions of species have thus been introduced over the 20th century, each of them using its own criteria—morphological, ecological, phylogenetic, biological, evolutionary, or genotypic. A major problem is that distinct markers do not diverge in time at the same rate [ 2 ]. For instance, in some taxa, morphological differences evolve faster than the expression of hybrid fitness depression, which in turn typically establishes long before genome-wide reciprocal monophyly [ 3 ]. In other groups, morphology is almost unchanged between lineages that show high levels of molecular divergence [ 4 ]. The erratic behavior and evolution of the various criteria is such that in a wide range of between-lineage divergence—named the grey zone of the speciation continuum—distinct species concepts do not converge to the same conclusions regarding species delineation [ 2 ].

Finally, we verified whether our inferences confirmed or contradicted the current taxonomy ( S1 table ). Our dataset comprises 26 pairs of recognized species and 35 pairs of populations (or subspecies) sharing a common binomen. Twenty-one pairs of recognized species belonged to the high-divergence zone (D a > 0.02). Of these, 16 were inferred to be currently isolated, 4 produced ambiguous results and 1 pair, Eunicella cavolinii versus E. verrucosa (gorgonian), was found to be connected by heterogeneous gene flow. Among the 5 remaining recognized pairs of species (with D a < 0.02), 2 were inferred as being fully isolated and 3 were inferred to be connected species: 2 pairs of semi-isolated species with heterogeneous gene flow (Mytilus galloprovincialis versus M. edulis and Macaca mulatta versus M. fascicularis) and the Gorilla gorilla versus G. beringei pair, which was found to be connected by homogeneous gene flow. Of the 35 pairs of recognized populations from the same species, 6 with D a > 0.02 were inferred to be isolated cryptic species. Genetic isolation has been previously suspected between northern and southern populations of Pectinaria koreni (trumpet worms) [ 38 ], between the blue and purple morphs of Cystodytes dellechiajei (colonial ascidians) [ 39 ], and between the L1 and L2 lineages of Allolobophora chlorotica (earthworms) [ 40 ], but genetic isolation is here newly revealed between Morrocan and European populations of Melitaea cinxia (Glanville fritillary), between Spanish and French populations of A. chlorotica L2, and between Mediterranean and tropical populations of Culex pipiens.

We investigated the influence of a number of ecological, geographical, phylogenetic, and life history variables on the posterior probability of ongoing gene flow. This was achieved under the heteroM_heteroN model using data from [ 31 ]. We detected no significant effect of species longevity or log-transformed propagule size (size of the developmental stage that leaves the mother and disperses) on the log-transformed probability of ongoing gene flow. In the same vein, marine organisms (n = 25) did not exhibit a higher propensity for ongoing gene flow than terrestrial ones (n = 36; r 2 below 0.01%). The log-transformed probability of ongoing gene flow was significantly higher (p-value = 0.002, r 2 = 0.14) in vertebrates (n = 20) than in invertebrates (n = 41), but the effect disappeared when the level of divergence was controlled for (net synonymous divergence <0.04: 17 vertebrate pairs, 22 invertebrate pairs, p = 0.32, r 2 = 0.03). This effect only reflects the paucity of pairs of vertebrate population/species with a high divergence in our dataset. Finally, we tested whether the current geographic distribution of species coincides with the establishment of genetic structure in our data by distinguishing pairs in which the two considered species/populations occur on the same versus distinct continents or oceans. We did not find any significant effect of this variable on the estimated values of P migration in either of the three divergence zones: D a < 0.5%, t test = –0.015269, df (degrees of freedom for the t-statistic) = 18.522, p-value = 0.988; 0.5% < D a < 2%, t test = –0.74229, df = 7.1996, p-value = 0.4814; 2% < D a , t test = 0.35512, df = 22.426, p-value = 0.7258.

We investigated the impact of assumptions about genomic heterogeneity in N e and M on the detection of current introgression ( S4 – S9 Figs). When both parameters were allowed to vary among loci, pairs of populations with D a exceeding 0.1% and showing strong statistical support for ongoing migration tended to obtain support for genomic heterogeneity in introgression rates. But when constant introgression rate was assumed (homoM_heteroN and homoM_homoN models), the importance of gene flow became underestimated in several divergent pairs of species, consistent with previous reports (e.g. [ 15 ]). When we compared models assuming homogeneous versus heterogeneous effective population size across loci, we found that the former tended to overestimate the prevalence of ongoing gene flow ( S8 Fig ), again in line with published analyses [ 19 ]. Analyses assuming homogeneous N e and M in many cases failed to support either isolation or migration, producing the highest number of ambiguous pairs ( S8 Fig ). The detected genomic heterogeneity in gene flow increased with D a until 2% of divergence. Finally, across the whole continuum, there was no significant effect of the divergence on the probability of supporting genomic heterogeneity in effective population size in our dataset.

Over the continuum of divergence, the 22 pairs with D a lower than 0.5% received a support for ongoing gene flow with a robustness ≥0.95 ( Fig 3 ). The first identified semipermeable barrier to gene flow was detected at D a ≈ 0.075%, a pair of Malurus (fairywren) species [ 37 ] for which ABC strongly supports heterogeneity in M. When the net divergence was between 0.5% and 2%, inferences about gene flow were variable and sometimes uncertain. In this grey zone, gene flow was strongly supported in 7 pairs, always with a strong support for genomic heterogeneity in introgression rates. Still, in the grey zone, ABC did not distinguish between isolation and introgression in 3 pairs of species and provided strong support for isolation in 2 other pairs. Finally, among the 27 most divergent pairs of species where D a was greater than 2%, we found 23 pairs with a strong support for current isolation and 4 ambiguous pairs ( Fig 3 ).

Each dot is for one observed pair of populations/species. x-axis: net molecular divergence D a measured at synonymous positions (log10 scale) and averaged across sequenced loci. y-axis: relative posterior probability of ongoing gene flow (i.e., SC, IM, and PAN models) estimated by ABC. Red dots: pairs with a strong support for current isolation. Grey dots: pairs with no strong statistical support for any demographic model (robustness <0.95). Blue dots: pairs with strong statistical support for genome-homogeneous ongoing gene flow. Purple dots: pairs with strong statistical support for genome-heterogeneous ongoing gene flow. Filled symbols: pairs with a strong support for genome-heterogeneous N e . Open symbols: genome-homogeneous N e . The light grey rectangle spans the range of net synonymous divergence in which both currently isolated and currently connected pairs are found (see S1 Data ).

For each of the 61 studied pairs of populations/species, we focused on synonymous positions and investigated the prevalence of ongoing gene flow by estimating the posterior probabilities of 16 different models under ABC. These 16 models represent the combinations of 5 demographic models (SI, AM, IM, SC, and panmixia) and four assumptions regarding the genomic heterogeneity in introgression (for AM, IM, and SC only) and drift rates (for all models; see above and Material and Methods ). The posterior probability P migration that the two populations currently exchange migrants was estimated by summing the contributions of the PAN, IM, and SC models ( Fig 1 ) and plotted against measures of molecular divergence ( Fig 3 ). D a , which can be understood as the per-site amount of neutral derived mutations being fixed in the different lineages, provided the best relationship ( Fig 3 ). Results with other measures of divergence and with the estimated age of the split (T split parameter under the IM model) are also shown ( S4 – S7 Figs).

We computed various measures of molecular divergence between species/populations: namely, D a , the relative average divergence, corrected for within-species diversity [ 36 ]; D xy , the absolute average divergence; and F ST , a classical measure of population differentiation. In our dataset, D a ranged from 5.10 −5 (French versus Danish populations of Ostrea edulis) to 0.309 (Crepidula fornicata versus C. plana) and F ST from 0 (between Anas crecca shemya and A. crecca attu) to 0.95 (between Camponotus ligniperdus and C. aethiops, S3 Fig ). As expected, D a was strongly correlated to F ST and less well to the absolute divergence D xy ( S3B Fig ). The across-loci variance in F ST was minimal for low and high values of D a ( S3B Fig ), which reflects an F ST homogeneously low at early stages of divergence, homogeneously high at late stages of divergence, and heterogeneous among genes at intermediate levels of D a ( S3 Fig ).

The posterior probability of ongoing gene flow was estimated in 61 pairs of species/populations of animals ( S1 Data ) showing variable levels of molecular divergence ( S1 Data ). Fifty pairs were taken from a recent transcriptome-based population genomic study [ 31 ], with two individuals per population/species being analyzed here. The datasets for the other 11 species pairs were downloaded from the NCBI ( S1 Data ). They correspond to sequences from published studies using either ABC, Ima [ 33 ], or MIMAR [ 35 ], for which 3 to 78 diploid individuals were analyzed.

In addition, the robustness of the ABC inference was only weakly dependent on the sample size when the number of loci was greater than 100: similar results were obtained when we simulated samples of size 2, 3, 25, or 50 diploid individuals ( S1 Fig ). Finally, and importantly, simulations showed that ABC is not accurate enough to discriminate between the IM and SC models. Datasets simulated under SC were assigned to SC with high confidence only when the period of isolation before secondary contact represents at least a proportion of about 60% of the total divergence time ( S2A Fig ). When shorter periods of isolation were simulated, the method either assigned the datasets to IM or did not provide an elevated posterior probability to any demographic model ( S2B Fig ).

Posterior probability P migration to support ongoing migration was estimated for a total of 116,000 simulated datasets across 16 models. A. P migration as a function of the net synonymous divergence D a . Dots represent datasets simulated under the IM, SC, and PAN models. The colors show datasets for which gene flow is correctly supported (green) or wrongly rejected (red). Grey dots represent datasets for which the robustness of the ABC analysis is <0.95. B. P migration as a function of the net synonymous divergence D a . Dots represent datasets simulated under the SI or AM models. The colors show datasets for which gene flow is correctly rejected (green; robustness ≥ 0.95) or wrongly supported (red; robustness ≥ 0.95). C. Proportion of true positives (green), false positives (red), and ambiguous analyses (grey) for different ranges of D a across IM, SC, and PAN datasets. Horizontal red line shows 5%. D. Proportion of true positives (green), false positives (red), and ambiguous analyses (grey) for different ranges of D a across SI and AM datasets.

Among the 58,000 simulated datasets in which current gene flow was assumed (IM, SC, and PAN; Fig 2A ), 99.462% were true positives (P migration > P isolation and robustness ≥ 0.95), 0.129% were false positives (P migration < P isolation and robustness ≥ 0.95), and 0.409% were ambiguous cases for which ABC did not provide any robust conclusion (robustness < 0.95). Among the 58,000 simulated datasets in which current isolation was assumed (SI and AM; Fig 2B ), 99.649% were true positives (P isolation > P migration and robustness ≥ 0.95), 0.002% were false positives (P isolation < P migration and robustness ≥ 0.95), and 0.34% were ambiguous cases (robustness < 0.95). When current gene flow was assumed, the rates of false positive and ambiguity were very low at every level of population divergence. When current isolation was assumed, a higher rate of ambiguity, but no elevation of the rate of false inference, was observed at low levels of divergence (D a < 0.01, Fig 2D ). This contrasts with the recent suggestion that the full-likelihood method developed in the IMa2 software [ 33 ] might be biased towards supporting current gene flow when isolation is recent [ 19 , 34 ]—our approach appears to be immune from this bias. To specifically address this point, we repeated the exact same simulations as in [ 34 ] and confirmed that our ABC approach has a reduced power (i.e., more ambiguous cases with robustness <0.95) when the split is recent but still a very low rate of false positive in these conditions (see S1 text ).

We explicitly tested the hypothesis of current gene flow by comparing the relative posterior probabilities of 16 models for 61 pairs of species distributed along a continuum of molecular divergence. In the ABC framework, the posterior probability of a model corresponds to its relative ability to theoretically produce datasets similar to the observed dataset, compared to a set of alternative models. Before analyzing datasets from the 61 pairs of animal species, we first assessed the power of the adopted ABC approach to correctly distinguish between models involving current isolation (SI + AM) versus ongoing migration (IM + SC + PAN). This was achieved by randomly simulating 116,000 datasets distributed over the 16 compared models and applying our ABC inference method to each of them. Specifically, we investigated which model had the highest posterior probability and assessed significance by estimating the associated robustness—the probability to correctly support a model given its posterior probability. A robustness greater than 0.95 can be interpreted as a p-value below 0.05 [ 32 ]. The analysis of simulated datasets allowed us to empirically measure a threshold value of 0.6419 for the posterior probability P migration (= P IM + P SC + P PAN ), above which the robustness to support ongoing migration is greater than 0.95. Similarly, a posterior probability P migration below 0.1304 implied a statistical support for the current isolation model with a robustness greater than 0.95.

The dominant assumption in published demographic inferences is the homoN submodel, in which it is assumed that most of the genetic variation in the genome is unaffected (or equally affected) by selection at linked sites. Here, homoN was simulated using a single value of effective population size shared by all loci across the genome, but the effective population size differed among populations. The heteroN submodel accounts for local genomic effects of directional selection (background selection, selective sweeps) by considering a variable effective population size among loci, here assumed to follow a rescaled beta distribution. The homoM submodel assumes that all loci share the same probability to receive alleles from the sister population (i.e., posits the absence of species barriers or of adaptively introgressed loci). Alternatively, the heteroM submodel accounts for the existence of local barriers to gene flow, of variable strengths, and of variable levels of genetic linkage to the sampled loci. HeteroM was here simulated by assuming that the effective introgression rate is beta distributed across the genome, thus intending to account for the combined effects of selection, recombination, and gene density. In principle, one could explicitly include information on local recombination rates and gene density, but no such data was available in the species analyzed here.

SI = strict isolation: subdivision of an ancestral diploid panmictic population (of size N anc ) in two diploid populations (of constant sizes N pop1 and N pop2 ) at time T split . AM = ancestral migration: the two newly formed populations continue to exchange alleles until time T AM . IM = isolation with migration: the two daughter populations continuously exchange alleles until present time. SC = secondary contact: the daughter populations first evolve in isolation (forward in time), then experience a secondary contact and start exchanging alleles at time T SC . PAN: panmictic model. All individuals are sampled from the same panmictic population. Red phylogenies represent possible gene trees under each alternative model.

Five demographic models differing by the history of gene flow between two diverging populations were considered ( Fig 1 ), namely strict isolation (SI), ancient migration (AM), isolation with migration (IM), secondary contact (SC), and panmixia (PAN). The latter three models involve ongoing gene flow between the two populations, whereas the former two do not. The five demographic models were subdivided into different genomic submodels that reflect alternative assumptions about the genomic distribution of indirect selective effects on the effective population size (homoN if homogeneous or heteroN if heterogeneous) and on the migration rate (homoM if homogeneous or heteroM if heterogeneous). Heterogeneous effective population size was considered in all the models, while heterogeneous migration rate was considered in models with gene flow (IM, AM, and SC). The SI and PAN models were divided into two submodels (homoN and heteroN), and the AM, IM, and SC models were divided into four submodels (homoN_homoM, homoN_heteroM, heteroN_homoM, and heteroN_heteroM).

Discussion

We performed a comparative speciation genomics analysis in 61 pairs of populations/species from various phyla of animals. Our ABC analysis, which takes into account the confounding effect of linked selection heterogeneity, provides a first global picture of the prevalence of gene flow between diverging gene pools during the transition from one to two species.

Accounting for Among-Locus Heterogeneity in Drift and Migration Rate Inferring the history of divergence and gene flow, which determines the rate of accumulation of species barriers, is of prime importance to understand the process of speciation [17]. This can be achieved by various methods, among which ABC approaches have proven particularly flexible and helpful to compare alternative evolutionary models. Our analysis of simulated datasets illustrates that ABC methods have the power to effectively discriminate recent introgression versus current isolation based on datasets of several hundreds of loci and a few individuals per species—typical of population genomic studies. Comparisons of alternative demographic models, however, can be strongly impacted by assumptions regarding the genomic distribution of effective population size (N e ) and introgression rate (M). Heterogeneities in N e and M are common in natural populations as a result of selective processes applying either globally (background selection [19,41,42]) or specifically against migrants (genetic barriers [12,43]). Following [13], we here introduced a framework in which each of the two effects, or both, can be readily accounted for. In our analysis, the number of pairs of populations/species for which ambiguous conclusions were reached was maximal when genomic heterogeneities of both migration and drift were neglected. Incorporating within-genome variation in N e tended to enhance the support for models with current isolation, as previously suggested [19]. The heteroN model makes a difference regarding inference of current gene flow between the highly divergent Ciona intestinalis and C. robusta species (see below). Conversely, incorporating heterogeneity in M doubled the number of pairs for which ongoing gene flow was supported when compared to analyses with homogenous M, in which most of these pairs exhibited ambiguous results. Our study therefore underlines the importance of accounting for genomic heterogeneities for both N e and M when comparing alternative models of speciation [14,15,19] and calls for prudence regarding the conclusions to be drawn from the analysis of a single pair. However, it is important to recall here that the action of natural selection on its molecular target and neighborhood is more complex than a simple reduction in N e . Our modeling of genomic heterogeneity in drift and selection by a beta distribution of N e throughout the genome is an approximation which cannot replace an explicit modeling of these processes. In our modeling, we assumed that a given locus i is independently affected by drift and selection in all of the simulated populations including the ancestral one. Our choice was motivated by the generality of this model. An alternative approach to model genomic heterogeneity in N e can be to assume that background selection is the main process shaping genomic landscapes of diversity. This can be approximated by assuming that a locus i is equally affected by drift and selection in all populations instead of assuming independent effects as in our study. Among models assuming ongoing gene flow, our ABC analysis of simulated and empirical data often failed to discriminate between the isolation with migration and secondary contact models. These two models yield similar signatures in genetic data, such that only relatively recent secondary contacts following long periods of interrupted gene flow can be detected with high confidence (S2D Fig) [44]. Similarly, among models excluding ongoing gene flow, distinguishing between strict isolation and ancient migration was not possible in a substantial number of cases. These are challenges for future methodological research in the field, with important implications regarding the debate about the requirement of geographic isolation to complete speciation [7,45]. Only two diploid individuals per population/species were used in this analysis for the sake of comparability between datasets (in many populations, no more than two individuals were available) and because of computational limitation. However, our evaluation of the effect of sample size on ABC-based demographic inference suggested that two diploid individuals per population were largely sufficient to capture the main signal when more than 100 loci are available (S1 Fig).

Prevalent Gene Flow between Slightly Diverged Gene Pools Although ABC analyses of particular pairs of populations can be affected by the choice of model of genomic heterogeneity, the overall relationship between net molecular divergence and detected ongoing gene flow was qualitatively similar among analyses. Pairs of populations diverging by less than 0.5% were found to currently exchange migrants. This includes populations that form a single panmictic gene pool and pairs of diverging populations/species connected by gene flow. The low-divergence area contains pairs of populations showing conspicuous morphological differences, such as eastern versus western gorilla or the cuniculus and algirus subspecies of rabbit (Oryctolagus cuniculus). No pair of populations in this range of divergence was supported to be genetically isolated or yielded ambiguous results. Simulations indicate that our ABC approach is not expected to yield false inference of gene flow in recently isolated populations, contrary to what was suggested with the full-likelihood approach of IMa2 [34]. The main risk is rather a false inference of isolation despite gene flow (Fig 2), which can be explained by the fact that the SI model is less parameterized than models assuming gene flow (IM and SC). ABC had a low false positive rate even when we simulated very recent splits, as has been done in previous papers [19,34]. This is probably because in strict isolation, shared polymorphisms are quickly sorted into private polymorphisms and fixed differences after population split, such that D a can hardly be very small in the absence of gene flow [46]. Our analysis therefore identifies D a < 0.5% as a good synthetic proxy to attest for the existence of gene flow. Other measures of divergence, although producing a qualitatively similar pattern, did not predict the existence of current gene flow as nicely as D a did. Pairs in the low range of divergence must correspond to populations that did not accumulate sufficiently strong and numerous genetic barriers, such that gene flow currently occurs at important rates. The detection of significantly heterogeneous introgression rates in a number of low-diverged pairs (D a < 0.5%) demonstrates the ability of our ABC approach to detect semipermeable barriers quite efficiently at early stages of speciation and supports the rapid evolution of Dobzhansky–Muller incompatibilities [47,48]. A majority of the pairs from the low-divergence area, however, did not yield any evidence for among-locus heterogeneity of introgression rate. Some might correspond to effectively isolated backgrounds that are missed by our method by lack of power when the signal of heterogeneity is too tenuous. It is quite plausible, however, that some pairs of populations/species in the low-divergence zone have differentially fixed mutations with major effects on hybrid fitness, whereas others do not because of mutational stochasticity and/or across-taxa differences in the genetic architecture of barriers—i.e., simple (two locus) versus complex incompatibilities and strength of associated selective effects [49].

Suppressed Gene Flow at High Sequence Divergence At the other end of the continuum, it appears that above a divergence of a few percent, barriers are strong enough to completely suppress gene flow: almost all pairs of species with D a > 2% were found to have reached reproductive isolation with strong support. This might result from impaired homologous recombination because of improper pairing of dissimilar homologous chromosomes at meiosis, which would reduce the fecundity of hybrids [50,51]. Of note, the upper threshold for reproductive isolation (D a = 2%, D x y = 5.5%) is of the order of magnitude of the maximal level of within-species genetic diversity reported in animals [31,52], somewhat consistent with the hypothesis of a physical constraint imposed by sequence divergence on the ability to reproduce sexually. Alternatively, the 2% figure may represent a threshold above which Dobzhansky–Muller incompatibilities are normally in sufficient number and strength to suppress introgression. The two hypotheses are not mutually exclusive but pertain to distinctive processes of genetic isolation; the former would be maximally expressed during F1 hybrid meiosis, while the latter would affect recombined, mosaic individuals carrying alleles from the two gene pools at a homozygous state. In the high-divergence area, no instance of among-locus heterogeneous migration was detected, indicating that introgression is blocked across the whole genome in these pairs of species. A number of highly divergent species pairs yielded support for among-locus heterogeneous N e , suggesting that the same regions of the genome are under strong background selection in the two diverging entities—presumably regions of reduced recombination and/or high density in functional elements. Neglecting the genomic heterogeneity in N e can lead to false inference of gene flow. For instance, allowing genomic heterogeneity in M but not in N e led to strong statistical support for a secondary contact between the highly divergent Ciona intestinalis (formerly C. intestinalis B) and C. robusta (formerly C. intestinalis A) species (S4 and S5 Figs), consistent with [14], but accounting for heterogeneity in both M and N e resulted in an ambiguous result without a sufficiently strong support for any models. The among-locus variance in differentiation between these two species, which was interpreted as mainly reflecting introgression at a few loci in [14], is shown here to possibly be the result of a more complex situation that our models failed to capture.

Intermediate Divergence Levels: The Grey Zone of Speciation The area of intermediate divergence from 0.5% to 2% of net synonymous divergence unveils the grey zone of the speciation continuum. In this grey zone, isolated pairs of populations/species coexist with pairs connected by migration, and the latter are mainly composed of semi-isolated genetic backgrounds, the situation under which taxonomic conundrums flourish. Cases of ambiguous conclusions about the demographic history also tended to be found in this intermediate zone, perhaps reflecting instances of complex divergence models that are not well predicted by our demographic models. Researchers should be ready to face problems regarding demographic inference—and therefore parameter estimation—when conducting a project of speciation genomics falling in the grey zone. Accounting for genomic heterogeneity of introgression and drift rates appears to be crucial for detecting current gene flow in this range of divergence (S4–S7 Figs). For instance, the mussel species M. galloprovincialis versus M. edulis and the gorgonian species Eunicella cavolinii versus E. verrucosa are the two most divergent pairs for which ongoing introgression was detected, but this only appeared when the genomic variation in M was accounted for—the homoM_homoN and homoM_heteroN models yielded ambiguous conclusions about these pairs of species, in which the existence of semipermeable barriers has previously been demonstrated [53,54]. Our analysis revealed significant among-locus heterogeneous migration in as many as thirteen pairs of populations/species (Fig 3). This illustrates the commonness of semipermeable genomes at intermediate levels of speciation, when some, but not all, genomic regions are affected by barriers to gene flow. Besides mussels and gorgonians, heterogeneous gene flow was newly detected between American and European populations of Armadillidium vulgare (wood lice) and Artemia franciscana (brine shrimp), between Atlantic and Mediterranean populations of Sepia officinalis (cuttlefish), and between the closely related Eudyptes chrysolophus moseleyi versus E. c. filholi (penguins) and Macaca mulatta versus M. fascicularis (macaques)—in addition to the previously documented mouse [55], rabbit [56], and fairywren [57] cases. The grey zone, finally, includes populations between which unsuspected genetic isolation was here revealed, such as the Moroccan versus European populations of Melitaea cinxia (Glanville fritillary) and the Spanish versus French populations of A. chlorotica L2 (earthworm), which according to our analysis correspond to cryptic species. Our genome-wide approach and proper modeling of heterogeneous processes therefore clarified the status of a number of pairs from the grey zone, emphasizing the variety of situations and the conceptual difficulty with species delineation in this range of divergence.