The Spatial Signature of Introgression After a Biological Invasion With Hybridization

Quilodrán, Claudio S.; Tsoupas, Alexandros; Currat, Mathias

doi:10.3389/fevo.2020.569620

ORIGINAL RESEARCH article

Front. Ecol. Evol., 29 September 2020

Sec. Population, Community, and Ecosystem Dynamics

Volume 8 - 2020 | https://doi.org/10.3389/fevo.2020.569620

This article is part of the Research Topic The Genomics of Biological Invasion View all 8 articles

The Spatial Signature of Introgression After a Biological Invasion With Hybridization

Updated

A correction has been applied to this article in:

Corrigendum: The Spatial Signature of Introgression After a Biological Invasion With Hybridization
1. Read correction

$\r\nClaudio S. Quilodrn,*$ Claudio S. Quilodrán^1,2*

Alexandros Tsoupas²

Mathias Currat^2,3

¹Department of Zoology, University of Oxford, Oxford, United Kingdom
²Department of Genetics and Evolution – Anthropology Unit, University of Geneva, Geneva, Switzerland
³Institute of Genetics and Genomics in Geneva (IGE3), University of Geneva, Geneva, Switzerland

The accumulation of genome-wide molecular data has emphasized the important role of hybridization in the evolution of many organisms, which may carry introgressed genomic segments resulting from past admixture events with other taxa. Despite a number of examples of hybridization occurring during biological invasions, the resulting spatial patterns of genomic introgression remain poorly understood. Preliminary simulation studies have suggested a heterogeneous spatial level of introgression for invasive taxa after range expansion. We investigated in detail the robustness of this pattern and its persistence over time for both invasive and local organisms. Using spatially explicit simulations, we explored the spatial distribution of introgression across the area of colonization of an invasive taxon hybridizing with a local taxon. The general pattern for neutral loci supported by our results is an increasing introgression of local genes into the invasive taxon with the increase in the distance from the source of the invasion and a decreasing introgression of invasive genes into the local taxon. However, we also show there is some variation in this general trend depending on the scenario investigated. Spatial heterogeneity of introgression within a given taxon is thus an expected neutral pattern in structured populations after a biological invasion with a low to moderate amount of hybridization. We further show that this pattern is consistent with published empirical observations. Using additional simulations, we argue that the spatial pattern of Neanderthal introgression in modern humans, which has been documented to be higher in Asia than in Europe, can be explained by a model of hybridization with Neanderthals in Eurasia during the range expansion of modern humans from Africa. Our results support the view that weak hybridization during range expansion may explain spatially heterogeneous introgression patterns without the need to invoke selection.

Introduction

The evolutionary history of many species has been influenced by range expansion and hybridization with related taxa (e.g., Nielsen et al., 2017; Duvernell et al., 2019). In addition, the range expansion of invasive species is a worldwide conservation issue due to human translocation of species and climate change, which modify migration patterns and the time and place of reproduction of numerous species (e.g., Muhlfeld et al., 2014; Larson et al., 2019).

Hybridization is recognized as a major factor influencing evolutionary change. The gene flow between taxa, leading to introgression, may increase genetic diversity and opportunities for adaptation to changing environmental conditions (Arnold and Martin, 2010; Taylor et al., 2015). There are multiple examples of adaptive introgression (Chhatre et al., 2018), as well as neutral genomic regions introgressed from past or current hybridization events (Schumer et al., 2016). During a biological invasion, a dynamic hybrid zone may arise, which is characterized by changes in the pattern of introgression in space and time (Billerman et al., 2019). Understanding the emergence of dynamic hybrid zones is important for elucidating the mechanisms that maintain divergence between lineages, despite gene flow, under changing environmental conditions (Mallet, 1986; Buggs, 2007; Taylor et al., 2015). This understanding is also important because hybridization may increase the extinction risk of species, for instance, due to outbreeding depression (Oakley et al., 2015), the introduction of maladaptive genes (Kidd et al., 2009), wasted reproductive effort (Quilodrán et al., 2014) or genetic swamping (Todesco et al., 2016).

The resulting level of introgression is influenced by the demographic dynamic of range expansion during a biological invasion. A well-documented pattern is the expected asymmetric introgression between local and invasive organisms, with small levels of introgression into the local organisms and considerable introgression in the invasive organisms (Currat et al., 2008; Excoffier et al., 2014; Quilodrán et al., 2019). There are many examples following this pattern in both plant and animal kingdoms (e.g., Duminil et al., 2006; Darling et al., 2014; Oswald et al., 2019). This pattern is influenced by three main factors. The first factor is demic diffusion, in which the invader gene pool is progressively diluted along the expansion due to recurrent admixture with the local organisms (Ammerman and Cavalli-Sforza, 1973). The second factor is different demographies, as the local organisms are at equilibrium while the invaders arrive with few individuals, increasing their likelihood of being involved in heterospecific mating (Hubbs, 1955). The third factor is allelic surfing (Edmonds et al., 2004; Klopfstein et al., 2006), in which genes introgressed at the wave front may reach very high frequencies in the newly colonized areas due to a combination of genetic drift, founder effects and demographic increases (Excoffier and Ray, 2008). A local allele introgressed at the wave front of the invasive range expansion therefore has less chance to be lost by drift when the invader population is rapidly growing (Excoffier et al., 2009). While exceptions appear when admixture rates are so low that hybridization events do not occur during the expansion phase but at demographic equilibrium only, this pattern is still observed when the invasive organism is at a higher density than the local organism or even when there is competition eventually driving the extinction of the local organisms (Currat et al., 2008). The relative densities of the two interacting taxa will influence the rate of admixture required to obtain the asymmetric pattern of introgression. If the size of the local taxon is larger than that of the invasive one, the latter is rapidly introgressed and asymmetry is obtained at a much lower rate of admixture. Conversely, if the invasive taxon’s density is larger than the local taxon, the invasive taxon will be less introgressed and more admixture will be required to achieve an asymmetric introgression pattern (Currat et al., 2008). In other words, a numerical advantage for the invasive taxon counterbalances, but does not completely prevent, the introgression of local genes through interspecific gene flow. These results hold for neutral genetic diversity and could be different when selection is influencing the system.

