Reconstructing the demographic history of divergence between European river and brook lampreys using approximate Bayesian computations

Sampling and genotyping

L. fluviatilis and L. planeri samples were collected from 2010 to 2014 in 6 population pairs from northern France (data from Rougemont et al., 2015 and Fig. 1). Three pairs coexist in sympatry (Aa, Bethune and Oir Rivers) and 2 pairs are not strictly sympatric because a small obstacle occurs on the first river (Hem), while in the second case, populations are located 8 km apart, on the same stream section (Bresle). The last pair is parapatric (Risle), but with a moderate F_ST value similar to what is observed in sympatric populations (Rougemont et al., 2015). We focused on pairs highly connected by gene flow previously identified in Rougemont et al., study because in disconnected pairs, inferences on the long term history of divergence may be biased by the recent effect of genetic drift in isolated L. planeri populations. In parapatric situations L. planeri populations were generally geographically isolated in upper parts of rivers due to anthropogenic barriers to migration with no opportunity for gene flow with L. fluviatilis. The sampling included temporal replicates on the Oir (2010, 2011 and 2014), Bresle (2011 and 2014), and Risle rivers (2011 and 2014). A set of 13 microsatellites was used to genotype a total of 727 individuals following the protocol described in Gaigher et al. (2013).

Summary statistics

Given the lack of genetic differentiation between samples collected in different years in the same river, they were merged together (Rougemont et al., 2015). Similarly, we pooled brook lamprey individuals sampled in upstream and downstream areas on the Aa and Hem river as these two groups were not significantly differentiated (Table S1). We also pooled river lamprey individuals from the Aa and Hem rivers given the lack of differentiation between these populations (Rougemont et al., 2015). Summary statistics were then computed for each pooled sample. As summary statistics used for comparison between simulated and observed datasets, we computed the average and standard deviation values of: the number of alleles (A), Allelic richness (Ar), observed and expected heterozygosity (Ho and He), allele size in base pairs, the Garza-Williamson index (GW, Garza & Williamson, 2001), G_ST (Nei1, 973) and δμ² (Goldstein et al., 1995). The GW index (M ratio) aims at assessing reductions in effective population size and is calculated as $M=\frac{A}{R+1}$ with A is the number of alleles for a given loci in each population and R is the allelic range. Following (Excoffier, Estoup & Cornuet, 2005) we add 1 to the denominator to avoided dividing by zero in case of a monomorphic sample. Goldstein et al. δμ² is a measure of genetic distance developed for microsatellite data defined as δμ2 = (μ₁ − μ₂)² where μ₁ and μ₂ represents the average number of allelic size differences within populations 1 and populations 2 respectively. Each statistics was computed within populations as well as globally except for the Gst and δμ² which are pairwise statistics. All statistics were computed using R scripts (R Development Core Team, 2011) available on github at https://github.com/QuentinRougemont/MicrosatDemogInference.

Testing alternatives demographic scenario

Model selection

Random forest model selection and cross-validations

In parallel to our ABC based model selection and cross-validation procedure we explored the ability of a random forest algorithm (Breiman, 2001) to discriminate the different models and to estimate which summary statistics were the most informative. Random forest (RF) is a machine-learning algorithm whose use has recently been advocated for model choice in ABC inference to circumvent curse of dimensionality problems and those linked to the choice of summary statistics (Pudlo et al., 2014). This approach is a non-parametric classification algorithm that uses bootstrapped decision trees to perform classification using a set (p) of defined predictor variables (here the summary statistics). Multiple (i.e., hundreds to thousands) decision trees are grown and merged together and the ensemble makes up the forest (Breiman, 2001). Simulations that are not used in tree building at each bootstrap (the so called out-of-bag simulation OOB) are then used to compute the OOB error rate, which provides a direct method for cross-validation (Breiman, 2001; Cutler et al., 2007). This method allows reducing the dimensionality of the data (Cutler et al., 2007) but also estimating the relative importance of variables (here the summary statistics) through rankings. Variable importance is measured by random permutations of the specified variable in OOB observations and new predictions are then obtained and compared to the original OOB data (Cutler et al., 2007). One particularly attracting feature of random forest is its insensitivity to strong correlations and high noise within data (Pudlo et al., 2014). In addition, the RF analysis has the advantage of being computationally far less extensive than the ABC cross validation procedure.

