Christelle Fraïsse

Cl Ément Mazoyer

Jonathan Romiguier

Étienne Loire

Alexis Simon

Nicolas Galtier

Laurent Duret

Nicolas Bierne

Xavier Vekemans

Camille Roux
email: camille.roux@univ-lille.fr

DILS : Demographic Inferences with Linked Selection

e present DILS, an online statistical analysis platform for conducting demographic inferences with linked selection from population genomic data using an Approximate Bayesian Computation framework. DILS takes as input single-population or two-population datasets and performs three types of analyses in a hierarchical manner, identifying: 1) the best demographic model to study the importance of gene flow and population size change on the genetic patterns of polymorphism and divergence, 2) the best genomic model to determine whether the effective size N e and migration rate N.m are heterogeneously distributed along the genome and 3) loci in genomic regions most associated with barriers to gene flow. Available via a web interface, an objective of DILS is to facilitate collaborative research in speciation genomics. Here, we show the performance and limitations of DILS by using simulations, and finally apply the method to published data on a divergence continuum composed by 28 pairs of M ytilus mussel populations/species.

Introduction

Population genomic data along with efficient computational methods are becoming increasingly available, paving the way to broad-scale application of modelbased inferences for understanding signatures of evolutionary processes [START_REF] Lohse | Come on feel the noise-from metaphors to null models[END_REF]. Neutral processes such as divergence, gene flow and changes in population size all shape patterns of genomic variation; and so demographic models attempting to reconstruct the past history of single populations or closely-related species can also serve as null models in genome scans for selection. Considering a single species, model-based inferences are especially suitable in domesticated crops for disentangling the effect of population size changes from selection on agronomic traits [START_REF] Gaut | Demography and its effects on genomic variation in crop domestication[END_REF]. Two-population models allow to tackle issues on speciation genomics, where this approach provides direct testing of distinct modes of speciation [START_REF] Sousa | Understanding the origin of species with genome-scale data: modelling gene flow[END_REF], with at the two extremes a model of allopatric speciation that occurs in complete isolation and a model where speciation is opposed by continuous gene flow. This is critical to build-up a unifying picture of the genic view of speciation by quantifying the reduction in gene exchange between lineages as a function of their molecular divergence [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF][START_REF] Peñalba | Current geography masks dynamic history of gene flow during speciation in northern Australian birds[END_REF]; and identify in silico genomic regions harboring speciation genes, given that their barrier effects can only be detected in the presence of ongoing gene flow [START_REF] Roux | Crossing the species barrier: genomic hotspots of introgression between two highly divergent Ciona intestinalis species[END_REF]. At a broader scale, model-based inferences can be applied to community ecology to infer, for example, the assembly history of trophically linked species [START_REF] Bunnefeld | Whole-genome data reveal the complex history of a diverse ecological community[END_REF].

Various methods have been proposed to extract such information from population genomic data. Site frequency spectrum (SFS)-based methods compute or approximate the likelihood of the allele frequency distribution from a demographic model using either the diffusion approximation [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF], the moment closure [START_REF] Jouganous | Inferring the joint demographic history of multiple populations: beyond the diffusion approximation[END_REF] or the coalescent [START_REF] Excoffier | Robust demographic inference from genomic and SNP data[END_REF]. While these methods are fast, they ignore linkage information which is informative about past demography [START_REF] Terhorst | Fundamental limits on the accuracy of demographic inference based on the sample frequency spectrum[END_REF]. Therefore, other methods rely on the block-wise SFS, i.e. the SFS of short non-recombining blocks of sequences [START_REF] Lohse | A general method for calculating likelihoods under the coalescent process[END_REF]. That way the genealogical information contained within each block is combined along the genome. Other multilocus methods can jointly infer recombination and demography, therefore capturing longer range linkage disequilibrium, but they are still restricted to simple demographic histories excluding migration [START_REF] Terhorst | Robust and scalable inference of population history from hundreds of unphased whole genomes[END_REF]. Still, the flexibility of simulation-based approximate Bayesian computation (ABC) enables including recombination within un-linked blocks in multi-locus inference of complex (and hopefully more realistic) evolutionary scenarios [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF]. Although more computationally expensive, the analysis of thousands of loci results in high-precision parameter estimation for most demographic scenarios [START_REF] Robinson | ABC inference of multi-population divergence with admixture from unphased population genomic data[END_REF][START_REF] Smith | Leveraging whole genome sequencing data for demographic inference with approximate Bayesian computation[END_REF].

In this paper, we present an ABC framework (DILS) building upon and extending current statistical machinery [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF][START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. Our method is flexible both in terms of the evolutionary scenarios that can be accommodated (allowing changes in population size over time, linked selection and implementing various models of migration), and type of data (SFSs and/or multi-locus sequences). A major improvement compared to most existing methods is decoupling the effect of linked selection and neutral history by relaxing the assumption that all loci share the same demography (see [START_REF] Sethuraman | Model-based assessments of differential introgression and linked natural selection during divergence and speciation[END_REF]Sousa et al., 2013 and Lohse pers. comm. for similar ideas). We model variation in the rate of drift among loci to account for linked selection effects due to background selection (i.e. purifying selection) and adaptive sweeps in low-recombining and genedense regions. And by explicitly modelling variation in migration rates among loci in two-population models, we can capture the effect of selection against migrants at neutral markers linked to species barriers, and so analyse further these candidate genomic regions for reproductive isolation [START_REF] Roux | Crossing the species barrier: genomic hotspots of introgression between two highly divergent Ciona intestinalis species[END_REF].

DILS offers an online platform for configuring demographic inferences based on genomic data of thousands of loci, performing them and visualizing the returned output. These advances are made possible by progress in simulator performance [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF], reduction in the number of simulations required to train prediction algorithms [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF] and development of computer clusters and tools facilitating parallelism [START_REF] Köster | Snakemake-a scalable bioinformatics workflow engine[END_REF]. Following other userfriendly ABC programs, DILS aims to ease the use of high-performance tools for non-experts in methodology [START_REF] Cornuet | Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation[END_REF][START_REF] Cornuet | DIYABC v2. 0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data[END_REF]. Importantly, as there is a limit to how much information can be extracted from genomic data, DILS also implements rigorous quality controls. Therefore, not only does the user receive 1) the best-supported model among those proposed (figure 1), 2) an estimate of the demographic parameters describing this model and 3) a locus-specific test to identify barriers to gene flow (when relevant); the user will also get feedback on whether the best model is relevant and to which extent the estimates are able to reproduce the observed data.