This general introgression pattern of asymmetry between local and invasive taxa in range expansion was first described by Currat et al. (2008). In addition, this study noticed the geographic heterogeneity of introgression levels, with introgression of local genes increasing in the invasive taxon with the increase in the distance from the source of the invasion. The existence of introgression gradients from the origin of a range expansion toward its margin was previously noted by Currat and Excoffier (2004, 2005). These gradients are often weak and occur for only a narrow range of low or intermediate interbreeding rates; otherwise, the invasive taxon is mostly introgressed by local genes, and no gradient is visible. Those pioneering studies focused on the degree of admixture between taxa that primarily affect the amount of introgression rather than the spatial patterns of introgression. Moreover, they did not evaluate in detail the spatial heterogeneity of introgression into the local taxa. Evaluating the potential effect of hybridization in generating spatial introgression gradients is important for elucidating the evolutionary history of the colonization of invasive organisms, as well as for assisting conservation management efforts by accurately projecting potential areas of future microevolutionary processes.

The aim of our study is therefore to explore in more detail the spatial distribution of introgression in both local and invasive taxa after a range expansion of the latter and the persistence of introgression patterns over time. In addition to the scenario originally investigated in Currat et al. (2008), which represents interactions between a local and invasive taxon throughout the whole range of the biological invasion, we also investigate a scenario where the local taxon occupies only a small portion of the area colonized by the invasive taxon. This scenario is a representation of real situations where an expanding population could have carried genes introgressed from a local population to areas where the local population never existed, such as Neanderthal genes carried by modern humans in East Asia (Green et al., 2010; Currat and Excoffier, 2011). We also investigate a third scenario in which two invasive taxa from different areas meet in a previously empty area, a situation that could happen, for instance, during postglacial recolonization (e.g., Hewitt, 2001). We compare the spatial patterns obtained under those three different scenarios and their persistence over time. We further evaluated the predictions of our theoretical approach with empirical observations of biological invasions in both plants and animals.

Finally, we use the reported case of hybridization between Neanderthals and modern humans to illustrate the influence of species spatial dynamics on the final level of introgression after a range expansion. Modern humans experienced a range expansion out of Africa and occasionally admixed with Neanderthals when colonizing Europe and Asia approximately 40,000–60,000 years before present (BP) (Currat and Excoffier, 2011; Nielsen et al., 2017; Skoglund and Mathieson, 2018). A level of introgression of approximately 2% has been estimated in current non-sub-Saharan human populations (Green et al., 2010). Moreover, this level of introgression is estimated to be approximately 12–20% higher in modern East Asian individuals than in modern European individuals (Meyer et al., 2012; Wall et al., 2013; Prüfer et al., 2017; Villanea and Schraiber, 2019). While selection has been proposed as an explanation for this spatially heterogeneous introgression pattern (Sankararaman et al., 2014; Villanea and Schraiber, 2019), we explore whether the species spatial dynamics may also explain the different levels of introgression observed in Asia and Europe without invoking selection. We based our simulations on the framework developed by Currat and Excoffier (2011) to estimate the rate of admixture that best explains the observed introgression from Neanderthal into current human populations (see also Excoffier et al., 2014). Both studies reported strong barriers to gene flow between the two taxa that could generate spatial heterogeneity in introgression levels across Eurasia.

Materials and Methods

Simulation Framework

We simulated biological invasions with various levels of interbreeding between two closely related taxa by using the program SPLATCHE3 (Currat et al., 2019). SPLATCHE3 is a spatially explicit program that allows the simulation of neutral genetic diversity. The simulations are performed in two consecutive steps: first, a forward simulation of the demography and migration of the taxa; second, a backward coalescent simulation that reconstructs the genealogy of a series of samples by using the demography and migration recorded during the previous step. The proportion of introgressed genes is measured by using this reconstructed genealogy.

Each simulated taxon is represented by a grid of demes arranged in a stepping-stone manner (Kimura and Weiss, 1964). Population densities are logistically regulated by using the parameters r (growth rate) and K (carrying capacity) in each deme. The gene flow within a grid of demes represents migration between neighbor demes belonging to the same taxon and is regulated by the parameter m (migration rate). The gene flow between the grid of demes represents hybridization between taxa and is regulated by the parameter γ (interbreeding success rate). This parameter represents the strength of the barriers to gene flow between interacting taxa. There is no interbreeding when this parameter is equal to zero, while both taxa behave as a single panmictic population when this value is 1. Mating between taxa is non-random with any intermediate value of γ. Hybridization can only occur between each pair of demes located in the same place in their respective grid. More information about the algorithms implemented in SPLATCHE3 is available in Currat et al. (2019). The executable file and input examples for simulations are available at www.splatche.com/splatche3.

Invasion in a Virtual Square World

Following the framework of Currat et al. (2008), a two-dimensional square space of 100 × 100 demes was used to simulate the range expansion and hybridization between taxa during 2,000 generations. We explored three different main scenarios of biological invasion, as presented in Figure 1 and described below:

FIGURE 1

