The divergence history of European blue mussel species reconstructed from Approximate Bayesian Computation: the effects of sequencing techniques and sampling strategies

Introduction

The biodiversity we inherited from the Quaternary was shaped by the process of species formation (Hewitt, 2000). A long-standing question concerns the timing and rate of gene exchange that occurred while populations diverged, during the incipient stages of speciation. Model-based inferences from genetic data have been used to investigate the history of gene flow (Beaumont et al., 2010). Special attention has been paid to the distinction between recent divergence in a strict isolation (SI) model, and older divergence with continuous migration (Nielsen & Wakeley, 2001), although of course, more complex scenarios are also possible (Marino et al., 2013, Sousa & Hey, 2013).

With next-generation sequencing technologies, thousands of SNPs throughout the genome can be used to infer the demographic histories of non-model species pairs (Sousa & Hey, 2013). One way of summarizing the information in these data is the unfolded joint site frequency spectrum (jSFS), that is, the number of copies of derived alleles found in each of the two sampled species. A recent and fast maximum-likelihood method based on the jSFS (Gutenkunst et al., 2009) has proven useful for distinguishing continuous migration from SI (e.g., in ragworts, Chapman, Hiscock & Filatov, 2013, and beach mice, Domingues et al., 2012). The method can also evaluate more complex scenarios (e.g., in sea bass, Tine et al., 2014; poplars, Christe et al., 2017; and whitefish, Rougeux, Bernatchez & Gagnaire, 2017), but it struggles to explore the parameter space in these cases. In addition, the method is not well suited for transcriptome data as model comparison by log-likelihood ratio tests assumes independence of SNPs. As a consequence simulations need to be conducted to evaluate competing models, and the computational speed advantage is lost.

As an alternative, Approximate Bayesian Computation (ABC) is a method based on simulations that avoids the need to explicitly compute the likelihood (Beaumont, Zhang & Balding, 2002). As such, histories of speciation characterized both by periods of SI and periods of gene exchange can easily be investigated, for example, the scenarios of ancient migration (AM) and secondary contact (SC). These scenarios can be extended by including two cycles of “isolation/gene exchange,” following climatic changes in the Pleistocene (Fig. 1). Methods have also been developed to include genome-wide heterogeneity in migration rates (Sousa et al., 2013; Roux et al., 2013). This is consistent with the “genic view” of speciation (Wu, 2001), whereby barriers to gene flow are often semi-permeable, varying in strength across the genome due to linked selection and recombination (Barton & Bengtsson, 1986). A major challenge in ABC, as compared to explicit likelihood methods, is the selection of summary statistics, which involves a trade-off between loss of information and reduction of dimensionality. Several methods have been suggested to select the most appropriate statistics for a given dataset and a set of models (Wegmann, Leuenberger & Excoffier, 2009; Nunes & Balding, 2010; Aeschbacher, Beaumont & Futschik, 2012); but most of these summarize the SFS (but see, e.g., Boitard et al., 2016) leading to a loss of information. Recently, several ABC studies have used the SFS directly to reconstruct the history of single populations (Boitard et al., 2016), or multiple populations (Xue & Hickerson, 2015; Smith et al., 2017). However, the number of statistics, that is, the number of classes in the spectrum, increases quadratically with the number of haploid genomes sampled in the case of a two-dimensional SFS, and even faster if more than two populations are considered. For this reason, recent studies have tested different ways of binning the SFS to circumvent the curse of dimensionality (Smith et al., 2017). While accumulating new data and developing new methods, studies that evaluate the impact of the sampling strategy and the inferential method to reconstruct species divergence history will be increasingly valuable (Li & Jakobsson, 2012; Robinson et al., 2014; Shafer et al., 2015; Cabrera & Palsbøll, 2017; Smith et al., 2017).

