Can secondary contact following range expansion be distinguished from barriers to gene flow?

Introduction

Hybrid zones are narrow regions in which genetically distinct populations meet, mate, and produce hybrids (Barton & Hewitt, 1985). Hybrid zones induced by secondary contact have often been observed in connection with the Quaternary glaciations (Hewitt, 2000). For instance, molecular markers suggest that the southern peninsulas of Europe were major ice age refugia of the European biota and that secondary contact occurred during the northward expansion which followed the last glacial maximum (Taberlet et al., 1998; Hewitt, 1999). With the advent of multi-locus molecular markers such as microsatellite or SNP data, hybrid zones can be investigated using various molecular signatures including gradients of allele frequency, admixture clines, and local increase of genetic differentiation (Nielsen et al., 2003; Adams, Lindmeier & Duvernell, 2006; Strand et al., 2012; Bermond et al., 2012). Molecular or morphological clinal patterns provide evidence for secondary contact in various plant and animal species such as Arabidopsis thaliana (Huber et al., 2014), Silene vulgaris (Keller & Taylor, 2010), the grasshopper Oedaleus decorus (Kindler, Arlettaz & Heckel, 2012), the European hare Lepus europaeus (Antoniou et al., 2013) or the parrotbill bird Paradoxornis webbianus (Qu et al., 2012) to name just a few examples.

However, typical molecular signatures of secondary contact zones can also occur under other evolutionary scenarios. For instance, admixture clines can be observed under pure isolation-by-distance models where nearby populations are connected through gene flow (Engelhardt & Stephens, 2010). Additionally, an increase of genetic differentiation can occur in isolation-by-distance models when there are barriers to dispersal (Barton & Bengtsson, 1986; Riley et al., 2006). With the advent of landscape genetics, the search for barriers to gene flow has attracted considerable attention (Manel et al., 2003; Storfer et al., 2010). Although secondary contact zones can occur at barriers to gene flow, the two models convey different evolutionary paradigms. Models of barriers to gene flow are usually based on isolation-by-distance settings where neighboring populations are connected through dispersal (Nagylaki, 1988; Safner et al., 2011; Blair et al., 2012). Around the barrier to gene flow, dispersal is lowered because of geographical or anthropogenic obstacles (Riley et al., 2006; Zalewski et al., 2009). By contrast, models of secondary contact include an initial phase of evolutionary divergence between two populations or between two sets of populations. The phase of evolutionary divergence is followed by a phase of gene flow between the two divergent units at the secondary contact zone (Murray & Hare, 2006; Durand et al., 2009). The fact that patterns of genetic differentiation can be attributed to different demographic factors is a recurrent problem when using molecular markers. Genetic structure may represent past or contemporary processes and it is notoriously difficult to disentangle between the two possible explanations (Epps & Keyghobadi, 2015).

Here, we use coalescent simulations to investigate to what extent molecular data provide information to distinguish secondary contact following range expansion from barriers to gene flow. We consider neutral simulations only and the barrier to gene flow is modeled by a reduced migration rate. Models of secondary contact of neutral markers are non-equilibrium models that converge to a migration-drift equilibrium (Endler, 1977; Bierne, Gagnaire & David, 2013), whereas locally adaptive loci would rather converge to a selection-migration balance (Barton, 1979). For both evolutionary scenarios, we simulate a one-dimensional stepping-stone model as shown in Figs. 1 and 2.

For comparing the molecular signal left by the two distinct scenarios, we consider statistical measures that have been developed to provide evidence for different demographic processes. The first set of summary statistics, which is used to detect hybrid zones, contains measures of individual admixture coefficients between the parental source populations (Nielsen et al., 2003; Durand et al., 2009), and a measure of linkage disequilibrium (LD) as an increase of LD is expected in admixed populations (McVean, 2002).

The second set of summary statistics contains measures to detect range expansion because secondary contact is frequently induced by geographical expansions of the ancestral populations (Hewitt, 2000). We consider the directionality index ψ as it is sensitive to the occurrence of recent range expansion and it should distinguish between equilibrium and non-equilibrium processes (Peter & Slatkin, 2013). The properties of the directionality index have not been studied yet when there are introgressive events. Furthermore, we include genetic diversity which has been shown to decrease along the expansion direction (Austerlitz et al., 1997).