Figure 1. Schematic representation of simulated scenarios of biological invasion and hybridization between two taxa. The “whole area” scenario (A): one taxon (T1) is already present in the whole area and is considered local, whereas a second taxon (T2) arrives with few individuals and starts an invasive range expansion from a corner of the distribution range of T1. The “restricted area” scenario (B): similar to the previous scenario, but the local taxon (T1) is distributed in a more restricted area (25% of the total area). The “two invasions” scenario (C): both taxa are invasive and start a range expansion from opposite corners of the simulated area. Empty demes are suitable to be colonized by both taxa in this last scenario.

(A) “Whole Area” Scenario

In the first scenario, 50 individuals of an invasive organism start to colonize the grid from the left corner of the virtual area, representing an invasive taxon under range expansion. The second taxon was already present throughout the area at demographic equilibrium and is considered to be local. This local taxon colonized the whole area 1,000 generations before the arrival of the invasive taxon (i.e., interaction between the two taxa starting at generation 1,000).

(B) “Restricted Area” Scenario

The second scenario is similar to the previous one, but the local taxon is restricted to a smaller area than the one suitable to be colonized by the invasive. This restricted area available for the local taxon represents 25% of the total area (25 × 25 demes in the left corner of the virtual world).

(C) “Two Invasions” Scenario

In the third scenario, both taxa are considered invasive organisms colonizing an empty area, each one starting from the opposite corner of the virtual world with an initial population size of 50 individuals. There is no pre-established local taxon under this scenario.

For each main scenario (A, B, and C), we explored different combinations of carrying capacity (K) and number of migrants between demes (km) by using five sub-scenarios nested within each main scenario (Table 1). These sub-scenarios are identical to the non-competition scenarios (NC) explored by Currat et al. (2008). These sub-scenarios are used to directly compare the three main scenarios of biological invasion (A, B, and C). Still following Currat et al. (2008), we explored values of interbreeding success rates (γ) between 0 and 6%. We used this range because under the genetic and demographic conditions explored in this virtual square world, higher values of γ lead to a complete introgression of the invasive gene pool by local genes after the simulated generations, and thus, no spatial heterogeneity is found. We performed 10,000 backward coalescent simulations, representing 10,000 independent neutral loci. For each locus, the proportion of introgression from taxon i to taxon j was computed by tracing back the genealogy of 1,000 gene copies sampled in 25 demes from taxon j that were equally spaced in the virtual world (40 genes per deme). The introgressed lineages represent the ones sampled in taxon j that were in taxon i at the onset of the invasive range expansion, and the reverse (lineages in taxon i that were in taxon j). We finally averaged the introgression proportions over the 10,000 simulated independent loci. The introgression levels were plotted with the function “filled.contour” in R (R Development Core Team, 2019).

TABLE 1

Table 1. Simulated sub-scenarios with different parameter values of demography and migration during the range expansion of an invasive taxon and hybridization with a local taxon (see Figure 1).

Neanderthal Introgression in Modern Humans

We simulated the range expansion of modern humans (Homo sapiens) and their hybridization with Neanderthals (Homo neanderthalensis) by following the same simulation framework used by Currat and Excoffier (2011). In the context of our modeling framework, we consider modern humans to be the invasive taxon and Neanderthals to be the local taxon. Here, we used only the best scenario estimated by Currat and Excoffier (2011), in which hybridization occurred several times over a large period and over a large area in western Eurasia. This scenario is supported by additional paleogenomic studies (Krause et al., 2007; Prüfer et al., 2014; Mafessoni, 2019), even if the number of hybridization events or pulses is under debate (e.g., Sankararaman et al., 2012; Vernot et al., 2016; Villanea and Schraiber, 2019). Under this scenario, modern humans start a range expansion at approximately 50,000 BP from northeastern Africa, while Neanderthals were already present in western Eurasia. This scenario is similar to the theoretical scenario B (“restricted area”) described in the previous section, where the area initially occupied by the local taxon (Neanderthals) is restricted compared to the area colonized by the invasive taxon (modern humans).

The continental surface of North Africa and Eurasia (Hammer–Aitoff projection) was represented by a grid of more than 6,000 demes, each of them representing an area of 100 km × 100 km. The carrying capacity of Neanderthals (K_N) and modern humans (K_H) were set to values representing 0.025 and 0.1 individuals per km², respectively (K_N = 200 and K_H = 800 in each deme). These values are compatible with previous estimations of the prehistorical population densities (Bocquet-Appel and Demars, 2000; Currat and Excoffier, 2004). The population growth rate parameter (r) was equal to 0.8 for both modern humans and Neanderthals. The migration rates were set to 0.2 for modern humans (m_H) and 0.1 for Neanderthals (m_N). Both migration and growth rates were chosen to have colonization times in Europe and Asia in accordance with archeological knowledge (Mellars, 2006).

The population densities in each deme are logically regulated by using a modified version of the Lotka–Volterra competition model (Volterra, 1928; Lotka, 1932). Similar to the classic Lotka–Volterra model, the modified version considers intraspecific and interspecific competition. However, the model used in this study differs in that the competition between species is not a fixed parameter; rather, it is a variable value that depends on the density of both species in each deme (Currat and Excoffier, 2004). The influence of one species over the other therefore increases with its density. This model of competition has already been used in previous simulations of interactions between Neanderthal and modern humans. This model also considers a lower carrying capacity for the Neanderthals than for the modern humans (Currat and Excoffier, 2004, 2011), thus assuming that modern humans have a competitive advantage over Neanderthals on exploiting the local resources (Klein, 2008).

Approximate Bayesian Computation Approach