Mytilus edulis (Linnaeus, 1758) and M. galloprovincialis (Lamarck, 1819) are two closely-related species that currently hybridize where their ranges overlap along the Atlantic French coasts (Bierne et al., 2003a) and the British Isles (Skibinski, Beardmore & Cross, 1983). Their inter-specific barrier to gene flow is semi-permeable, and it has been shown to involve multiple isolating mechanisms, both pre-zygotic (e.g., assortative fertilization and habitat choice; Bierne et al., 2002, Bierne, Bonhomme & David, 2003b), and post-zygotic (hybrid fitness depression, Simon, Bierne & Welch, 2017). Other evidence of ongoing gene flow between M. edulis and M. galloprovincialis comes from footprints of local introgression of edulis-derived alleles into a population of M. galloprovincialis enclosed within the Atlantic hybrid zone (Fraïsse et al., 2014). Another study (Fraïsse et al., 2016) revealed that the Atlantic population of M. galloprovincialis was more introgressed than the Mediterranean population on average. At some specific loci, however, the Mediterranean population was found to be fixed for edulis alleles while the Atlantic population was not introgressed at all, suggesting that an ancient contact between M. edulis and Mediterranean M. galloprovincialis occurred during glacial periods. Finally, direct model comparisons have been conducted with the IMa method of Hey & Nielsen (2007), and with an ABC framework, and shown that M. edulis and M. galloprovincialis have experienced a complex history of divergence punctuated by periods of gene flow in Europe (IMa: Boon, Faure & Bierne, 2009, ABC: Roux et al., 2014, 2016).

Here, we analyzed coding sequence datasets from this well-known pair of sister species in Europe, and systematically reconstructed its speciation history by ABC, using different sampling strategies. By different sampling strategies, we mean that we varied (i) the number of individuals sampled (two, four or eight), (ii) the number of SNPs (which were obtained by different sequencing techniques, “exome capture” vs. “rna-seq”), and (iii) the number of bins in the jSFS (binning the spectrum into four, seven or 23 classes). We then evaluated the influence of these choices on model selection, using eleven distinct scenarios of speciation. Our results show that the influence on inferences of the number of individuals and loci sampled is surprisingly limited, while the different ways of binning the jSFS strongly affect the results. This suggests that demographic reconstructions are nowadays more limited by the way data are analyzed than by their volume. Moreover, we find that an history of periodic connectivity, with both ancient and contemporary introgression, and a semi-permeable barrier to gene flow, best fit these data, arguing for the development of flexible inference methods to better describe complex divergence histories (Simon & Duranton, 2018).

Materials and Methods

Data analysis

We restricted our analysis to loci assembled in all individuals, and longer than 300 bp after filtering positions containing missing data, or more than two segregating alleles when the three species were aligned. Only synonymous positions were used. The total number of loci and SNPs retained in each dataset is given in Table 2.

Site frequency spectrum

We first computed the jSFS for each dataset (Fig. 2 and Fig. S1). Each derived allele, oriented by treating the outgroup sequence as a fixed ancestral allele, was assigned to one cell of the jSFS, depending on its frequency in each of the two populations. From the full jSFS, different classes of polymorphism were extracted and used as summary statistics (Tellier et al., 2011). Specifically, we used the four Wakeley–Hey classes (jsfs = 4 in Fig. 2): fixed differences, Sf; private polymorphisms for each species Sx₁ and Sx₂; and shared polymorphisms, Ss (Wakeley & Hey, 1997). We also considered a summary in which the Sf and Sx classes were split depending on whether the derived allele was fixed or absent in the other species (jsfs = 7 in Fig. 2; Ramos-Onsins et al., 2004). The third decomposition of the jSFS contains 23 classes of polymorphisms because singletons and doubletons in each population were included as new classes (jsfs = 23 in Fig. 2). This corresponds to the full jSFS with n = 2 diploid individuals.

Estimators of polymorphism and divergence