A long-term aim of DILS is to facilitate collaborative research in speciation genomics. The degree of reproductive isolation appears to follow a quasi-shared molecular clock among animals, depending on the level of net genomic divergence between lineages [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. However, for the same level of divergence, two opposite situations coexist in the so-called grey zone of speciation with, on the one hand, semi-isolated pairs capable of exchanging genetic material and, on the other hand, pairs of species that are fully reproductively isolated. Many hypotheses have been advanced to explain such a reproductive barrier contrast within the same range of molecular divergence, including differences related to life history traits (internal versus external fertilization), ecology (marine versus terrestrial organisms), reproductive systems (e.g. in plants: self-incompatibility versus self-fertilization), genome size and recombination landscape, functional redundancies in genomes, etc... Speciation is such a multi-factorial process that it seems impossible for a single research group to study these different components. Consequently, our aim is to include in DILS a collaborative science option, allowing to feed in real time the relationship between molecular divergence and genetic isolation between lineages. Available as a choice, the sharing of the inferences made in DILS, associated with the expertise that users have about their biological model, will contribute to a long-term collaborative study aiming to better understand the speciation process. This objective is illustrated here with the analysis of 28 new pairs of mussel populations whose transcriptomes were recently published, revealing ongoing gene flow for levels of divergence greater than 2%.

In this study, we have four objectives: 1) providing a flexible and powerful demographic inference method with linked selection to analyse genome-scaled dataset in single and two-population models; 2) presenting a user-friendly web-platform that implements this approach and paves the way for collaborative science; 3) testing the performance and limitations of the method by using simulations; 4) and applying it to an empirical dataset of M ytilus mussels.

Results

In the current version of DILS, evolutionary scenarios can be investigated for sampling schemes involving one or two populations. For both types of analysis, the first step in DILS is to compare the demographic models described in figure 1. With a single population, DILS will examine the changes in size over time. With two populations, such variations in population size are also implemented, but DILS will additionally compare alternative temporal patterns of introgression.

An innovative feature of DILS is to include linked selection, either through the effect of background selection (and selective sweeps) that modulates the effective population size along the genome, or through the effect of selection against migrants that reduces locally the effective introgression rate in genomic regions locked to gene flow. Therefore, all demographic models exist under two alternative genomic models regarding the effective population size (homo-Ne versus hetero-Ne), and, in models with migration, the introgression rate (homo-N.m versus hetero-N.m), depending on whether these parameters are homogeneous or heterogeneous among loci.

Comparisons of demographic and genomic models

In this section, we present DILS performance to compare demo-genomic models involving one (section 2.1.1) or two (section 2.1.2) populations. These evaluations were performed by analyzing pseudo-observed datasets simulated under specified models, in order to assess the efficacy of our approach to correctly support the true model. In both analyses, a given demographic model corresponds to the set of all its genomic submodels. All model comparisons are performed using random forest algorithms [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF][START_REF] Fraimout | Deciphering the routes of invasion of Drosophila suzukii by means of ABC random forest[END_REF].

Single population models

For studies where a single population is sampled, three demographic models are compared describing either i) a Constant population size N e current , ii) recent demographic Expansion or iii) Contraction. Demographic changes are assumed to be instantaneous, with a transition from N e past to N e current occurring T dem generations ago (figure 1-A).

For a given dataset, we first estimate the best-fitting demographic model among those depicted in figure 1 by carrying out 10, 000 simulations under each genomic alternative sub-model (homogeneous versus heterogeneous Ne). These simulations produce reference tables, i.e., a set of simulated summary statistics used to train a random forest algorithm to predict which of the proposed models best explains observed data. Thus, in the comparison Expansion versus Constant versus Contraction (algorithm 1), each model was simulated 10, 000 times under the homogeneous Ne model, and 10, 000 times under the heterogeneous Ne model.

To test the performance of DILS in model comparisons (among demographic models and among genomic models), we simulated 30, 000 pseudo-observed datasets of 100 loci by using random combinations of parameters. Throughout the manuscript, pseudoobserved datasets are used to evaluate the performance of the random forest and do not contribute to its training. These simulated datasets are equally distributed between the demographic models (10, 000 for each of the three Expansion/Constant/Contraction models) and between the genomic models (5, 000 for each of the two homo/hetero Ne alternative genomic models for a given demographic model). Then, for each of these pseudo-observed datasets, we apply step 1 and step 2 of the algorithm 1 in order to obtain for each model M, the proportion of pseudo-observed datasets that is i) correctly and strongly captured by the random forest approach, ii)falsely and strongly captured and Data: A single fasta containing all genes sequenced in all individuals sampled in the studied population Result: Posterior probabilities for the best 1) demographic and 2) genomic models Conversely, the Isolation with Migration (IM) and Secondary Contact (SC) models describe two populations that are currently connected by introgression events at rate N.m. The two-population models shown here are of constant size, but DILS optionally incorporates alternative versions of the same four models where effective size can change independently in both daughter populations between the present time and T split .

iii) ambiguous, i.e, associated to an insufficiently high posterior probability (figure 2).

Figure 2: Performance of DILS for hierarchical comparison of single-population models

The required data is a single fasta file containing the sequences obtained for different genes, in different individuals. DILS will first perform a comparison of the three demographic models (Expansion versus Constant versus Contraction). In a second step, DILS will compare two genomic models (homogeneous versus heterogeneous genomic distribution of N e) for the best supported demographic model. The pie charts designate for each model the proportion of simulations performed under the corresponding model that is strongly and correctly captured (correct: blue), strongly and incorrectly captured (wrong: yellow) and without strong statistical support for any of the studied models (ambiguous: purple). The performance of DILS for each comparison was based on 10, 000 pseudo-observed datasets for each of the Expansion/Constant/Contraction demographic models. Each of these 10,000 simulated datasets are evenly distributed between the two genomic models, homo and hetero Ne. The parameters used for the simulated datasets are randomly drawn from uniform laws, with N e in [1-1,000,000] individuals and T dem in [1-2,000,000] generations. Each simulated dataset consists of 100 loci of length 1, 000 nucleotides.

For a given dataset, the model comparison will provide two pieces of information: i) What is the best model among all those arbitrarily proposed in the analysis? ii) With which posterior probability is the best model supported? It is from the latter probability that we determine whether an inference is strong or ambiguous. For this we use a probability threshold beyond which an inference is considered strong. This threshold is determined recursively on the basis of the false positive rate which decreases monotonically by increasing the value of the threshold. From datasets randomly simulated under different models, we establish a threshold value such that the false positive rate is less than or equal to 1%. With this approach, the falsepositive rate remains consistently low, but the relative proportions of true positives versus ambiguous cases vary according to the power of the ABC to discriminate among an arbitrary set of models. Thus, among the 30, 000 pseudo-observed datasets simulated under the Expansion, Constant and Contraction models, if a threshold is applied that keeps the error rate below 1%, the proportions that are correctly supported by our approach are 86%, 89% and 99% respectively, while the proportions that are ambiguous are 13%, 10% and 0.2% (figure 2).