The last set of summary statistics pertains to isolation-by-distance and barriers to gene flow. First, we include the decay of correlation between allele frequencies as a function of distance as it provides evidence of isolation-by-distance (Hardy & Vekemans, 1999). To detect barriers to gene flow, many numerical methods such as “wombling” identify zones of sharp changes in allele frequencies (Barbujani & Sokal, 1990; Manel et al., 2003). Here, we use local F_ST defined as F_ST per unit of spatial distance for this purpose. Local F_ST can be provided with georeferenced data by the software LocalDiff and we expect them to be larger around the barrier to gene flow (Duforet-Frebourg & Blum, 2014). The development of the software LocalDiff questioned the possibility of distinguishing between the two evolutionary scenarios under consideration. When studying patterns of differentiation in various alpine plants with LocalDiff, we found a region of larger local genetic differentiation in the Italian Aosta valley that was shared across alpine plant species. Both secondary contact zone following postglacial expansions or a barrier to gene flow in a equilibrium stepping-stone process were putative explanations (Duforet-Frebourg & Blum, 2014).

Methods

Results

To test our simulations, we first compare our results to theoretical predictions from the diffusion model. Reassuringly, we find a close agreement. Average clines are accurately predicted by neutral diffusion expectations (Fig. S4). In the equilibrium case without a barrier, pairwise F_ST increases approximately linearly with slope 1/(4Nσ²) (Fig. S5), as predicted by Rousset (1997).

With the validity of the simulations confirmed, we compare the secondary contact model with no barrier at the contact zone (m_B = m) to the stepping-stone model with a barrier. We plot the summary statistics either as a function of time since secondary contact or of the intensity of gene flow across the barrier (Fig. 3; see Supplementary Section S5 for the underlying values). For summary statistics computed per deme (genetic diversity Δ, LD) or per pair of neighboring demes (F_ST), we consider the pattern along the whole range of demes. Many important features are captured by the ratio between the values at the barrier or the suture zone, respectively, and the values to the left and right of it (see Figs. S1–S3 for examples of the pattern along the whole range of demes).

First, we consider the average admixture coefficient for the five populations that are located on the left-hand side of the barrier (demes 45–50). For the isolation-by-distance model with a barrier, these five populations are found to be admixed to an extent depending on the barrier permeability: when increasing the barrier permeability, admixture coefficients of individuals on the left-hand side of the barrier approach 50%. As expected, the populations are also found to be admixed in the secondary contact model without barrier (between 35% and 50%) except for the scenario where data is collected just before secondary contact occurs (t_C = 0).

The ratio between LD at the center (demes 49–52) and on both sides of the range (demes 24–27 and demes 74–77; demes closer to the edge of the range are skipped to avoid the edge effect; see Figs. S1–S3) shows that LD is homogeneous along the whole range of demes for different barrier permeabilities in the stepping-stone model with a slight increase at the barrier for low values of m_B/m. However, in the secondary contact model, LD is considerably increased in the secondary contact zone, at the scale of mixed ancestry predicted by diffusion theory, $\pm 2\text{σ}\sqrt{t_C}$ (Figs. S1 and S2). The excess of LD ranges from a more than twofold to an approximately 1.3 fold increase and decreases as time since secondary contact increases.

Apart from random fluctuations, the directionality index ψ is constant for the stepping-stone model with constant migration rate (m_B/m = 1) as well as for old secondary contact (t_C = 5). More recent secondary contact results in a U-shaped pattern. The pairwise statistics ψ_26,i, i = 27,…,50 first decreases as expected when moving away from the origin of the expansion, but increases again toward the location of secondary contact. For the barrier model with a moderate or strong barrier (m_B/m ≥ 0.1), ψ_26,i remains constant for most of the range, but decreases slightly close to the barrier.

In the stepping-stone simulations the number of pairwise differences (Δ) stays approximately constant over the range of demes (apart from an edge effect) and is hardly affected by the barrier (see also Fig. S3), as expected under the invariance principle. Conversely, in the secondary contact model Δ increases in the suture zone, again at the spatial scale predicted by diffusion theory (Figs. S1 and S2). Only when secondary contact has not occurred yet (t_C = 0), the statistic Δ captures the effect of range expansion, and decreases when moving away from the origin.