We used an approximate Bayesian computation approach (ABC) to estimate the interbreeding rate parameter (γ) that best explains the level of introgression from Neanderthals to modern humans (Beaumont et al., 2002). The ABC allows us to estimate a 95% confidence interval associated with the estimated value of γ, which was further used to compare the resulting range of introgression, which is expected to be higher in Asia by 12–20% than that in Europe. Following Currat and Excoffier (2011), we compared the observed with the simulated introgression for one sample located in France and another sample located in China. We used the average introgression level of 2% as the observed value in both France and in China to select the best set of simulations. The higher introgression level estimated in Asia compared to Europe was therefore not included as prior information in our analysis, unlike Excoffier et al. (2014). We used a neural network with a tolerance level of 2% for parameter estimation. The tolerance level is the proportion of simulations retained for the estimation of parameter values, which represents the simulated introgression values that are closest to the observed introgression. This condition means that the best 8,000 out of 400,000 simulations were retained for parameter estimation. The resulting estimations of the interbreeding rate are similar when using tolerance levels of 1 and 5% (Supplementary Table S1, supporting information). This analysis was performed with the R package “abc” 2.1 (Csilléry et al., 2012).

Literature Survey

We performed a literature survey of cases of invasive taxa in range expansion and hybridization with related local taxa. We included examples of hybridization at the intraspecific and interspecific levels. We used the online database ISI Web of Knowledge with the following terms: (A) “range expansion” AND admixture AND hybrid^∗; (B) “range expansion” AND introgression; (C) “increasing introgression” OR “increasing admixture” OR “increasing hybridization”; (D) introgression AND invasion OR invasive; (E) admixture AND invasion OR invasive. We included literature published from January 1995 to March 1, 2020. We searched for studies reporting the spatial distribution of introgression for one or both interacting taxa but also for studies in which the spatial distribution of admixture proportions was present. Indeed, the admixture components of the hybrid population were used as a proxy of the level of introgression in many studies instead of proper introgression levels. All studies reporting admixture proportions used the software STRUCTURE (Pritchard et al., 2000). We excluded all cases in which either the spatial distribution of introgression or admixture was not reported. We included examples of any kind of genetic marker, such as autosomal chromosomes, sexual chromosomes, mitochondrial DNA and chloroplasts, for animal and plant taxa, respectively.

Results

Spatial Pattern of Introgression

We explored the resulting spatial pattern of introgression under various scenarios of interactions between two taxa interbreeding with a level of 5% (γ = 0.05) in a virtual square world. After 2,000 generations of interactions, a spatial pattern of introgression was observed in both the invasive and the local taxon (Figure 2). The “whole area” (A) scenario shows that local introgression in the invasive gene pool becomes higher with increasing distance from the source of the biological invasion, while the opposite is found for invasive introgression in the local gene pool (e.g., introgression becomes lower with increasing distance from the source of the invasion). However, the gradient may be so weak that it is hardly visible in graphs showing absolute introgression values instead of relative values (i.e., sub-scenarios NC2 and NC5 in Figure 2 and Supplementary Figure S1). The scenario “restricted area” (B) shows identical patterns in the area where both taxa occur together (bottom left quarter of the area), whereas it shows a different pattern in areas never occupied by the local taxon. Along the main axis of the invasion (diagonal from the source toward the upper-right corner), the introgression of local genes in the invasive gene pool resembles that obtained when the two taxa coexist but not in the periphery (upper-left and lower-right of the area), where introgression is much lower. This result occurs because along the main axis of the expansion, the opportunity to admix is higher (more deme to cross) than that along the shorter x and y axes, where fewer local demes are crossed by the invasive taxon. In the scenario “two invasions” (C), in which both taxa are invasive, the introgression level in each taxon increases with the distance from their respective source of expansion (Figure 2). Because both taxa are colonizing from opposite corners of the simulated area, there is a strong asymmetric introgression between them (Figure 2 and Supplementary Figure S1).

FIGURE 2

Figure 2. Spatial pattern of introgression after biological invasions and hybridization. Two taxa were simulated with an interbreeding success rate of 5% (γ = 0.05). The different colors denote the amount of introgression in each of the interacting taxa (%). Three scenarios of biological invasion (see Figure 1) and five demographic sub-scenarios (see Table 1) were simulated. The demographic parameters of sub-scenarios for taxon 1 and taxon 2 are also presented (taxon 1 | taxon 2). Taxon 1 is considered to be local under the “whole area” and “restricted area” scenarios but invasive under the “two invasions” scenario. Taxon 2 is always considered to be an invasive taxon. The values of introgression are averaged over 10,000 simulations.

The same general trend is observed among all five sub-scenarios of each main scenario, with some variation depending on the population sizes and numbers of migrants (NC1, NC2, NC3, NC4, and NC5 in Figure 2). For instance, when increasing the population densities of both taxa, the spatial gradient is less pronounced than that with smaller population sizes (compare NC1 and NC2). However, the gradient is more pronounced when higher densities send a greater number of migrants to surrounding areas (NC2 vs. NC3). The pattern still occurs, even when both populations have different population sizes, when either the invasive or the local taxon has a higher density than the other (NC4 vs. NC5). In this last sub-scenario, when the local taxon has a much higher density than the invasive taxon and occurs over the whole area, the spatial gradient of introgression is weak in both taxa (Supplementary Figure S1, supporting information). The spatial gradient is more easily observed when plotting the relative difference of introgression for each simulated area (Supplementary Figure S1, Supporting information).

Interbreeding Success Rate

We explored a range of interbreeding success rates to evaluate its influence on the spatial gradient of introgression. We explored the effect of the interbreeding rate on the level of introgression along the main axis of the invasion (diagonal from the source toward the upper-right corner, Figure 3). A minimum γ threshold is needed for the emergence of the spatial gradient of introgression depending on the explored sub-scenario and scenario of biological invasion. A small gradient following the expected trend is observed when plotting the relative differences of introgression along the main axis of invasion, in which case a minimum threshold is also needed (Supplementary Figure S2, supporting information), even where the gradient is extremely weak, as revealed by the absolute values (Figure 3). Note that no gradients, as well as a heterogeneous spatial level of introgression, may be obtained when values of interbreeding rates are very small (γ; Figure 3 and Supplementary Figure S2, supporting information). This finding supports the observation that the spatial gradient of introgression can only be found for a narrow range of interbreeding rates.