The second step of the hierarchical comparison, which classifies models with genomic variation in effective size (hetero) or without variation (homo), is evaluated using the same procedure as step 1 (figure 2). For this purpose, the proportion of ambiguous, correct and false inferences are measured among 5, 000 pseudo-observed datasets simulated for each of the 6 combinations [demographic models] x [genomic models]. For the Expansion and Constant demographic models, the correct recapture rates of homo and hetero Ne genomic models range from 75% (Expansion + hetero Ne; 24% of ambiguity) to 86% (Constant + homo Ne; 13% of ambiguity). Finally, while recovering the Contraction demographic model is a very robust analysis with 99% of inferences that are both correct and associated with a high posterior probability, it is more complicated to distinguish "Contraction + homo Ne" from "Contraction + hetero Ne". About 41% of the pseudo-observed datasets simulated in the "Contraction + homo Ne" model are correctly captured by the random forest, and only 26% for the "Contraction + hetero Ne" model. The occurrence of a recent bottleneck tends to reduce the genomic variance of Ne to levels that generate apparent homogeneity.

Two population models

The two-population models are grouped into two supermodels: with current isolation (Strict Isolation + Ancient Migration; see figure 1) and ongoing migration (Isolation Migration + Secondary Contact). The first step of the hierarchical comparisons performed by DILS therefore aims to determine which supermodel best explains the data observed in the two sampled populations (see algorithm 2 and figure 3). This is achieved by labeling as "isolation" all reference simulations performed under the two SI models (with homo-Ne or hetero-Ne) and the four AM models (homo-Ne or hetero-Ne, and homo-N.m or hetero-N.m). All other models are labeled as "migration" supermodel.

To test the power of DILS to correctly recapture an isolation or migration model, the same simulationbased evaluation as in section 2.1.1 is performed here. We evaluated the performance of DILS for 60, 000 pseudo-observed datasets simulated under the "isolation" supermodel (10, 000 for each combination of [SI; AM], [homo-Ne; hetero-Ne] and [homo-N.m; hetero-N.m] for the AM model only) and 80, 000 under the "migration" supermodel (10, 000 for each combination of [IM; SC], [homo-Ne; hetero-Ne] and [homo-N.m; hetero-N.m]). As shown in figure 3, 95% of the datasets simulated under the supermodel "isolation" with random combinations of parameters from large priors are correctly recaptured by the random forest approach with a high probability (4% ambiguity and 1% error if we apply a posterior probability threshold of 0.84; table S1)). Similarly, 98% of the pseudo observed datasets under the "migration" supermodel are strongly recaptured (with 1% of ambiguity and 1% of error for a threshold of 0.665). Models with migration are globally more efficiently recaptured by DILS, relying on a lower threshold probability to be robustly supported. The second step of the hierarchical model comparisons is to classify the demographic models described in figure 1-B, within the supermodel that was best supported in the previous step (figure 3-B). The results of the performance analyses first show that a pseudo-observed dataset simulated under an SI model (homo-Ne and hetero-Ne) is very unlikely to be strongly supported in an SI versus AM comparison. Out of 20, 000 simulations, only 1% are correctly recaptured by DILS (98% ambiguity and 1% error for a threshold of 0.845). The AM model is more robustly supported than SI (56%), but the 10, 000 inferences made under each of the AM sub-models lead to weak support (43% ambiguity, 1% error for a threshold of 0.705).

Data

On the contrary, the two models making the "migration" supermodel (IM and SC) are more efficiently distinguished by DILS. The 40, 000 pseudo-observed datasets randomly simulated under the IM model are captured at 79% with a high probability in the IM versus SC comparison (20% ambiguity, 1% error for a threshold of 0.885). Similarly, 70% of the 40, 000 pseudo-observed datasets from the SC model are correctly recaptured by DILS (29% ambiguity, 1% error for a threshold of 0.915).

We then evaluate DILS performance for discriminating among alternative models of genomic distribution for the Ne (figure 3-B; table S2) and N.m (figure 3-C).

Concerning the effective population size, DILS systematically recaptures the homogeneous model more easily than the heterogeneous model for each of the four demographic models tested. The most complicated model to recapture is the genomic heterogeneity of Ne in an SC model (≈ 64% true positives), while homo-Ne under an SI model is the most straightforward.

Concerning the genomic model of introgression, the main parameter determining the quality of inference is the relative duration of gene flow versus speciation time T split (figure 1). This leads to a higher robustness for the IM model (probabilities of correctly supporting homo-N.m and hetero-N.m of ≈ 86% and ≈ 76% respectively) compared to the SC model (≈ 71% and ≈ 60%).

Detection of barriers to gene flow

One additional purpose of DILS is to identify genomic regions which are associated to a local reduction in the introgression rate N.m. This analysis will only be carried out by DILS if the observed dataset is better explained by 1) a demographic model with ongoing migration (IM or SC) and 2) a genomic model with gene flow heterogeneity (hetero-N.m). To achieve this purpose, DILS will infer the parameters under the model that best explains the data. Then, for each gene, DILS performs a comparison between two models that differ only for the parameter N.m: 1) the migration model corresponds to the whole set of parameters estimated under the best supported model; 2) the isolation model corresponds to the previous model whose N.m has been set to zero, because a barrier gene impedes gene flow locally along the chromosome. Therefore, such a gene should be supported by the isolation model if the barrier effect is strong. This approach therefore seeks to approximate a continuous variable, N.m, by a dichotomous choice of model: region with a local-isolation versus local-migration.

In order to evaluate the performance of DILS, the locus-specific model comparison was applied for locus simulated under an IM model with different values of N.m in [0, 1] (figure 4). A value of zero means no exchange during the divergence process from one population to another. A value of 1 means that there is one immigrant individual on average every genera-tion. We simulated genes 10, 000 times for different combinations of N.m and T split under an IM model. Then, for each simulated dataset, we applied the locusspecific model comparison to finally record for each locus which model is the best between local-migration and local-isolation. Ideally, we aim that DILS considers 100% of the simulations with N.m = 0 as localisolation, and 100% of the simulations with N.m = 1 as local-migration. Values of N.m greater than 1 were not explored because the comparison between "high migration" and "very high migration" is not relevant here.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 N.m prop. of inferred barrier genes T split q q q q q 0.05 0.5 1 5 10 As shown in figure 4, in the case of two populations of 100, 000 individuals separated only 5, 000 generations ago (T split = 0.05), DILS will support gene flow for ≈ 50% of the genes that have a N.m migration equals to zero (figure 4). For N.m = 1, the proportion of genes inferred as local-isolation is of similar mag-nitude, indicating that DILS is not at all designed to detect barriers to gene flow in the genomes of populations that have separated very recently. Our recommendation, therefore, is to disregard the results of DILS if the studied populations are extremely recent. However, as soon as barrier regions have enough time to differentiate (T split ≥ 0.5; figure 4), then ≈ 100% of loci with N.m = 0 are correctly inferred as local-isolation, and only few loci with N.m = 1 are incorrectly supported by the model of local-isolation. The performance of DILS therefore depends directly on the true history of the studied populations/species, not on the amount of data. The ideal case for identifying which genes in the genome are linked to barriers occurs when the patterns of polymorphism and divergence at such genes differ greatly from the rest of the genomic background (figure 4). An ideal demographic scenario for identifying barriers with DILS would be :