When considering the decay of allele frequency correlation as a function of distance, we find a sharp decrease around the suture zone or around the barrier, respectively. In the barrier model, the correlation decreases linearly with distance and drops sharply at the barrier. In the secondary contact model we observe a more sigmoid shape. For older secondary contact, the sigmoid decay converges toward the linear decay of the equilibrium stepping-stone model predicted by Rousset (1997).

Pairwise F_ST between neighboring demes is increased at the barrier to gene flow; the less permeable the barrier, the larger the ratio of F_ST at the barrier compared to the rest of the range. In the secondary contact model, local F_ST is increased at the center when measured just before secondary contact (t_C = 0), but it remains constant along the range of demes when secondary contact is already established (t_C ≥ 1).

To assess the robustness of these results, we performed simulations of less extreme scenarios. Also with more moderate founder events, a lower expansion speed and higher migration rate between demes, we find that the pattern of LD, genetic diversity and the directionality index remain distinctive summary statistics (Supplementary Section S3). However, we also observe that the footprint of secondary contact is more difficult to detect for very slow expansions or high migration rates between neighboring demes (Figs. S7 and S8). But even with these parameter settings, the directionality index ψ remains a discriminant statistic.

If secondary contact occurs at a barrier to gene flow, the difficulty of detecting the secondary contact from molecular data increases. We consider additional simulations where secondary contact occurs in a region where gene flow is reduced by a factor of 10 (m_B/m = 0.1). In many respects, we see an intermediate pattern between the two previously considered scenarios, yet, genetic diversity and the directionality index ψ still provide evidence for secondary contact and LD is even more increased at the contact zone (Fig. 4).

Smaller subsets of the data are used to study the impact of the number of loci on the variance of the summary statistics. Analyses of smaller genomes with 10 and 20 unlinked sequences instead of 100 are presented in Figs. S9 and S10, respectively. The confidence intervals of all summary statistics are considerably larger, while the spatial pattern of the directionality index and the correlogram is preserved and visible, even through larger random fluctuations. The summary statistic least affected by the data reduction is local F_ST. LD and genetic diversity Δ provide a distinctive pattern for very recent secondary contact (t_C = 1) when the data consists of at least 20 loci, but not for 10 loci.

We also consider other sampling schemes where the amount of sampled data is reduced. When sampling four genomes per deme instead of 20, our results remain unchanged (Fig. S11). When sampling only from every fifth deme instead of every deme, they provide similar results for all but one summary statistic (Fig. S12). When reducing spatial sampling, it becomes more difficult to detect the spatial pattern of the directionality index ψ.

Discussion

We find that admixture coefficients alone do not provide sufficient evidence for secondary contact in the presence of isolation-by-distance. Some other summary statistics of genetic data such as local values of F_ST or decay of correlation with distance were not more informative to identify the occurrence of secondary contact either.

By contrast, both an excess of LD and of genetic diversity at the suture zone are found to be unique signatures of secondary contact and also to be informative about the timing of secondary contact. In our simulations we observe an increase of these statistics in areas with mixed ancestry as predicted by the diffusion model ($\pm 2\text{σ}\sqrt{t_C}$). Several of these findings are well supported by previous theoretical considerations, as outlined above. In particular, admixture LD after secondary contact has been thoroughly studied before (Sedghifar et al., 2015). We do not detect such a marked increase of LD near the barrier to gene flow in our stepping-stone model simulations, even for very strong barriers. It seems that the increase of differentiation across the barrier (which increases admixture LD) is roughly balanced by reduced migration across the barrier (which reduces admixture LD). This effect is not well documented in the literature, and warrants further theoretical investigation.

Although a peak of diversity can occur in glacial refugia (Hewitt, 2000), it has been observed before that the genetically most diverse populations were not located in southern Europe but at intermediate latitudes resulting from the admixture of divergent lineages that had expanded from separate refugia (Petit et al., 2003). Many hybrid zones have been formed by such postglacial secondary contact. In some cases, the divergence between hybridizing populations has evolved over very long timescales, and may have arisen already in primary contact. We still expect increased LD and admixture near the suture zones as unique signals of recent secondary contact then, since these signals are generated by allele frequency differences immediately before secondary contact. They do not depend on the evolutionary force that created this divergence.