We then computed a set of genetic statistics across loci to make use of the coalescent information contained within each sequence, instead of considering each SNP separately. Following previous studies (Fagundes et al., 2007; Ross-Ibarra et al., 2008; Roux et al., 2011), we used the following statistics: (i) nucleotide diversity, π₁ and π₂ (Tajima, 1983); (ii) Watterson’s θ_w1 and θ_w2 (Watterson, 1975); (iii) total and net inter-specific divergence, div and netdiv; (iv) between-species differentiation, FST, computed as 1 − π_S/π_T, where π_S is the average pairwise nucleotide diversity within species and π_T is the total pairwise nucleotide diversity of the pooled sample across species. We also included the four Wakeley–Hey’s classes as explained above (Sf, Sx₁, Sx₂, and Ss). Finally, we assessed departure from mutation/drift equilibrium using Tajima’s D₁ and D₂ (Tajima, 1989a, 1989b). The average and standard deviation across the loci of these statistics were calculated with the program MScalc (available from http://www.abcgwh.sitew.ch/, see Roux et al., 2011), and their values are given in Table 2 (“mscalc”).

Results

Discussion

NGS data give us the opportunity to capture the diversity of coalescent histories across loci, and so to reveal the complexity of the speciation process (Sousa & Hey, 2013). Recently, important efforts have been made to develop statistical methods of inference making use of population genomics data. Computing the jSFS is an efficient way of summarizing the demographic information contained in NGS data because anonymous SNPs can be used (e.g., produced by RAD-sequencing) and it does not rely on phased data. However, the jSFS obtained from low-coverage sequencing data (typically <10× per position and per individual) can be biased toward a deficit of rare variants; and this is of particular concern when investigating the demographic history of populations (Nielsen et al., 2012; Han, Sinsheimer & Novembre, 2013). A first category of maximum-likelihood methods uses forward diffusion theory to compute numerical solutions to the jSFS under complex models (see Gutenkunst et al., 2009; Lukic & Hey, 2012). A second category of methods estimates the expected jSFS under any demographic model, as simulations are used to approximate the composite-likelihood of the data (see Naduvilezhath, Rose & Metzler, 2011; Excoffier et al., 2013). Here, we used a likelihood-free method (ABC) on the jSFS of coding sequences, and we evaluated the influence of different sampling strategies on the inference of speciation models in mussels.

Our ABC-based model comparison shows little qualitative effect of individual and SNP sampling on the outcomes. First, inferences based on n = 2, n = 4, or n = 8 individuals supported the same model of PSC (Table 3), consistent with previous results that relatively few individuals are sufficient to make robust inferences on the way divergence occurs between lineages (Robinson et al., 2014). This is because most coalescence events occur recently in the population history; so increasing the number of individuals is only helpful to characterize recent demographic events rather than the past divergence history. Second, we showed that inferences from the two sequencing techniques (“exome capture” vs. “rna-seq”) were qualitatively the same, despite the very different number of SNPs that were sampled, again supporting the PSC scenario (Table 3). This implies that neither the sequencing technique, nor the number of informative sites have a substantial effect on the inference. In fact, the jSFS calculated from the different data sets were very similar (Fig. S1); we only found a deficit of divergent SNPs in the “exome capture” data that may be due to the difficulty of mapping highly divergent alleles onto a single reference (on the contrary, the “rna-seq” reads from the different species were independently assembled and mapped). From a theoretical perspective, adding more loci (or longer loci) provide information about deep coalescent events which are important for shedding light on the divergence history of closely-related species (Wang & Hey, 2010). In fact, previous simulation studies showed an influence of locus length and locus number on the ABC performance in model choice, and highlighted a threshold effect above which adding more loci did not significantly improve inferences (Li & Jakobsson, 2012; Robinson et al., 2014; Shafer et al., 2015). Thereby, we argue that the number of SNPs sampled in our study was sufficient to accurately represent the diversity of coalescent histories in the Mytilus genome, and consistently support the same speciation model.

In most studies, functions of the jSFS are used as summaries of the data. For example, the likelihood method of Nielsen & Wakeley (2001) uses four classes of polymorphisms to estimate migration rates and divergence times in the IM scenario. In the ABC approach, choosing a suitable set of summary statistics is difficult because it implies a trade-off between loss of information and reduction of dimensionality. Accordingly, practical methods to identify approximately sufficient statistics have been developed (Wegmann, Leuenberger & Excoffier, 2009; Nunes & Balding, 2010; Aeschbacher, Beaumont & Futschik, 2012). Here, we compared ABC inferences based on 23-binned jSFS (jsfs = 23) with inferences using a subset of polymorphic classes (jsfs = 4 and jsfs = 7). Model checkings through simulations pointed out the loss of information when only four or seven classes in the jSFS were considered; particularly for the inference of recent migration events. By decomposing the jSFS into 23 classes, we could reveal the excess of shared polymorphisms that are at high frequency in one species and low frequency in the other, a pattern produced by recent migrants (Table 3). This is in agreement with the simulation study of Tellier et al. (2011) that showed a significant improvement in the estimation of the timing of gene flow when these additional polymorphic classes were considered; and the study of Alcala et al. (2016) that showed an excess of high frequency derived alleles is the characteristic footprint of secondary contacts. Across all models, we showed that the probability to correctly infer the true model (accuracy rate) increases with the number of classes considered, and the ambiguity rate correspondingly decreases (Table S3). Smith et al. (2017) similarly found that the statistical power of their ABC model selection increases with the number of classes until reaching a plateau of error rates (but above which computation efforts continues to increase).

Inferences based on the 23-bin jSFS consistently support a model of PSC with genomic heterogeneity in gene flow (PSC.hetero) in both the “exome capture” and “rna-seq” mussel datasets, and simulations showed that this model has a high accuracy rate, and one of the lowest ambiguity rates (Table S3). Although its relative posterior probability was moderately higher than that of the other models with migration (SC and IM, Table S1), we showed that the method has some power to distinguish PSC from these models (Fig. S3A). These results are in agreement with previous ABC-based studies revealing clearly a secondary contact history between the two mussel species (Fraïsse et al., 2014; Roux et al., 2014), although they relied on eight nuclear loci and on a different set of summary statistics (those that we called “mscalc” here). The periodic connectivity models (PSC and PAM) were not included in these previous studies because of the lack of power to test for intermittent gene exchange since secondary contact. In the present study, we directly compared the use of the jSFS vs. “mscalc” statistics by performing additional ABC inferences on the “mscalc” statistics presented in Table 2. Results were similar to those obtained with jsfs = 23 (Table S1), that is, the best supported models included one or two SCs, except for two datasets (“exome capture” data with n = 2 and n = 8 individuals) for which SI was chosen with a posterior probability of 37% and 34%, respectively. However, the goodness-of-fit of the SI model was quite poor for the standard deviation of the number of fixed sites, “sf_std,” and very poor for the standard deviation of the FST, “fst_std” (“exome capture” data with n = 2 individuals, Fig. S2D). Both statistics were underestimated by the model suggesting that a history without gene flow cannot produce the observed variation of genetic divergence along the genome. Across models, inferences based on the “mscalc” statistics showed similar accuracy rate to those based on jsfs = 23; however, the ambiguity rates for models with current migration were somewhat higher (PSC = 0.52 vs. 0.30, SC = 0.61 vs. 0.31, IM = 0.40 vs. 0.25, respectively). These supplementary analyses suggest that extracting summary statistics from the jSFS can lead to a substantial loss of information.

Conclusion

In this work, we have shown that two high-throughput sequencing datasets (“exome capture” and “rna-seq”), imply the same history of divergence in mussels, regardless of the number of individuals or SNPs sampled, but conditional on the inclusion of informative classes in the jSFS. Thus, genome-wide data coupled with flexible inference methods allow us to test for more complex scenarios of divergence, by providing a comprehensive picture of the gene histories across the genome. Here, we incorporate in our ABC framework heterogeneity in migration rates among loci, to account for the semi-permeability of the barrier to gene flow between recently diverged species (Barton & Bengtsson, 1986). This variation in the rate of effective migration results in variation of genetic divergence along the genome. As shown elsewhere (Roux et al., 2013, 2014, 2016), failing to account for this heterogeneity can mislead inferences. In a similar way, the effect of background selection (i.e., purifying selection at linked loci) and genetic hitchhiking (i.e., positive selection at linked loci) in regions of low recombination can now be incorporated in demographic inferences by including heterogeneity in effective sizes among loci (Sousa & Hey, 2013, Roux et al., 2016, Aeschbacher et al., 2017). Combining modelling approaches that account for both sources of genomic heterogeneity (e.g., Roux et al., 2016) may provide further insight into the complex interplay between linked selection and resistance to introgression during speciation with gene flow (e.g., Duranton et al., 2018).

Supplemental Information