1. a divergence that is old enough to allow the neutral regions linked to barriers to have at least one position with two variants that are differentially fixed between the two populations/species. 2. a migration rate in the genomic background high enough to counteract the effect of differentiation in non-barrier regions.

Parameter estimates

In this section, we describe performance tests for estimating the parameters of different demographic models. The same procedure was applied for single and twopopulation models: first, simulating pseudo-observed datasets (10, 000 for the three single-population models, 2, 000 for the 14 two-population models) and then ABC estimation of the parameters to test DILS ability to recapture the parameter values used. We only detail here the results obtained for the demographic parameters, i.e., those describing the mean effective population size (N e current and N e past ), time of split (T split ), the date for the cessation of gene flow (T AM ), the age of the secondary contact (T SC ) and the migration rates (N.m).

Single-population models

The current effective population size is by far the parameter that is most accurately recaptured, especially in a constant and homogeneous model with a meansquared error (MSE) close to zero (MSE≈ 0.005; figure 5-A; table S2).

The introduction of a recent demographic change reduces the quality of inferences for N e current , more for the expansion model than for the contraction model.

The inference of ancestral size conducted for 4 x 10, 000 pseudo-observed datasets shows globally a low error rate on the raw values of N e past with a M SE max ≈ 0.09 (table S2). Errors depend very closely on the relative values between N e current and N e past (figure 5-B). Hence, while the estimate of a N e past is reliable for a change in size by a factor of 10, it becomes less accurate as the contrast with N e current is increasingly sharp. Figure 6: Parameter estimation for two-population models 2, 000 pseudo-observed datasets are simulated by taking random parameter values under the 14 models and analyzed using the same procedure as to produce the figure 5, but for the effective size of the current population N e current (A), the ancestral population size N e past (B), the time of split T split (C, in generations), the times of demographic transitions T AM and T SC (D) and the migration rate N.m (E).

In a similar manner, the quality of inferences of the age of demographic change T dem is highly dependent on its relative value with N e current (figure 5-C). Any change more recent than 0.1N e current generations ago will be dated with poor precision. Conversely, the age of events older than 0.1N e current appears more accurately recaptured by our ABC approach.

Two-population models

To estimate the accuracy of DILS in recapturing model parameters, we followed the same procedure based on the analysis of simulated datasets as in the previous paragraph. Note that the two-population models comprise two additional parameters that affect patterns of divergence: the time of split (T split ) and the migration rate (N.m).

The error rate in the estimation of the parameter N e current is of the same order of magnitude as in models with a single population (figure 6-A; table S3). However, the imprecision increases with ongoing migration and tends to underestimate N e current . Estimates of N e current are thus more accurate for the SI model, than for the AM model, and the worst for IM and/or SC. This negative effect of ongoing migration on the accuracy of parameter estimation is more pronounced for the ancestral population size N e past (figure 6-B). Hence, the ongoing migration implemented in the IM and SC models will lead to the overestimation of very low N e past values and underestimation of large N e past values.

The precision of the estimate of T split for a given model is of the same order of magnitude as for the ancestral size, with the exception of an accentuated imprecision of T split in the IM model when the migration is homogeneous along the genome (figure 6-C; table S3).

The AM and SC models both have an additional parameter describing the time of the demographic transition between two periods (with and without migration). In the AM model, T AM describes the number of generations during which the two current populations remain genetically isolated after a period of ancestral migration. Conversely, in the SC model, T SC describes the number of generations where the two current populations are connected by gene flow during a secondary contact occurring after a past period of isolation. For the AM model, T AM is better estimated than T split unlike T SC under the SC model (figure 6-D; table S3).

Finally, the performance of DILS to estimate the migration rate N.m (expressed in number of individuals immigrating per generation) is reported on the figure 6-E. The poor estimation accuracy for N.m contrasts sharply with the reliable inferences obtained when comparing the 'ongoing migration' versus 'current isolation' supermodels (paragraph 2.1.2). Indeed, it is straightforward to discriminate between these two cat-egories of supermodels while an accurate estimate of the migration rate is more challenging to obtain (figure 6-E). We were unable to reach a reliable measure of N.m for the AM model. More accurate inferences are obtained for both the IM and SC models. Hence, accuracy is reported to increase for genomic models where N.m is homogeneous (table S3).

Illustration of DILS with RNA-seq data from 28 pairs of Mytilus populations