FIGURE 3

Figure 3. Effect of the interbreeding rate (γ) and distance along the diagonal axis from the source of the invasion with a variable amount of interbreeding rate. The distance along the diagonal axis is always measured from the bottom left corner. Two taxa were simulated. The different colors denote the amount of introgression in each of the interacting taxa (%). Three scenarios of biological invasion (see Figure 1) and five demographic sub-scenarios (see Table 1) were simulated. The demographic parameters of sub-scenarios for taxon 1 and taxon 2 are also presented (taxon 1 | taxon 2). Taxon 1 is considered to be local under the “whole area” and “restricted area” scenarios but invasive under the “two invasions” scenario. Taxon 2 is always considered to be an invasive taxon. The values of introgression are averaged over 10,000 simulations.

The minimum γ threshold that results in a spatial gradient of introgression is also influenced by population sizes and the number of migrants. Indeed, this threshold is very small, and the gradient is steep when the population sizes of both taxa are large (NC2 vs. NC1), with an intermediate pattern when high-density populations send out a greater number of migrants (NC3). When the invasive taxon has a much higher density than the local taxon, the γ threshold necessary to obtain a gradient of introgression is relatively similar to the case where both taxa have equal and small densities, but the gradient is steeper in the local population than in the invasive population (NC4 vs. NC1), except for the “restricted area” scenario. When the local taxon has a higher density than the invasive taxon, there is almost no gradient under the “whole area” scenario (but see Supplementary Figure S2, supporting information, for relative values). This result is because, in such conditions, very small values of interbreeding result in an almost complete genetic replacement in the invasive taxon and an almost null introgression in the local taxon (Figure 3). However, similar patterns to equal small populations are observed under the “restricted area” scenario and “two invasions” scenario. For the “restricted area” scenario, there are almost no differences among the various sub-scenarios, except when both densities (local and invasive) are large (NC2). In that latter case, the gradient in the invasive taxon decreases at short distances and then increases again at long distances. For the “two invasions” scenario, in which both taxa behave as an invasive taxon in range expansion, there are always gradients of high introgression when the distance from the expansion sources increases.

Temporal Effect

The spatially asymmetric introgression between the source of the biological expansion and the more distant region, in which both taxa are present (extreme demes), appears during the colonization of the invasive taxon (Figure 4). Under all sampled times, the invasive introgression in the local gene pool is inversely proportional to distance, while the opposite is observed for the invasive gene pool (Figure 4). In the third scenario of biological invasion (“two invasions”), both taxa are invasive, and introgression is directly proportional to the distance from the expansion source for each taxon under all sampled times.

FIGURE 4

Figure 4. Effect of time since the onset of biological invasion on the level of introgression. Two interacting taxa are represented with different colors (yellow and gray). The level of introgression (%) is estimated in both the source of the expansion and the extreme corner in which both taxa are present. The arrows indicate the direction of the invasive taxon colonization. Three scenarios of biological invasion (see Figure 1) and five demographic sub-scenarios (see Table 1) were simulated. The demographic parameters of sub-scenarios for taxon 1 and taxon 2 are also presented (taxon 1 | taxon 2). Taxon 1 is considered to be local under the “whole area” and “restricted area” scenarios but invasive under the “two invasions” scenario. Taxon 2 is always considered to be an invasive taxon. The values of introgression are averaged over 10,000 simulations.

The range of variation in the value of introgression along the main axis of the expansion is at its maximum when both taxa are invasive (scenario “two invasions”). This range seems smaller for the “restricted area” scenario than for the “whole area” scenario, but the distance between the two sampled locations is not equal in those two scenarios. We used the demes located as far away from each other as possible along the main axis of expansion (the diagonal linking the bottom left and top right corners of the area) in which both taxa were present, meaning that the distance is shorter for the “restricted area” scenario (B) than for the “whole area” scenario (A). This difference is also the reason for which the invasive taxon seems to colonize the range of the local taxon earlier under the “restricted area” scenario than under other scenarios. The range of variation between extreme demes is similar when both taxa have equal population sizes and migration rates (NC1, NC2, and NC3). However, this variation tends to be larger when the invasive taxon has a higher population size than the local taxon (NC4) and smaller when the local taxon has a higher density (NC5).

Looking at the evolution of introgression levels through time, we found that the general trend is toward a gradual but very slow decrease in the difference in the level of introgression between the two extremities of the axis of expansion (source and periphery). The higher the rate of migration is, the faster this standardization is achieved (compare sub-scenario NC3 – with high migration, to sub-scenario NC2, Figure 4). Our results show that in all investigated situations, spatial heterogeneity tends to decrease over time, in both the local and invasive taxa, but this is a slow process. If hybridization is stopped, the level of spatial heterogeneity may still vary for a while until it reaches a value that remains constant over time (Supplementary Figure S3, supporting information).

Literature Review

We found 16 examples of hybridization during a range expansion in which the spatial distribution of introgression, or alternatively admixture components, is reported. Among those examples, two are in plants and 14 in animals. The observation made based on our simulations is supported by 14/16 of those examples (Table 2). This result means that our observation of an asymmetrical spatial gradient of introgression within interacting taxa, which was expected to be directly proportional to distance for the invasive taxon and inversely proportional to distance for the local taxon, is supported by a series of empirical examples observed in field conditions.

TABLE 2

Table 2. Examples of spatial gradients of introgression or admixture after range expansion of biological invasions.

Implications on Hybridization Between Modern Humans and Neanderthals