We first constructed 6 random forests (one by river) using the randomForestSRC package in R (Ishwaran et al., 2008; Ishwaran & Kogalur, 2015) allowing for parallelization and fast computations. We grew 1,000 trees on subsets of 50,000 simulated datasets (5%) that were used as a training set. Prior analysis using different numbers of trees and training set sizes indicated that the OOB errors reached stationarity using between 500–1,000 trees (see also Fig. 3), so we did not grow bigger forests that would have required extensive computations. All summary statistics were included to get an estimation of the importance of each variable. This allowed us to estimate the OOB error rate for each comparison, which is similar to a prior error rate in ABC inference (Pudlo et al., 2014). Ultimately our forest was used as a prediction tool to compute the probability that our observed data belongs to one of the 5 alternatives models.

Parameter estimation and cross-validation

Parameter estimation was performed for the best models using nonlinear regressions. We first used a logit transformation of the parameters on the 2,000 best replicate simulations providing the smallest Euclidian distance δ (Csilléry, François & Blum, 2012). We then jointly estimated parameters’ posterior probability using the neural network procedure implemented in the abc package. We obtained the best model by weighted nonlinear regressions of the parameters on the summary statistics using 50 feed-forward neural networks and 15 hidden layers. We performed posterior predictive checks for cross-validation in an attempt to check the ability of our parameter estimates to generate data summary statistics close to the observed summary statistics. For each best model, we selected 10,000 samples drawn from the posterior distribution obtained after parameter estimation (from the abc package) and simulated 10,000 new datasets by using again ms and custom R scripts. We then again plotted the distance between our observed original values and our new simulations and computed the p-value for each statistic as the proportion of values that were larger or smaller than the observed value.

Population diversity and divergence

A total of 6 populations pairs (727 individuals in total) were analysed using 13 microsatellites markers. As already observed (Rougemont et al., 2015), the averaged allelic richness of river lamprey populations was significantly higher than that of brook lampreys (Ar_Lf = 3.43, Ar_Lp = 3.116, p = 0.0010, Table 2). However, there was no significant difference in expected heterozygosity between river lamprey (He _Lf = 0.507) and brook lamprey (He _Lp = 0.46). (p = 0.208). The global genetic differentiation between river and brook lampreys was F_ST = 0.061 (99%CI [0.044–0.079]) and pairwise F_ST ranged from 0 to 0.192 (Table S1). The differentiation among river lamprey populations was significantly lower (F_ST = 0.002) than among brook lamprey populations (F_ST = 0.109) (p = 0.003). No river lamprey population differed significantly from the others whereas all brook lamprey populations were significantly differentiated (Table S1). The pairwise differentiation between river lamprey and brook lamprey within each river ranged from 0.028 (Oir and Bethune) to 0.091 in the Bresle River and was significant in all cases but the Bethune River (Table S1).

Model comparison and misclassification

The classical ABC model-choice and random forest approaches generally yielded similar results as detailed in Table 3. In all population pairs, the model of strict isolation (SI) and of ancient migration followed by a period of strict isolation (AM) were clearly rejected. In two population pairs (Aa and Bresle), the best supported model by both model-choice approaches was the SC model. In the Bethune River, the best supported model by both methods was the IM model. In two cases (Hem and Risle), none of the methods was able to accurately discriminate between the two scenarios (SC and IM model). Finally, in the Oir River the two methods gave incongruent results with the model of panmixia (PAN) being the best supported model under the ABC framework, while the RF failed to distinguish between the IM and SC models.

Cross-validations using the ABC framework indicated that the SI and PAN models were correctly classified in all rivers. The PAN model had very low type II errors across all rivers (Table S2) as well as low type I errors in pairwise model comparisons (Table S3). Type II errors for SI and AM comparisons were high (Table S2 and Table S3) but this was a minor concern here as these models were not supported by our data. Cross-validations comparing the two most probable models (IM and SC) indicated a high risk of selecting the IM model instead of the SC model in all rivers (averaged type II error = 0.4025). In contrast, the risk of selecting the SC model instead of the IM was low (averaged type II error = 0.057).

To gain further insights into the performance of our model selection procedure we tested if a random forest approach could help confirming the validity of the empirically rejected models and distinguishing between the IM and SC models. The RF results confirmed that the models of strict isolation, ancient migration and panmixia were classified with a high accuracy (Table 4, Fig. 3 and Fig. S1). The overall error rate (28.79%) hid very different accuracies depending on models. Average OOB errors in pairwise analyses of IM versus SC models were as high as 45% demonstrating that it was generally not possible to correctly classify simulated data in their correct categories (see details in Table 4 and Fig. S1). The estimation of variable importance (Fig. 3 and Fig. S1) indicated that the most informative variables were systematically the mean and variance of G_ST and of δμ², generally followed by the estimators of allelic richness and expected heterozygosity in each population and globally (Fig. 3, Fig. S1 and Table S4).

Difficulty in distinguishing between ongoing migration and secondary contact