We now illustrate the potential of DILS to contribute to the study of speciation among 8 populations of a complex of four Mytilus mussels species (figure 7-A; [START_REF] Bierne | Introgression patterns in the mosaic hybrid zone between Mytilus edulis and M. galloprovincialis[END_REF][START_REF] Popovic | Twin introductions by independent invader mussel lineages are both associated with recent admixture with a native congener in Australia[END_REF]. The relationship between molecular divergence and genetic isolation was established over 28 pairs of Mytilus populations using DILS. The aim was to identify which pairs of populations, characterized by different levels of molecular divergence (net synonymous divergence between 0.003% and 6.705%), are inferred to be currently connected by ongoing gene flow (figure 7-B). This large scale analysis within the same genus was made possible by the use of a large RNA-seq dataset recently published by [START_REF] Popovic | Twin introductions by independent invader mussel lineages are both associated with recent admixture with a native congener in Australia[END_REF] in 44 individuals of Mytilus from 8 populations (one M. edulis, five M. galloprovincialis, one M. planulatus and one M. trossulus; figure 7-A). Out of 28 pairs of Mytilus populations that have been tested for ongoing gene flow, 9 pairs receive support for models with current isolation, while models with migration suggest a better fit to the patterns of polymorphism and divergence observed in the remaining 19 pairs (figure 7-B). Within the group composed of M. galloprovincialis and M. planulatus, the 15 possible pairs are characterized by levels of divergence ranging from ≈ 0.003% (Crique -Herceq) to ≈ 0.673% (Primel -Spring). All of them are supported by models with ongoing gene flow. Our ABC analysis provides support for gene flow for a single interspecific pair with high divergence level ≈ 6% (Darling -Lighthouse). This pair is the only one to be genetically connected by migration among the 8 pairs in our analysis that have a level of net synonymous divergence higher than 2%, most probably due to gene flow between West Atlantic trossulus populations and edulis ones. Such gene flow also contrasts with an analysis realized for 61 pairs of animal species that had been studied without a priori on their speciation history, and for which a threshold of ≈ 2% had emerged above which all interspecific pairs were currently isolated (figure 7-B; [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF].

Discussion

Our statistical analysis platform, DILS, goes beyond simple summary statistics by explicitly testing evolu- https://cutt.ly/ErrDuwj). A median number of 1, 407 coding genes was used to perform demographic inferences after filtering the data (min = 144 genes; max = 2, 899 genes; depending on the pair of Mytilus considered). B: The x-axis shows the net divergence measured at the synonymous positions of the sequenced genes. The y-axis shows the probability provided by DILS of models with ongoing migration (IM + SC). Grey dots correspond to the 61 pairs of populations/semi-isolated species/species studied in [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. The coloured dots correspond to the 28 pairs of newly analyzed Mytilus. The colours refer to the labelled species from which the partners of each designated pair originate.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
tionary scenarios in a model-based inference framework. This approach is especially time-wise in speciation genomics as comparative studies of closelyrelated species are accumulating (e.g. in butterflies: [START_REF] Cong | Genomic determinants of speciation[END_REF][START_REF] Edelman | Genomic architecture and introgression shape a butterfly radiation[END_REF][START_REF] Martin | Recombination rate variation shapes barriers to introgression across butterfly genomes[END_REF]in birds: Peñalba, Joseph, and Moritz, 2019;in fishes: Malinsky et al., 2018;in plants: Stankowski et al., 2019); and so there is a strong demand for efficient and powerful inference tools. Ap-plied to genomic data (SFSs and/or multi-locus sequences), DILS will first identify the best demographic model to test for changes in effective size and migration rate over time, then it will identify the best genomic model to test for genome-wide heterogeneity in these parameters, and finally it will identify loci most associated with genomic regions locked to gene flow.

Performances of DILS

For single-population models, DILS is highly efficient at distinguishing the three demographic models (Expansion versus Constant versus Contraction). It fairly discriminates among the two genomic models (homo-Ne versus hetero-Ne) for the Expansion and Constant models, but has too much ambiguity to distinguish them in a Contraction model. Current population sizes (N e current ) are accurately estimated under all three demographic models, as well as the time of size change (T dem ) provided that it is not too recent. However, the past population size (N e past ) is increasingly overestimated in an expansion model (respectively, underestimated in a contraction model) when the contrast with N e current is increasingly sharp.

We also found that in two-population models, DILS very accurately discriminates between the two supermodels classically tested in speciation ("current isolation" versus "ongoing migration"), and it discriminates reasonably well among models with ongoing migration (IM versus SC), but quite poorly among models without (SI versus AM). Within each demographic model, the two Ne-genomic models are fairly discernible; and the same is true for the two N.m-genomic models (homo-N.m versus hetero-N.m) in scenarios of ongoing migration. Parameters are reasonably well estimated in all models (i.e. population sizes N e current and N e past , and times T dem , T AM and T SC ); except for the migration rate (N.m). It is poorly estimated in ongoing migration models, and cannot be evaluated at all when migration happened in the past (AM).

Therefore it is critical for users to be aware of the limits of DILS; especially when it does accurately discriminate among models and estimate parameter values, and when it does not. For example, the best scenarios for identifying barriers to gene flow is when the genetic signal for these genes strongly contrasts with the rest of the genome, i.e. when speciation time is long enough to build-up divergence at barrier regions, and migration rate is high enough to homogenize the genomic background between species. In general, DILS fails to make accurate inferences when divergence or changes in population size have occurred very recently.

Collaborative research

DILS was designed with the objective to facilitate collaborative research in speciation. One major question in the field is to understand how fast reproductive isolation builds-up with divergence between lineages, and so how fast introgression decreases along a continuum of molecular divergence. This relationship has been investigated in 61 pairs of animals [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]) only providing a partial picture. Here, we extended this work by analyzing genomic data of 28 species/populations of M ytilus mussels. Within this specific clade, we found a pattern of non-linear decrease of migration probability with the neutral molecular divergence, similar to what was observed in [START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. However, we also documented ongoing migration between two highly divergent mussel species, hence pushing the grey zone of speciation threshold beyond 2% of net synonymous divergence, maybe due to the outstanding life history traits of mussels (i.e. broadcast spawning, high-dispersal larvae, large effective population sizes and living in a highly connected marine environment).

DILS offers the possibility for the users to participate to this enterprise, and record on the web-platform where their biological model falls within this global speciation picture. Such a global picture of transition from gene flow to no gene flow is necessary for the central problem of species delineation (Hey and Pinho, 2012). Although a universal criterion for delineating species seems impossible, as exemplified by the mussel dataset, the idea of defining a grey zone by taxonomic system is promising [START_REF] Galtier | Delineating species in the speciation continuum: A proposal[END_REF]. Thus, our collaborative approach option included in DILS will allow in the future to establish a relationship between molecular divergence and genetic isolation for different taxa, i.e. vertebrates, terrestrial plants, algae, etc ..., and thus will provide delimitation rules by system.

Non-detailed features

The raw data can be easily visualized with DILS as a site frequency spectrum and summary statistics across loci. DILS produces comprehensive results for each inference step: (1) the global model comparison to estimate the best demo-genomic model; (2) the locusspecific model comparison to identify barrier loci; and

(3) the estimation of parameter values for the best model. To help users interpreting these results, DILS produces a series of goodness-of-fit tests to the data. These tests are performed by simulating under the best model each population genetic statistic calculated in section 4.1.1 (genomic mean and variance of π, θ, F ST , etc.), as well as for each bin of the SFS (or jSFS for two-population models). In addition to an individual test for each summary statistic, a test is also performed from statistics transformed by a PCA following [START_REF] Cornuet | Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation[END_REF][START_REF] Cornuet | DIYABC v2. 0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data[END_REF]. DILS also provides values for each locus of: 1) each summary statistic; 2) the approximated recombination rate calculated based on the four-gamete rule [START_REF] Galtier | Codon usage bias in animals: disentangling the effects of natural selection, effective population size and GC-biased gene conversion[END_REF][START_REF] Hudson | Statistical properties of the number of recombination events in the history of a sample of DNA sequences[END_REF]; and 3) the posterior probability of being genetically linked to a barrier to gene flow (for two-population models only). These results are outputted as interactive graphics on the web-platform.

Our method is implemented in a user-friendly webplatform allowing the configuration of the ABC analysis via a graphical interface, its execution and the visualization of the results. Detailed information for how to use DILS is provided in the manual. The released version of DILS is currently hosted by the French Institute of Bioinformatics (XXX). To ensure full reproducibility and portability on any server, DILS is packaged in a singularity container freely available at https://github.com/popgenomics/DILS_web. The complete analysis of a dataset (model comparison + parameter estimates + locus-specific tests + goodnessof-fit tests) on the host server takes 4h30 on average.

Prospects for the future

With the improvement of computational methods, it is now possible to simulate entire chromosomes under the full ancestral process of coalescence and recombination [START_REF] Kelleher | Efficient coalescent simulation and genealogical analysis for large sample sizes[END_REF]. Combining this type of coalescent simulators with haplotype-based statistics in our ABC framework would be very promising to improve estimates of the timing and extent of gene-flow after secondary contact [START_REF] Harris | Inferring demographic history from a spectrum of shared haplotype lengths[END_REF]. The architecture of DILS has been designed to easily add simulators other than ms and its modified versions [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF]. Thus, it would be readily achievable to use forward-in-time simulations including direct selection [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF], and therefore making inferences for any selective scheme while taking into account the demographic history of the sample, without changing the pipeline upstream or downstream of the simulator.

Materials and methods

ABC

Summary statistics