We illustrate the influence of range expansion on the spatial distribution of introgression by exploring the case of hybridization between modern humans and Neanderthals, which interbred during the range expansion of the former across Eurasia (Figure 5A). The ABC approach estimated the most likely values of interbreeding ranging between 0.4 and 1.3% (γ = [0.004–0.013]), with a mode estimated to be 0.8% (γ = 0.008) and a maximum value of 2% (γ = 0.02), to explain the level of introgression observed in current human populations (Figure 5B).

FIGURE 5

Figure 5. Case study on hybridization between Neanderthals and modern humans. (A) Simulated framework used in our analysis. The red area denotes the range occupied by Neanderthals. The sky-blue dot represents the starting area of modern human range expansion [chosen arbitrarily in east Africa, see Currat and Excoffier (2011)]. The two yellow dots are sampling locations where the amount of introgression is estimated. (B) Estimated value of the interbreeding rate between Neanderthals and modern humans estimated through an approximate Bayesian computation approach. The dotted line denotes the prior distribution of explored values. (C) Estimated augmentation of introgression in East Asia compared to Europe by using the most likely value of interbreeding rate (minimum and maximum bonds as well as mode). The average augmentation of simulated values is presented with a 95% confidence interval. The dotted lines delimit the range of observed augmentation values (12–20%).

For each of the modes and the two bonds of the most likely values of interbreeding success rate (i.e., γ = 0.004, 0.008, and 0.013), we computed the difference in the introgression level between East Asia and Europe using simulated samples from France and China. We thereafter call the ratio of introgression in East Asia compared to that in Europe “augmentation”. The level of augmentation decreases with increasing interbreeding rate (γ = [0.004–0.013]). The variance in augmentation is larger with the two extreme values of γ (0.004 or 0.013) compared to the estimated mode (0.008). For any likely value of interbreeding, the simulated augmentation overlaps the one observed in real data (12–20%; Figure 5C). This analysis shows that a combination of low interbreeding rate and spatial distance can explain the spatial pattern of Neanderthal introgression in humans without the need to invoke the effect of selection.

Discussion

Spatially explicit simulation is a very useful tool to evaluate the effect of range expansions on the genetic diversity of organisms (Benguigui and Arenas, 2014). In particular, this method has revealed a pattern of asymmetric introgression between a local and an invasive taxon during a biological invasion (Currat et al., 2008). This general pattern for neutral genes has been shown to be robust to the mode of dispersal of organisms, with some variation depending on the density-dependent dispersal effect (Quilodrán et al., 2019) and long-distance dispersal (Amorim et al., 2017). The emphasis of previous studies was placed on the effect of the amount of admixture between taxa, while the exploration of the spatial distribution of introgression within a taxon has been relatively neglected, especially for local organisms (Currat and Excoffier, 2004, 2005; Currat et al., 2008). Based on the same theoretical framework as previous studies, we explored in more detail the spatial patterns of introgression after different scenarios of biological invasions, both in the local and invasive taxa, as well as the persistence of these patterns over time.

Our results show that the amount of interbreeding combined with the range expansion of one or both interacting taxa has a combined influence on the spatial introgression patterns. According to our simulations, spatial heterogeneity in introgression can be found above a minimum rate of interbreeding, which depends on the demographic and migration parameters. Below this threshold of interbreeding, introgression is uniformly distributed over the whole area, but weak localized gradients may also be observed. There is also a maximum threshold to obtain a spatial pattern, above which introgression is almost complete over the whole area, as already observed by previous simulation studies (Currat and Excoffier, 2004; Currat et al., 2008). The general trend is that introgression from a local taxon to an invasive taxon in range expansion increases with the increase in the distance from the source of the expansion. The opposite trend is projected for introgression of invasive genes in the local gene pool, meaning that introgression decreases with the increase in the distance from the source of the biological invasion. The pattern is still observed when both taxa are invasive and even when hybridization is stopped after some time.

Our simulations also show that this pattern is fully correct in only the area where both taxa coexist but not necessarily in areas where only the invasive taxa occurs. Indeed, in areas not occupied by the local taxa, the amount of local genes carried by the invasive taxon can be heterogeneous. In that case, introgression is highest in zones colonized by a population that had the greatest opportunity to hybridize with the local taxon during its spread. In our simulated area, this result translates to more elevated levels of introgression along the main axis of the expansion and less introgression on the sides of the expansions (top left and bottom right corners in Figure 2, Scenario “restricted area”).

We therefore propose a double effect of range expansion and hybridization on the level of introgression. Introgression is expected to be asymmetrical “between” interacting taxa according to the observations of Currat et al. (2008) but also heterogeneous “within” a given taxon, showing a spatial gradient along the main axis of the expansion, as shown by our simulations. While our literature survey is restricted to a limited number of examples in which the spatial distribution of introgression (or equivalent proxy) is reported, the empirical evidences available suggest that the pattern revealed by our simulations could be relatively common. We thus suggest the patterns projected by our simulations as a general expectation for neutral genomic introgression, which may be further evaluated as a null hypothesis in cases of spatial range expansion and hybridization. Indeed, there are a number of factors that may cause deviation from this neutral expectation, such as selection (e.g., Johannesen et al., 2006), sex-biased mating preferences (e.g., While et al., 2015), or low hybrid fitness (e.g., Bundus et al., 2015).

Shape of the Spatial Gradient of Introgression