In spite of the availability of large amount of genetic data and computer resources, few studies have explicitly tested alternative models of divergence (e.g., Ross-Ibarra, Tenaillon & Gaut, 2009; Duvaux et al., 2011; Roux et al., 2013; Roux et al., 2014; Butlin et al., 2014). While populations may diverge (and eventually become reproductively isolated) under various demographic scenarios, our results indicate that distinguishing between primary differentiation (divergence with gene-flow) versus allopatric divergence followed by secondary contact remains difficult when using a limited number of neutral markers, even with advanced computational tools. Indeed, ABC as well as RF cross-validation clearly showed that the two models were wrongly classified almost half the time. The SC model tended to display a larger proportion of simulations wrongly classified into the IM model, a result that can be explained by the higher complexity of this model that displays one supplementary parameter and is inherently more difficult to infer. In contrast, even though the OOB error rate was high in the IM model, it tended to display fewer simulations wrongly classified into the SC model. Given the inherent difficulty to correctly classify the SC model even when it is true, our support for this model in some cases seems conservative and may suggest that it could be the true model under which lampreys have diverged.

Our inability to distinguish between a scenario of isolation with migration and secondary contact is in accordance with theoretical expectations from Bierne, Gagnaire & David (2013). Using a simple modelling approach they showed how genetic environmental associations at neutral markers such as microsatellites can be quickly lost after secondary contacts and then reach migration/drift equilibrium together with a pattern of isolation by distance, which is the pattern observed in the studied lamprey populations (see Rougemont et al., 2015). This was indeed expected under the relatively high mutation rate of microsatellite loci that allows neutral equilibrium to be rapidly attained. Bierne, Gagnaire & David (2013) applied their model to the well-studied freshwater/marine stickleback system (e.g., Colosimo et al., 2005; Hohenlohe et al., 2010; Hohenlohe et al., 2012), which shares several characteristics with the lamprey system: a single nearly panmictic marine population and small and potentially independent freshwater populations. The application to the stickleback model showed that introgression proceeded independently between the different streams and was strongly asymmetric from the migratory to the resident populations, which is similar to the pattern we observed here (Table 5) under both the isolation with migration model and the secondary contact model. Although we cannot accurately infer the symmetry of migration, the migration appeared slightly higher from river lamprey to brook lamprey in some rivers (Table 5).

The difficulty in distinguishing the two models is further explained by the fact that when allopatric divergence is short and secondary contacts are very long, the SC model converges to the same signal as that obtained under the IM model. The same difficulty is expected when comparing a model of ancient migration in which migration lasts a short amount of time and is followed by a long period of divergence without gene flow to a model of strict isolation. Accordingly, these models could not be clearly distinguished in our cross-validations (e.g., Table 4).

The failure to reject panmixia in the Oir River can also be explained in the light of Bierne, Gagnaire & David (2013) conclusions. It could be linked to the low genetic divergence observed (F_ST = 0.028) especially given the small number of markers we used, but this pattern of nearly panmixia can also be attributed to a stronger introgression in this system than in all other investigated streams. In this river the mean size of river lampreys (225 mm, n = 134) was much smaller than the size observed in other sites (mean = 303 mm, n = 389). Assuming that size difference is the most important cause of reproductive isolation (Beamish & Neville, 1992) a smaller size difference may facilitate mating of the two ecotypes and subsequent genome swamping. In both cases, inferences from this system based on neutral markers are necessarily difficult as this pattern may be explained by strong gene flow after an isolation period as well as by a very early stage of ongoing divergence.

Demographic parameter estimations and new insights on lamprey history

We observed asymmetries in effective populations size between river lamprey (averaged median = 1,448; 95% CI [890–2,300] under SC, 1,158 and 95% CI [788–1,710] under IM) and brook lamprey (averaged median = 612; 95% CI [356–1,090]; and 494 95% CI [318–840] under SC and IM respectively). The greater estimates of effective population size in river lamprey reflect the idea that homing is moderate in this species and population size as well as gene flow are large (Rougemont et al., 2015; Spice et al., 2012) resulting in a situation similar to that observed in the panmictic populations of marine sticklebacks (Hohenlohe et al., 2010; Hohenlohe et al., 2012). We also observed reductions in effective population size of each resident brook lamprey population as compared to the ancestral populations (averaged median = 1,428; 95% CI [962–1,886]; and median = 1,748 with 95% CI [1,490–2,026] under SC and IM respectively). Such a reduction is expected following the independent river colonization by founding individuals from the resident ecotype (Bierne, Gagnaire & David, 2013).