Since ABC is a category of inferential method based on the comparison between statistics summarizing simulated and observed datasets, we first describe here the statistics computed in our framework.

We assume that users are interested in carrying out inferences from multilocus datasets. For single population models, DILS calculates for each locus: i) the pairwise nucleotide diversity (π) [START_REF] Tajima | Evolutionary relationship of DNA sequences in finite populations[END_REF]; ii) Watterson's θ [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF] and iii) Tajima's D (Tajima, 1989). In addition to these three statistics, the sitefrequency spectrum (SFS; [START_REF] Fischer | The genetical theory of natural selection[END_REF][START_REF] Wright | Evolution in Mendelian populations[END_REF]Wright, 1938) is also used to summarize the data, i.e, the number of single-nucleotide polymorphism (SNP) where the derived allele is present in [2, ..., nseq -1] copies in the studied population/species, where nseq represents the number of copies sampled for a given locus. If the SFS is folded by the absence of an outgroup, then the SFS will be described by the number of SNPs where the minor allele is present in [2, ..., nseq/2 -1] copies. Finally, for single population models, multilocus inferences are based on 6 multilocus summary statistics which are the means and standard deviations of π, θ and Tajima's D, to which we add [nseq -2] individual statistics corresponding to the SFSs summed over all loci. In absence of an outgroup, the SFS will be represented by [nseq/2 -2] individual statistics.

For models with two populations/species, π, θ and Tajima's D statistics are also calculated for each of the two samples. These are supplemented with statistics approximating the joint SFS (jSFS; [START_REF] Ramos-Onsins | Multilocus analysis of variation and speciation in the closely related species Arabidopsis halleri and A. lyrata[END_REF]: 1) the fraction of sites showing a fixed difference between the populations/species (Sf ), 2) the fraction of sites showing an exclusive polymorphism to a given population/species (Sx A and Sx B ) and 3) the fraction of sites with a polymorphism shared between the population/species (Ss). Statistics describing the divergence between the two populations/species are also calculated, including the raw (D xy ; [START_REF] Nei | Molecular Evolutionary Genetics Columbia University Press New York[END_REF] and the net (D a ; [START_REF] Nei | Mathematical model for studying genetic variation in terms of restriction endonucleases[END_REF] divergence between the population/species, and their relative genetic differentiation measured by F ST (Wright, 1943). Finally, all bins in the jSFS (except singletons) can optionally be used as an additional vector of summary statistics. If the jSFS is unfolded, then this vector has a length of [(nseq A +1) * (nseq B +1) -4] available statistics (minus 4 to remove the two bins corresponding to singletons and the two bins corresponding to the fixation of the derived or ancestral allele in both samples), and a length of [(nseq A + 1/2 * nseq B + 1/2 -3] if the jSFS is folded. The use of the jSFS as a vector of summary statistics is an option that the user can switch on or off, to avoid cases where the jSFS is composed by a large number of bins.

Simulations

In the current version of DILS, all simulations are performed using the msnsam coalescent simulator [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF][START_REF] Ross-Ibarra | Patterns of polymorphism and demographic history in natural populations of Arabidopsis lyrata[END_REF]. Each simulated multilocus dataset takes properties from the observed datasets (same number of genes, lengths and sample size). Since the summary statistics used to perform the ABC inferences are averages and standard deviations measured over all the surveyed genes, then for model comparisons and parameter estimations we randomly subsample 1, 000 genes if more loci are present in the total dataset. The purpose of this subsampling is to avoid unnecessarily long simulation times because the values of statistics for a given locus will not impact the used summary statistics over 1, 000 loci.

If an outgroup is specified by the user, then it will be used for each locus (or contig, or gene) to correct its mutation rate µ i to µ. divi div where µ is the neutral mutation rate assumed by the user, div i is the raw divergence between the focal population/species and the outgroup measured at a given locus i, and div is the average raw divergence between the focal population/species and the outgroup measured over all loci. The other implication of using an outgroup will be to orientate the mutations and consequently to unfold the jSFS. Finally, the loci are assumed to be genetically independent, and a ρ θ ratio value has to be specified by the user where ρ is the average population recombination rate 4.N e.r (r in number of recombination events per generation and per nucleotide).

Model comparisons

Here, when used alone, the term model means a given combination between a demographic and a genomic model. All comparisons are performed by using the abcrf function of the eponymous R package [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF]. The comparison is a two-step process.

First, grow the random forest with the abcrf function. This requires one reference table per model for the training. The reference table of each model is produced by 10, 000 multilocus simulations whose parameters correspond to random combinations sampled from priors. They are composed of one row per multilocus simulation and one column for each summary statistic described in section 4.1.1. When categories of models are compared following the hierarchical approaches (figures 2 and 3), the reference tables of the different models in the same category are merged together. For instance, in the comparison between Current isolation and Ongoing migration (figure 3), 60, 000 multilocus simulations are used for the training of the super-model Current isolation, and 80, 000 multilocus simulations for the training of Ongoing migration. Each forest is made up of 1, 000 grown decision trees regardless of the comparison made throughout the hierarchical approach.

The last step is the prediction of the best model among those proposed by passing the observed data through the trained random forest.

DILS reports the model supported by the largest number of decision trees in the random forest and its associated posterior probability.

Parameter estimates

Two strategies are applied simultaneously to estimate the parameters describing the best-supported model among those compared:

1) a joint estimation of the set of parameters using a rejection/regression method [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF]. Estimation is based on the 5, 000 multilocus simulations producing the statistics closest to the observed dataset among 1, 000, 000 simulations. We then correct for imperfect matches between observed and retained values of statistics. The parameter values of the selected simulations are weighted by their Euclidean distance and corrected according to a nonlinear regression method using a neural network. 10 trained neural networks with 10 hidden networks are used in the regression.

2) an individual estimation of each parameter by constructing a random forest of 1, 000 trees per parameter [START_REF] Raynal | ABC random forests for Bayesian parameter inference[END_REF].

The results from both approaches are returned to the user. There is no evidence to further support a method over the other in terms of estimation accuracy. Within the framework of the models currently compared in DILS, both approaches produce similar estimates when tested on pseudo-observed datasets. However, joint parameter estimation has the advantages of including parameter co-variations as well as providing a probability density. This is achieved at a computational cost that is ≈ 100 times greater regarding the number of multilocus simulations, since 10, 000 are required for a parameter estimation using random forest versus 1, 000, 000 when using the rejection/regression algorithm.

Locus-specific model comparison