We show that the shape of the spatial gradient of introgression is affected by the population sizes of the invasive and local taxa as well as their migration rates. Invasive populations that are smaller in size than the local populations tend to be highly introgressed by local genes throughout the whole colonized area, reducing the length of the spatial gradient of introgression. This phenomenon occurs because a founder population of small size has a higher chance of being involved in heterospecific mating (Hubbs, 1955). In contrast, invasive populations that are larger than the local populations tend to produce a wider gradient of high intensity due to a more equal demographic balance between the two taxa. High migration rates result in wide introgression gradients due to the continuous supply of genes from conspecific populations located in the core of the expansion, which counteract the effect of interspecific gene flow (Petit and Excoffier, 2009). Similarly, long-distance dispersal can accelerate the colonization time, diminishing massive introgression into the invasive taxon (Amorim et al., 2017). At the opposite, Quilodrán et al. (2019) showed that positive density-dependent dispersal leads to a slower colonization time, which lengthens the cohabitation at the leading edge and maximizes introgression into the invasive taxon. They also showed that under the same amount of interbreeding, fast colonization times and less introgression are expected under negative density-dependent dispersal (i.e., higher dispersal probability toward lower density areas). All together, the type of the species dispersal may affect the quantity of spatial asymmetric introgression within the taxon highlighted by our simulations.

Persistence of the Spatial Gradient of Introgression

Our results suggest that heterogeneous patterns of introgression within a taxon are quite stable if the population is structured. High levels of migration may erase this heterogeneity by standardizing the spatial distribution of introgression, but this is a very slow process and the pattern may persist as long as the population is not panmictic. The stability of the spatially heterogeneous level of introgression within a taxon is even more pronounced if hybridization ends after biological invasion. In this case, the asymmetric introgression between taxa also persists with time, reaching an almost steady state.

Generation of Spatial Introgression Patterns

The process that creates gradients of introgression of local alleles into the invasive taxon could be explained as follows. First, at the wave front, the population size of an invasive taxon is small, and introgressed alleles have a relatively high chance of increasing their frequency in the population by drift. Second, the population at the front passes through a demographic increase, and alleles present at the front have thus more chance to reach high frequency than expected under a stationary stage due to the production of many descendants per individual during this growing phase. Third, the pioneers of the next colonization step are taken from the population at the front, thus dispersing introgressed alleles even further along the expansion. These three processes have been jointly described as “allele (or mutation) surfing” within theoretical frameworks (Edmonds et al., 2004; Klopfstein et al., 2006). It describes the fate of a mutation that increases in frequency by “surfing” on the wave of advance of a population expansion. Indeed, a mutation that appears in the wave front of a range expansion has more chance to surf than if it appears in the core of the expansion. This process has then been confirmed empirically under laboratory conditions (Hallatschek et al., 2007) and by historical records (Moreau et al., 2011). In the case of an expansion with hybridization, introgressed alleles are continuously introduced at the front of the invasion if the admixture rate is big enough, leading to a progressive dilution of the invasive gene pool with the local gene pool. In such conditions, introgressed local alleles thus have a relatively higher probability of “surf” over the expansion wave. In contrast, the introgression of invasive alleles in the local taxon does not benefit from this demographic and migratory “push”. This is because the local population sizes are usually at a demographic equilibrium, or are at least higher than the invasive density at the wave front. This phenomenon means that introgressed alleles from the invasive taxon always enter the local taxon at low frequency and cannot benefit from a population demographic increase to reach higher frequencies. In addition, because the invasive taxon is progressively more introgressed by local genes due to allelic surfing, the local taxon has less chance to be introgressed by invasive genes when the distance from the source of expansion increases. Indeed, far from the origin of the expansion, a large portion of the genetic pool of the invasive population is made up of introgressed genes from the local population, and thus, a backcross of introgressed alleles into its original taxon would not be detectable.

Expansion Facilitation

Theoretical studies have demonstrated that considering the dynamic of range expansion can dramatically change the interpretation of evolutionary processes and that a range expansion behaves differently than a pure demographic expansion (Excoffier et al., 2009). While both processes usually result in an increased total population size, range expansion may also have important evolutionary consequences (Excoffier and Ray, 2008). One of these consequences is that allelic surfing of deleterious mutations may induce a fitness decrease at the expanding front (Travis et al., 2007; Peischl et al., 2016), which is known as “expansion load” (Excoffier and Ray, 2008). However, the same principle applied to hybridization may produce the opposite output. Indeed, individuals at the leading edge are expected to contribute disproportionally to the genetic diversity of the invasive wave front. Those individuals are introgressed with genes from the local taxon that have had a long period of time to adapt to the environmental conditions of the area. While our simulations are purely neutral, we may speculate that the chance of increased fitness due to adaptive introgression may also increase at the expansion front. In other words, the invasiveness of an exotic species during range expansion may be enhanced by hybridization with the local organisms (i.e., “expansion facilitation”). While there are many documented cases of biological invasions that have been facilitated by adaptive introgression (e.g., McMullan et al., 2015; Whitney et al., 2015), the observation of the resulting spatial gradient would potentially be difficult because adaptive introgressed alleles may reach fixation in very few generations throughout the colonized area.

Literature Review

The available literature, where the spatial introgression pattern between taxa during expansion is reported, is scarce. The spatial pattern of admixture (which can be indicative of the levels of introgression) is reported more often than the spatial pattern of introgression, but few studies on the matter are available. However, the few empirical examples seem mostly in accordance with the pattern explored in our study. By putting into the spotlight a largely unexplored phenomenon, we hope that our study will incentivize researchers to investigate the spatial patterns of introgression in future studies regarding biological invasions or encourage them to use the already available data on such an analysis.

Implication for the Case of Hybridization Between Neanderthals and Modern Humans