Estimates of divergence times suggested that the two ecotypes may have split around 201,000 years ago (95% CI [148,000–250,000]) under the IM model and around 257,000 years ago (95% CI [195,000–344,000]) under the SC model (assuming a generation time of 5 years, Potter & Potter (1971)). Such estimates are rather similar to what was observed in Dicentrarchus labrax (Tine et al., 2014). Consequently, it seems unlikely that the divergence was initiated rapidly, following the recent glacial retreats around 10,000–15,000 years ago (Bernatchez & Wilson, 1998; Aldenhoven et al., 2010) as often proposed to explained ecotypic divergence in various aquatic species (e.g., Schluter & Nagel, 1995; Espanhol, Almeida & Alves, 2007; Bracken et al., 2015; Mateus et al., 2016). Importantly, under the SC model, the secondary contact would have started around 92,000 years ago, representing more than one third of the total divergence time between the two species. Such an ancient secondary contact implies that the genetic signature of historical geographic isolation carried by neutral markers may have been lost. In these conditions, neutral markers can converge to the same state than the one observed under primary differentiation explaining again the difficulty of discriminating the two models (Barton & Hewitt, 1985; Charlesworth, Nordborg & Charlesworth, 1997; Bierne, Gagnaire & David, 2013). The SC model implies the accumulation of some Dobzhansky-Muller incompatibilities in allopatry when the two ecotypes started to diverge. While both theory (Orr, 1995) and empirical evidence (Moyle & Nakazato, 2009; Matute et al., 2010; Wang, White & Payseur, 2015) predict that DMI should accumulate faster than linearly in time, our result suggest a limited amount of isolation. The period of isolation was certainly too short to allow for sufficient DMI to accumulate and to develop strong barriers to gene flow. This would likely explain the low differentiation observed for mtDNA (Espanhol, Almeida & Alves, 2007; Blank, Jürss & Bastrop, 2008; Bracken et al., 2015) and is fully compatible with the observation of viable F1 hybrids (Hume et al., 2013; Rougemont et al., 2015).

Another question we started to address is the origin of parallel divergence observed between river lamprey and brook lamprey. The process of recent and independent postglacial divergence with gene flow (parallel divergence) either from standing genetic variation or de novo mutations seems a conceptually simple explanation often used to explain divergence between freshwater resident and anadromous (or marine) stickleback (Schluter & Nagel, 1995). Similar scenarios have been proposed in lampreys (e.g., Espanhol, Almeida & Alves, 2007; Bracken et al., 2015; Mateus et al., 2016). However, models of divergence with gene flow did not receive higher support than models of secondary contact and divergence time estimates appeared older than the onset of glacial retreats. Alternative scenarios involve either multiple independent secondary contacts between populations inhabiting different refugia or secondary contact between anadromous parasitic and freshwater resident populations each having diverged in a different ancestral place. Under this scenario a spatial re-assortment of alleles involved in non-parasitism would allow the recolonization of neighboring rivers (i.e., the transporter hypothesis; Schluter & Conte, 2009; Bierne, Gagnaire & David, 2013; Welch & Jiggins, 2014). Considering (i) the relatively small scale of investigation and (ii) the nearly panmictic situation in river lamprey; the hypothesis of a spatial re-assortment of ancestral variation by migration between adjacent rivers appears more parsimonious. In this case, the history of divergence would reflect a re-interpretation of the transport hypothesis under a scenario of secondary differentiation rather than primary differentiation. The major conceptual difference is the existence of the brook lamprey background before the recent colonization of rivers. This latter process may have arisen either through hybrid genotypes colonizing rivers or transport of alleles broken up by recombination and at low frequency in the river lamprey background. While the global divergence implies secondary contact, it is debatable whether the spatial re-assortment process constitutes divergence with gene flow or secondary contact hence the difficulty in distinguishing between the two models. Disentangling these scenarios of primary versus secondary differentiation under the transporter hypothesis is challenging and may be addressed by combining genome wide data with historical modelling among multiple pairs of populations. Our results also suggest that the RF method provides a valuable complement to the standard ABC model comparison (Robert et al., 2011; Pudlo et al., 2014; Marin et al., 2014). The two methods provided similar outcomes in terms of model choice and subsequent cross-validations except in one population pair (Oir River). Our ability to distinguish between SC and IM was low in both cases. In line with Pudlo et al. (2014) we find that the RF approach possesses a series of advantages over the ABC approach such as (1) fast model choice procedure with simultaneous cross validation through OOB computations, and (2) considerable reduction of computational time. Estimating variable importance can be particularly interesting when a large set of variables are used without prior knowledge about the pertinence of the summary statistics used. The choice of summary statistics is an important process in ABC methodology (Csilléry et al., 2010). RF may provide such an objective tool that may be complementary to conventional ABC model choice and cross validation procedure. Note however, that the neural network method provided in the abc package, performed very well and provided similar results to the RF model.