To identify barriers to gene flow among a set of sequenced DNA fragments (genes for instance), we adopt the same procedure as in Ciona intestinalis [START_REF] Roux | Crossing the species barrier: genomic hotspots of introgression between two highly divergent Ciona intestinalis species[END_REF] and Mytilus [START_REF] Roux | Can we continue to neglect genomic variation in introgression rates when inferring the history of speciation? A case study in a Mytilus hybrid zone[END_REF] but by replacing the neural network with a random forest to divide the computational cost by 100. This step is performed by DILS only if 1) observed data better fits models with ongoing migration (IM or SC; figure 1) and 2) genomic models of N.m. variation explain the data better than homogeneous models. We first estimate the parameters of the best model from the multilocus dataset. Based on this estimation, two models are compared at each locus: 1) local-migration: the multilocus estimated model with the non-zero migration rate estimated over the whole genome. 2) local-isolation: the multilocus estimated model with a migration rate set to zero. A random forest of 1, 000 trees is then trained to recognize combinations of summary statistics specific to each of the two evaluated models. This forest allows to return for each sequenced DNA fragment the locus-specific model that best explains the statistics observed, and its posterior probability.

Pseudo-observed datasets

In this study we distinguish two types of simulations. Simulations carried out to build reference tables used to train a random forest, or from which a small proportion will be sub-sampled in a rejection/regression algorithm based on the Euclidean distance with the observed data. Then a second type of simulations produces pseudo-observed datasets. These are not used for training, but to evaluate the inferential power of the ABC approach, and test whether it can recapture the parameters used to simulate the pseudo-observed datasets. To assess the reliability of model comparisons

( θ i -θ i ) 2 ,
where θ i and θ i represent the estimated and the true parameter values respectively

Figure 1 :

 1 Figure 1: Demographic models currently implemented in DILS A: compared single-population models. Demographic changes occurring T dem generations ago are modeled backwards in time by sudden transitions from N e current to N e past , either for expansions or contractions. B: compared two-population models. The Strict Isolation (SI) and Ancient Migration (AM) models are characterized by an absence of ongoing migration.Conversely, the Isolation with Migration (IM) and Secondary Contact (SC) models describe two populations that are currently connected by introgression events at rate N.m. The two-population models shown here are of constant size, but DILS optionally incorporates alternative versions of the same four models where effective size can change independently in both daughter populations between the present time and T split .

Figure 3 :

 3 Figure 3: Performance of DILS for hierarchical comparison of two-population models Two-population analyses are performed in three steps 1. Testing the general demographic models "Current isolation" versus "Ongoing migration" (panel A) 2 Testing the demographic sub-models (panel A): 2.a if the best model is 'current isolation', then DILS tests SI versus AM 2.b if the best model is 'ongoing migration', then DILS tests IM versus SC 3. Testing the genomic models for variation of: 3.a effective population size, N e (panel B) 3.b migration rate, N.m (panel C) The letters 'o' and 'e' in panels B and C indicate simulations performed under genomic homogeneity and heterogeneity models, respectively The pie charts designate for each model the proportion of simulations performed under the corresponding model that is strongly and correctly captured (correct: blue), strongly and incorrectly captured (wrong: yellow) and without strong statistical support for any of the studied models (ambiguous: purple).

Figure 4 :

 4 Figure 4: Detection of barriers to gene flow x-axis : 11 explored values of the locus-specific N.m migration rate under an IM model. y-axis : proportion of simulations supported by DILS as being linked to a barrier to gene flow. The colors designate five different divergence times of the IM model (T split , figure 1). The unit time is in Ne generations where Ne is the number of haploid individuals making up the population. If N e = 100, 000 individuals, then T split = 5 means a divergence time of 500, 000 generations under the IM model. If Ne is the number of diploids, then T split must be multiplied by two to find the same relationship. Each combination of T split and N.m was independently simulated 10, 000 times and analyzed by DILS to get the proportion of model-assignation for a given combination of parameters. The estimated points are connected by dotted lines for visibility.

  Figure 5: Parameter estimation for singlepopulation models 10, 000 pseudo-observed datasets are simulated by taking random parameter values (x-axis) under the 6 models. These parameters are estimated using DILS (y-axis). The lines represent the loess regressions between exact and estimated parameter values for each of the six models. The fields represent the 99% confidence interval and the dotted line represents x = y. Estimation of the effective size of the current population N e current (A), of the ancestral population N e past (B) and the time of demographic changes T dem (C). N e past and T dem are both expressed here in N e current individuals.

10⁶

Figure 7 :

 7 Figure 7: Application of DILS to a Mytilus RNA-Seq dataset A: Transcriptomes were obtained from Popovic et al., 2019 from 44 individuals sampled from 4 labelled species (throughout 8 localities), providing 28 possible pairs of Mytilus localities analysed to test for a genetic connection through migration events (SRP218536;https://cutt.ly/ErrDuwj). A median number of 1, 407 coding genes was used to perform demographic inferences after filtering the data (min = 144 genes; max = 2, 899 genes; depending on the pair of Mytilus considered). B: The x-axis shows the net divergence measured at the synonymous positions of the sequenced genes. The y-axis shows the probability provided by DILS of models with ongoing migration (IM + SC). Grey dots correspond to the 61 pairs of populations/semi-isolated species/species studied in[START_REF] Roux | Shedding light on the grey zone of speciation along a continuum of genomic divergence[END_REF]. The coloured dots correspond to the 28 pairs of newly analyzed Mytilus. The colours refer to the labelled species from which the partners of each designated pair originate.

simulations to train the random-forest: 10

	• Data cleaning:
	forall genes i making the dataset do
	.discard from the alignment of gene i the
	individuals with too many Ns;
	.discard gene i if there are not enough
	retained individuals;
	.discard gene i if it doesn't contain enough
	positions without an aligned null allele
	carried out to build reference tables used to
	train a random forest, or from which a small
	proportion will be sub-sampled in a
	rejection/regression algorithm based on the
	Euclidean distance with the observed data.
	Then a second type of simulations produces
	pseudo-observed datasetitN;
	end
	• Reference , 000 multilocus datasets
	under each combination of [demographic
	models] x [genomic models];
	• Model comparisons:
	Random Forest (RF) comparisons of:
	step 1. Expansion [homo + hetero] versus
	Constant [homo + hetero] versus Contraction
	[homo + hetero];
	step 2. best demographic model supported in the
	previous step with homogeneous Ne versus
	heterogeneous Ne;
	Algorithm 1: Single-population hierarchical model
	comparisons

  : A single fasta containing all genes sequenced in all individuals sampled in the two studied populations Result: Posterior probabilities for 1) ongoing

	hetero-Ne
	IM + hetero-Ne SC];
	step 4. if migration then [homo-N.m IM +
	homo-N.m SC] versus [hetero-N.m IM +
	hetero-N.m SC];
	Algorithm 2: Two-population hierarchical model
	comparisons

migration, 2) [SI, AM] (in case of current isolation) or [IM, SC] (in case of ongoing migration), 3) Ne and N.m (in case of ongoing migration) • data cleaning: forall genes i do forall population/species j do .discard from the alignment of gene i the individuals with too many Ns; .discard gene i if there are not enough retained individuals in population j; .discard gene i if it doesn't contain enough positions without an aligned null allele N; end end • Reference simulations to train the random-forest: 20, 000 multilocus datasets under each combination of [demographic models] x [genomic models]; • Model comparisons:; Random Forest (RF) comparisons of:; step 1. isolation ([all SI + all AM]) versus migration ([all IM + all SC]); step 2. if isolation then [all SI] versus [all AM] else [all IM] versus [all SC]; step 3. if isolation then [homo-Ne SI + homo-Ne AM] versus [hetero-Ne SI + hetero-Ne AM] else [homo-Ne IM + homo-Ne SC] versus [