We expect that excess of LD and of genetic diversity due to admixture would be similarly observed in a two-dimensional model. Qualitative patterns of these signals will be affected by the movement along the second dimension, but the underlying theoretical considerations due to admixture LD and increased diversity due to deeper coalescence time remain the same. The edge effects observed here, that is, reduced diversity and excess of LD because of shorter coalescence times (Wilkins & Wakeley, 2002), will be most pronounced at the four corners of the species range then (Wilkins, 2004).

We also find that the directionality index ψ conveys a signature of secondary contact following expansion. Under range expansion, the ψ statistic is a monotonous function of the distance from the origin of the expansion (Peter & Slatkin, 2013) (see also Fig. 3 for t_C = 0). When secondary contact follows range expansion, it has a distinctive U-shaped pattern (Fig. 3). The distinctive U-shaped pattern is found only for recent enough secondary contact (t_C = 1) but is robust to a wide range of bottleneck intensity and expansion speed values (Figs. S6 and S7). This statistic adds to the toolbox of population geneticists and provides a promising attempt to distinguish between equilibrium and non-equilibrium spatial processes.

The simulation setting was designed to mimic the evolutionary history of species that have undergone a population split during the Quaternary glaciations with subsequent expansion and secondary contact. Assuming a generation time of 1 year and 1,000 diploid organisms per deme, it includes the time frame of expansion and secondary contact after the last glacial maximum in Europe. Species that had spent the last glacial period in southern refugia started to expand northward around 16,000 years ago, and subsequently, many plants established a stable distribution around 6,000 years ago (Hewitt, 1999). We assume the ancestral population split up and started diverging 38,000 years ago (t_S = 19) and the onset of the expansion varies from 16,000 to 6,000 years ago and lasted 6,000 years. Finally, secondary contact is established on the range of 10,000 years ago (t_C = 5) to present time (t_C = 0; in this setting, both populations have expanded, but no gene-flow has occurred yet). Our simulations show that the molecular signal of secondary contact vanishes after approximately 10,000 years. To apply our results to a specific organism, parameters like effective population size, time of divergence and expansion rate need to be calibrated.

In additional simulations, we found that the same summary statistics are distinctive for a wide range of parameter values, and also in datasets with fewer samples. However, reducing the number of loci increases the variance of the summary statistics considerably as it has been observed in coalescent models of single populations or isolation-migration models (Felsenstein, 2006; Wang & Hey, 2010).

Our findings are relevant when investigating modes of speciation using computational approaches (Becquet & Przeworski, 2009). Secondary contact following divergence without gene flow (allopatry) is often compared to models of speciation where species start diverging while exchanging migrants (sympatry or parapatry) (Becquet & Przeworski, 2009; Duvaux et al., 2011; Roux et al., 2013). The different frameworks to study speciation are based on isolation-and-migration models, which do not account for the spatial and potentially continuous repartition of individuals (Pinho & Hey, 2010). As shown in our simulation study, accounting for spatial processes provides additional information that can partly be caught with the ψ directionality index, which has power to reveal evolutionary events such as secondary contact and range expansions.

The fact that isolation-by-distance affects the ascertainment of population structure has already been documented (Novembre & Stephens, 2008; Frantz et al., 2009). Accounting for space is a general recommendation that also stands when studying admixture between divergent populations of the same species (Patterson et al., 2012). Although isolation-by-distance is usually perceived as a confounding factor (Meirmans, 2012), the spatial sampling of individuals is in fact a chance to develop more powerful statistical approaches in evolutionary biology. Accounting for continuous populations should also be possible when performing simulations to choose the most probable scenario of speciation (Duvaux et al., 2011). Numerical simulators of genetic variation that account for the spatial repartitions of individuals are now available (Ray et al., 2010; Kelleher, Barton & Etheridge, 2013). We hope that these developments will encourage researchers to study speciation models that reflect the complex spatio-temporal dynamics of realistic species’ evolutionary histories (Alvarado-Serrano & Hickerson, 2015).

Supplemental Information