A recent observation of the result of hybridization in the course of human evolution is the slightly greater Neanderthal introgression in modern Asian populations than in European populations (Meyer et al., 2012; Wall et al., 2013; Prüfer et al., 2017; Villanea and Schraiber, 2019). Different levels of purifying selection have been invoked to explain this different level of introgression (Sankararaman et al., 2014; Villanea and Schraiber, 2019). This observation has also motivated a debate about the number of admixture events between Neanderthals and modern humans. Either single or multiple hybridization pulses have been proposed, both possibly linked with further dilution of introgression in Europe due to interbreeding with a “ghost” unadmixed population (Green et al., 2010; Sankararaman et al., 2012; Vernot and Akey, 2015; Villanea and Schraiber, 2019). Our simulation framework does not need different levels of selection or complex demographic scenarios of dilution to explain the higher introgression level in East Asia than in Europe. This phenomenon may simply be explained by the range expansion of modern humans, acting as an invasive taxon and colonizing an area already occupied by Neanderthals. The difference in introgression levels in East Asia and Europe would be explained by the longer distance where interbreeding could have occurred in Asia than in Europe. If modern humans crossed a large distribution zone of Neanderthals in Western Asia, they would have many opportunities to incorporate Neanderthal genes through interbreeding. These genes would then have been carried until East Asia.

Because the 12–20% difference in introgression between the current populations in Western and Eastern Eurasia (Meyer et al., 2012; Wall et al., 2013; Prüfer et al., 2017; Villanea and Schraiber, 2019) is a genome-wide average estimate, we believe that the consequences of the range expansion of modern humans outside Africa is a more parsimonious explanation than a differential selection pressure in these two areas acting on the genome as a whole. This expansion effect is not exclusive to selection pressure at certain loci, such as genes related to the immune system (Mendez et al., 2012; Quach et al., 2016), skin pigmentation (Vernot and Akey, 2014) and adaptation to altitude (Huerta-Sánchez et al., 2014). However, a regional selective pressure that would have favored Neanderthal introgression on the whole genome in East Asia seems to us more difficult to understand, notably because the Neanderthal was probably more adapted to the European climate, where its specific traits appeared, than to the East Asian climate.

By following the same simulation approach of Currat and Excoffier (2011), we show that the most likely value of interbreeding rate also results in an increased introgression level in East Asia that overlaps with the observed values. The value of interbreeding estimated by our simulation (< 2%) is similar to that estimated by Currat and Excoffier (2011) but slightly lower than the value < 3% estimated by Excoffier et al. (2014). This difference is because the later study included the higher value of introgression in East Asia than in Europe as a prior for the estimation of the interbreeding rate, which was not included in our simulations or in Currat and Excoffier (2011). Our estimation of interbreeding is therefore consistent with previous estimates, which suggests a strong avoidance of interspecific mating, low hybrid fitness, or both. Our results are consistent with multiple events of admixture, which is in line with recent studies (Sankararaman et al., 2012; Wall et al., 2013; Kuhlwilm et al., 2016; Vernot et al., 2016; Villanea and Schraiber, 2019). By using the same framework, Currat and Excoffier (2011) have already favored this scenario over a single or few events of admixture in a restricted area. These events of admixture were potentially distributed in different areas of Europe and Asia, and our results suggest that they were slightly more numerous toward the East than the West to explain the difference in introgression level across Eurasia.

Current knowledge on the introgression of Denisovan genes into modern humans is not incompatible with recurrent episodes of hybridization in East Asia or Oceania during the expansion of the latter outside Africa (Jacobs et al., 2019). However, current data on the geographical area occupied by the Denisovans is extremely limited (Reich et al., 2010; Chen et al., 2019) and various hypotheses could explain the fact that remains of Denisovans are located more than 8,000 km away from the zone where the maximum introgression in modern populations has been found, around New Guinea (Reich et al., 2011).

Conclusion

Our simulations show spatially heterogeneous introgression patterns between and within interacting taxa resulting from a biological invasion with hybridization in structured populations. Empirical observations of these patterns could thus be indicative of a biological invasion with a limited amount of hybridization between taxa, in which the orientation of the patterns could help recover the source of the invasion. However, the interpretation of patterns in terms of population dynamics is not trivial, as similar configurations could be obtained by different means (equifinality). Contrary to demographic effects, introgression during a biological invasion does not affect all loci but only a limited number of them (those introgressed), which makes it difficult to distinguish introgression from selection (Klopfstein et al., 2006; Currat et al., 2008; Excoffier and Ray, 2008; Excoffier et al., 2009). This finding means that spatially heterogeneous levels of introgression observed within a taxon after a biological invasion do not necessarily imply selective effects, and we argue that the consequences of the demographic history of species should be considered before invoking selection

Although range expansion and hybridization have occurred in the history of various species (Mitchell et al., 2019) and are potentially occurring at an increasing rate due to climate change (Oyler-McCance et al., 2016), investigation of their genomic consequences deserves more attention. We believe that the theoretical framework we used to study the genomic effect of biological invasion provides a useful tool in this field.

Data Availability Statement

All the data necessary to repeat the simulations are described in the main text and can be reproduced by using SPLATCHE3 (www.splatche.com/splatche3). The setting files and custom scripts are deposited in Zenodo: http://doi.org/10.5281/zenodo.4034797.

Author Contributions

CQ and MC conceived, designed, and coordinated the study. CQ carried out simulations and analyses. AT compiled and critically reviewed the data. All authors interpreted the results, contributed to the writing of the manuscript and gave final approval for publication.

Funding

This study was financed by grants from the Swiss National Science Foundation no 31003A_182577 to MC and P400PB_183930 to CQ.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

All computations were performed using the High-Performance Computing (HPC) cluster at baobab.unige.ch. We thank Sonya Clegg for her careful reading of a previous version of our manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.569620/full#supplementary-material

References

Ammerman, A. J., and Cavalli-Sforza, L. L. (1973). “The explanation of culture change: models in prehistory,” in The Explanation of Culture Change, ed. C. Renfrew (London: Duckworth), 343–357.