New South Wales Crique, Les Issambres Darling Marine Station, Maine Herceg Novi Lighthouse Park, British Columbia Primel Spring Bay, Tasmania Sydney Harbour, New South Wales

	•	•	•	•	•
			Batemans Bay,	• • •

Supplementary Table 1: Performance of DILS model comparisons

  For each model, 10, 000 sets of pseudo-observed data were analysed. These datasets were simulated by taking random combinations of parameters from large prior distributions. The table reports for each model comparison the proportions among these simulations that lead to correct, ambiguous or wrong inferences according to a threshold set to keep the rate of wrong inferences below 1%.Estimation errors are calculated on N = 2, 000 simulated datasets for each model, by using random combination of parameter values. The reported values are measured as follows 1 N * variance(θ i )

		comparison				target	correct	ambiguous	wrong	threshold
		Single-population models
	Expansion versus Constant versus Contraction	Expansion	0.8627	0.1275	0.0098	0.919
						Constant	0.891	0.099	0.01	0.856
					Contraction	0.9937	0.0024	0.0039	0.501
		Expansion homo-Ne versus hetero-Ne			homo-Ne	0.8248	0.1652	0.01	0.83
						hetero-Ne	0.7492	0.2408	0.01	0.92
		Constant homo-Ne versus hetero-Ne			homo-Ne	0.8562	0.134	0.0098	0.835
						hetero-Ne	0.8336	0.1564	0.01	0.89
		Contraction homo-Ne versus hetero-Ne			homo-Ne	0.4074	0.5828	0.0098	0.906
						hetero-Ne	0.2634	0.7266	0.01	0.901
		Two-populations models
		migration versus isolation			migration	0.97922	0.01078	0.01	0.665
						isolation	0.95053	0.03968	0.00978	0.84
		SI versus AM				SI	0.01465	0.9754	0.00995	0.845
						AM	0.81482	0.18518	0	0.705
		IM versus SC genomic model N e	N e (current)	N e (past)	IM SC T split	0.791 0.69958 T AM	0.19908 0.29042 T SC N.m	0.00992 0.01 number barriers of	0.885 0.915 α	β
	AM	SI homo-Ne versus SI hetero-Ne homo homo 0.02114 0.09633 0.08024 0.02503 homo-Ne 0.9688	0.0213 1.01278	0.0099	0.8
		hetero	hetero-Ne 0.02832 0.10649 0.09502 0.02797 0.8173	0.1741 1.0092	0.0086	0.96 0.39736 0.96193
		AM homo-Ne versus AM hetero-Ne hetero homo 0.02134 0.09238 0.08401 0.04687 homo-Ne 0.95495	0.0355 1.02054 1.01096 0.00955	0.82
		hetero	hetero-Ne 0.03793 0.21534 0.22518 0.07039 0.8156	0.17445 1.00063 0.94219 0.46799 0.97759 0.00995 0.95
	IM	IM homo-Ne versus IM hetero-Ne homo homo 0.08547 0.30056 0.60409 homo-Ne	0.7006	0.29005 0.37401	0.00935	0.855
		hetero	hetero-Ne 0.09235 0.35789 0.62885	0.6469	0.3433 0.40637	0.0098	0.93 0.452	0.97489
		SC homo-Ne versus SC hetero-Ne hetero homo 0.05707 0.31596 0.27658 homo-Ne	0.7063	0.284 0.6166	0.0097 0.52701	0.85
		hetero	hetero-Ne 0.07003 0.29829 0.27734	0.64645	0.34445 0.65794 0.55112 0.55623 0.97603 0.0091 0.925
	SC	IM homo-N.m versus IM hetero-N.m homo homo 0.08764 0.23384 0.19008 homo-N.m	0.85915	0.13105 0.33885 0.50572	0.0098	0.845
		hetero	hetero-N.m 0.09052 0.26734 0.21924	0.759	0.23265 0.35622 0.50668	0.00835	0.95 0.46972	0.9829
		SC homo-N.m versus SC hetero-N.m hetero homo 0.07016 0.25806 0.14854 homo-N.m	0.7144	0.27585 0.30055 0.78375 0.49972 0.00975	0.855
		hetero	hetero-N.m 0.07937 0.26373 0.15774	0.59875	0.39165 0.32815 0.77411	0.0096 0.5055	0.915 0.49075	0.9738
	SI	homo	0.0183	0.04164 0.01158
		hetero	0.02918 0.05768 0.01458				0.38479 0.95003
							N
							i=1

and parameter estimates, for single and two population models, we simulate pseudo-observed datasets consisting of 100 loci, of length equal to 1, 000 nucleotides, sampled from 10 diploid individuals in each population/species and a mutation rate of 5.10 -8 mutations/nucleotide/generation. These datasets are simulated according to demographic histories using random combinations of parameters from the priors.

Analysis of the Mytilus dataset

We downloaded the raw RNA-seq data deposited to the NCBI sequence read archive (BioProject ID: PR-JNA560413; https://cutt.ly/OtQN1Y0) by [START_REF] Popovic | Twin introductions by independent invader mussel lineages are both associated with recent admixture with a native congener in Australia[END_REF]. The raw data consist in a total of ≈ 145Gb from the sequenced transcriptomes of 47 mussel individuals. Three individuals from the M. californianus species were removed as they do not belong to the M. edulis complex. The reference transcriptome used for the mapping is made up of 16, 151 CDS, for a total length of ≈ 23M b. The reference was indexed using bowtie2 (version 2.2. 6 Langmead and Salzberg, 2012). For each individual, reads were aligned to the reference with bowtie2, and cleaned using samtools with a mapping quality threshold of 20 (version 1.3.1; [START_REF] Li | The sequence alignment/map format and SAMtools[END_REF]. Individual genotypes were called using reads2snp [START_REF] Tsagkogeorga | The population genomics of a fast evolver: high levels of diversity, functional constraint, and molecular adaptation in the tunicate Ciona intestinalis[END_REF] at positions covered by at least 8 reads. We then ran DILS for each of the 28 possible pairs of localities by tolerating up to 20% of missing data, rejecting genes with less that 100 codons without a missing data, and by keeping 6 copies per genes within each population/species. Simulations were conducted by exploring uniform priors for effective population sizes between 0 and 500, 000 diploid individuals, times of different demographic events (split, secondary contact, arrest of migration) between 0 and 1, 750, 000 generations. Presentation of the results were carried out with R [START_REF] Wickham | Welcome to the tidyverse[END_REF][START_REF] Chang | shiny: Web Application Framework for R. R package version 1[END_REF][START_REF] Sievert | [END_REF][START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF].

Acknowledgements

XXX

Supplementary Table 2: Mean-squared error in parameter estimations for single-population models