Molecular and ecological signatures of an expanding hybrid zone

Abstract Many species are currently changing their distributions and subsequently form sympatric zones with hybridization between formerly allopatric species as one possible consequence. The damselfly Ischnura elegans has recently expanded south into the range of its ecologically and morphologically similar sister species Ischnura graellsii. Molecular work shows ongoing introgression between these species, but the extent to which this species mixing is modulated by ecological niche use is not known. Here, we (1) conduct a detailed population genetic analysis based on molecular markers and (2) model the ecological niche use of both species in allopatric and sympatric regions. Population genetic analyses showed chronic introgression between I. elegans and I. graellsii across a wide part of Spain, and admixture analysis corroborated this, showing that the majority of I. elegans from the sympatric zone could not be assigned to either the I. elegans or I. graellsii species cluster. Niche modeling demonstrated that I. elegans has modified its environmental niche following hybridization and genetic introgression with I. graellsii, making niche space of introgressed I. elegans populations more similar to I. graellsii. Taken together, this corroborates the view that adaptive introgression has moved genes from I. graellsii into I. elegans and that this process is enabling Spanish I. elegans to occupy a novel niche, further facilitating its expansion. Our results add to the growing evidence that hybridization can play an important and creative role in the adaptive evolution of animals.

indicate that such de novo sympatric zones are commonly being created across a wide range of species, including invertebrates (Mallet, Wynne, & Thomas, 2011), vertebrates (Kelly, Whiteley, & Tallmon, 2010), and plants (Gómez, González-Megías, Lorite, Abdelaziz, & Perfectti, 2015), efforts to explore the long-term consequences of secondary sympatry on the genetic diversity and realized niche use are largely lacking. Indeed, recent research suggests that the known examples merely represent the tip of the iceberg and that future work will show that introgressive hybridization following range shifts is much more widespread than has been anticipated.
The outcome of species mixing via introgressive hybridization directly affects the genetic diversity of populations and is hence a critical determinant of the long-term evolutionary population processes and patterns (Baack & Rieseberg, 2007;Ryan, Johnson, & Fitzpatrick, 2009). While such species mixing can be an important source of new genetic variation (Kays, Curtis, & Kirchman, 2009;Seehausen, 2004), it can also lead to genetic swamping and thus the genetic extinction of a population or even a species (Rhymer & Simberloff, 1996). Hybridization can also result in geographically limited hybrid zones with low levels of gene flow across otherwise stable population boundaries or may lead to reinforcement of reproductive barriers due to low-hybrid fitness (Arnold, 1997). The majority of work thus far on these processes has focussed on plants (Rieseberg & Brunsfeld, 2012) and fishes (Keller et al., 2013;Seehausen, 2004), because hybridization in these groups is particularly common, but to a much lesser extent on invertebrates, with the exception of a few well-documented cases of hybridizing non-native species (Sánchez-Guillén, Cordoba-Aguilar, Hansson, Ott, & Wellenreuther, 2016).
In addition to effects on the genetic diversity, introgressive hybridization can also affect the ecological population dynamics of species via adaptive trait introgression, that is, the introduction of new alleles that increase the competitiveness of individuals (Arnold, 2004). Lewontin and Birch (1966) proposed that such adaptive trait introgression can facilitate the rate of adaptive evolution and gathered empirical evidence for this by studying the fruit fly Drosophila tryoni in Australia. This fly species has greatly extended its geographic range into warmer areas over the last 100 years, and Lewontin and Birch (1966) were able to demonstrate that the increase in upper temperature tolerance was caused by introgressive hybridization with Drosophila neohumeralis. Since then, it has become increasingly established that even though the initial hybridization is often deleterious, subsequent introgression can provide a unique opportunity for species to tap into genetic variation present in a closely related species and potentially take advantage of beneficial alleles (Baack & Rieseberg, 2007;Huerta-Sanchez et al., 2014;Mallet, 2005;Rieseberg, 2009;Seehausen, 2004). It is noteworthy that alleles with widespread effects can spread particularly easy across species boundaries when species isolating barriers are incomplete (Coyne & Orr, 2004) to produce the signature of a trans-species selective sweep (Brand, Kingan, Wu, & Garrigan, 2013).
Here, we study the hybrid zone of the blue-tailed damselfly Ischnura elegans (Odonata: Coenagrionidae) in southern Europe, where it is overlapping with its sister species, the Iberian bluetail Ischnura graellsii. These species show only incomplete reproductive barriers, enabling chronic introgression (Sánchez-Guillén et al., 2012). Specifically, we investigate the molecular and ecological signatures of this hybrid zone to (1) quantify the pattern and extent of molecular population divergence and introgression across populations inside and outside the hybrid zone, and 2) test if trait introgression has affected the ecological niche use of both species and its hybrids.

| Species natural history and ecology
Ischnura elegans has become an eco-evolutionary model species because of its enigmatic female-limited genetic color polymorphism and has been studied extensively for its population morph dynamics (Cordero & Andrés, 1996;, the genetics of color Chauhan et al., 2014;Sánchez-Guillén, Van Gossum, & Cordero-Rivera, 2005), and the processes maintaining the color polymorphism in nature (Sánchez-Guillén et al., 2017;Svensson, Abbott, & Hardling, 2005;Takahashi, Yoshimura, Morita, & Watanabe, 2010). Ischnura elegans and I. graellsii are closely related (Sánchez-Guillén et al., under revision) and share many phenotypic traits, including many preference traits for habitats. Ischnura elegans extends from Ireland to the Mediterranean and Japan, while I. graellsii extends across the western Mediterranean area including Iberia and Maghreb. Both species are common along running and standing water bodies, and both are tolerant to brackish waters . Whenever the habitat is suitable, both species form highly abundant populations and are commonly the dominant species in ponds and lakes. Previous studies have demonstrated a lack of isolation by distance in I. elegans and I. graellsii populations .
In the last years, I. elegans has also received attention because it is undergoing a range expansion in northern Europe and in Spain (Lancaster et al., 2016;. Despite its good dispersal ability, the distribution of I. elegans in Spain is discontinuous: It is common in north-central Spain and along the Mediterranean coast, from Pyrenees to Murcia, where it is the dominant species. However, several populations of I. elegans were found in the early 1990s in north-western Spain, more than 400 km away of the central Spanish populations (Ocharan Larrondo, 1987). Central Spain Earliest evidence for this hybridization came from work by Monetti, Sánchez-Guillén, and Cordero-Rivera (2002) who showed that allopatric I. graellsii are smaller in body size and have narrower and shorter wings and shorter tibiae, whereas I. elegans have a narrower space between the branches of each cercus and a greater distance between the branches of each paraproct. In sympatry, however, these morphological traits were intermediate, strongly suggesting that introgressive hybridization affects the overall phenotype (2002).
Our preliminary work on the molecular population signatures in this hybrid zone showed evidence for extensive and asymmetric introgression (because of strong asymmetry in reproductive barriers, see Sánchez-Guillén et al., 2012), with I. elegans replacing I. graellsii from many areas in Spain .

| DNA extraction and microsatellite genotyping
Genomic DNA for microsatellite typing was extracted from 14 populations (208 I. elegans individuals and 29 I. graellsii, Table 1) using a standard phenol/chloroform-isoamyl alcohol protocol (Sambrook, F I G U R E 1 Map of sampled populations. Diamonds denote allopatric Ischnura elegans; squares sympatric, introgressed I. elegans; upside triangles allopatric Ischnura graellsii; and downside triangles: I. graellsii sympatric with I. elegans. Solid symbols correspond to populations sampled for molecular analyses

| Microsatellite data analyses
Twelve microsatellites were used to assess the population level of genetic diversity in terms of allele frequencies, observed heterozygosity (H O ), and allelic richness in the program GenAlEx 6.503 (Peakall & Smouse, 2012). Pairwise population differentiation was assessed by calculating F ST (Weir & Cockerham, 1984) and Dest (Jost, 2008) in the R v3.1.2 (R Development Core Team 2016) package diveRsity (Keenan, McGinnity, Cross, Crozier, & Prodöhl, 2013) using 1,000 bootstraps for the calculation of 95% confidence intervals. Structure version 2.2.3 (Pritchard, Stephens, & Donnelly, 2000) was used to quantify admixture proportions of each population for F 1 -hybrids and backcrosses (Beaumont et al., 2001;Hansen & Mensberg, 2009;Sanz, Araguas, Fernández, Vera, & García-Marín, 2009;Vähä & Primmer, 2006). We conducted admixture analyses in Structure to assign Ischnura individuals (N = 156) from the parapatric (Vigueirat and Menorca) and sympatric populations (Doniños, Arreo, Baldajo, Alfaro, Estanyo de Europa, Amposta, and Marjal del Moro) into each of two genetic clusters, one representing I. graellsii genotypes and one I. elegans. We used the "prior population information" option to facilitate the clustering process of the reference individuals (i.e., allopatric I. elegans, and allopatric I. graellsii) and to calculate the admixture proportions (and ±90% credible regions) of each individual in parapatric and sympatric populations. This approach was used to measure the degree of introgression of I. graellsii genetic material into the genome of I. elegans in Spain. The model was run for K = 2, where each cluster corresponded (to a very high extent) to I. elegans and I. graellsii, respectively. We used the "population flag" option to exclude parapatric and sympatric I. elegans as reference individuals, which implied that the clustering process was based on only I. graellsii samples and I. elegans samples collected outside of Spain (i.e., outside the hybrid zone). The model was run for 100,000 MCMC replicates, after an initial burn-in period of 50,000 replicates, using the admixture model and correlated frequencies for one iteration.
TA B L E 1 Population and country names that were included in the microsatellite work. Species are denoted as Ie for I. elegans, and Ig for I. graellsii. For each population, the sampling year (Year), Ecology (sympatric or allopatric), number of individual samples (N), number of alleles observed across all loci (Alleles#), the observed heterozygosity (H O ), and expected heterozygosity (H e ) are given Simulated hybrid genotypes were generated to compare them with the levels of admixture detected in natural populations. That comparison will allow us to confirm the pattern of hybridization in each population. To generate simulated genotypes of hybrids and backcrosses, we applied the program Hybrid-Lab (Nielsen, Bach, & Kotlicki, 2006) using the genotypes of 29 individuals of I. graellsii and 52 genotypes of I. elegans collected outside of Spain as initial genotypes. We generated 50 genotypes of each of the following crosses:  were constructed setting most parameters to default ("Auto features," convergence = 10 −5 , maximum number of iterations = 500, background = 10,000 points) and varying the prevalence (0.5, 0.6, and 0.7) and the regularization value β (1, 2, and 3) to find which combination generated the best outcomes (estimated using the highest area under the curve, or AUC value) while minimizing the number of model parameters, as well as producing "closed," bellshaped response curves.

| Niche analyses
We subsequently performed niche identity and background similarity tests, comparing the actual ecological niche models with random models generated from 100 pseudoreplicate data sets (Warren, Glor, & Turelli, 2008) using Schoener's D and a modified Hellinger distance's I statistics as implemented in the package "phyloclim" (Heibl & Calenge, 2014) for R v3.1.2 (R Development Core Team 2016). The values of the two statistics vary between 0 (no niche overlap) and 1 (identical ecological niche). The background similarity test is particularly pertinent in this study because we seek to test whether the niches from the two allopatric populations are more different than would be expected given the underlying environmental differences between the regions in which they occur.

| Molecular population structure
Microsatellite analyses showed that the populations exhibited a high degree of genetic variation, as shown by the pronounced genetic diversity at each locus (Table 1). The total number of alleles over all loci per population ranged between 32 (Menorca) and 112 (Arreo).  (Table 2). This is consistent with the idea that this zone is characterized by ongoing introgression between the species. Further, populations close to the Mediterranean appear to be more closely related to I. graellsii compared with populations from central Spain (Table 2) Figure 2a). The introgressed populations were classified as a mixture of these two clusters (Table 3, Figure 2a).

Most individuals of the European I. elegans outside Spain (72%) and
I. graellsii (72%) were assigned with a certainty of at least 75% to each of these clusters. However, only 54% of the parapatric and 8% of the I. elegans of the sympatric populations were assigned with a certainty of at least 75% to the I. elegans cluster (Table 3; Figure 2a).
The rest of the sympatric I. elegans were intermediate between the two clusters, suggesting a significant degree of introgressed I. graellsii alleles (Table 3; Figure 2a).
The admixture proportions of the artificial hybrids and backcrosses ranged between 19 and 89% (Table 4; Figure 2b). The F 1 showed admixture proportions to I. elegans between 45% and 68%; the three I. elegans backcrosses (1-3EB) between 54% and 89%, and the three I. graellsii backcrosses (1-3GB) between 19% and 61% (Table 4; Figure 2b). Three conservative assignment groups were defined based on the results from the artificial hybrids: F 1 and BCE-1 and BCG-1 ranging between 55% and 50%; backcrosses with I. elegans ≥75; and backcrosses (only BCG-2 and BCG-3) with I. graellsii <25%. The parapatric and sympatric I. elegans genotypes (Table 3; The genotypic linkage disequilibrium analyses between pairs of loci showed significant associations for allopatric I. elegans populations. Of the 55 pairs of loci that could be tested, four deviated significantly from equilibrium (Table S1). Furthermore, when testing the sympatric I. elegans populations 10 of the 66 pairs of loci deviated significantly from equilibrium (Table S1). In contrast, no significant deviations from equilibrium were detected for allopatric I. graellsii populations (21 pairs of loci were tested) (Table S1).

| Niche analyses
We tested if the environmental niche models of the two species were more different than expected by chance (identity test), and found that all pairwise comparisons of sympatric and allopatric populations TA B L E 3 Summary of the results from the K = 2 admixture models in Structure for I. elegans and I. graellsii showing the number of individuals per population assigned to different admixture proportion categories (based on individual assignment to the red cluster in Figure 2) of I. elegans and I. graellsii were less similar than expected by chance (Table 5). This indicates that all populations (I. elegans sympatric and allopatric and I. graellsii sympatric and allopatric) are found to occupy different niches, despite often occurring in close geographic proximity.
We also carried out background niche tests to determine whether the species' ranges were more different from one another than expected based on the environmental background differences (Figure 3a). Our background niche analyses using Hellinger's I demonstrated that both species from their allopatric distributions show decreased niche similarity (Figure 3b). Interestingly, I. elegans in the sympatric portion shows decreased niche similarity with I. elegans in the allopatric portion ( Figure 3c) and I. graellsii in the allopatric portion of their range (Figure 3d) as both similarity scores were lower than the generated null hypotheses, while I. graellsii from the sympatric portion was not statistically lower or higher (than what would be expected by random) compared with I. graellsii from allopatry (Figure 3e) or I. elegans from allopatry ( Figure 3f).

Our background niche analyses also demonstrated that I. elegans
in the sympatric portion showed increased niche similarity with I. graellsii in the sympatric portion of their range (Figure 3g), as indicated by the similarity score (red line) being higher than the null hypothesis of niche equivalency. Together this indicates that it has been a niche displacement in I. elegans resulting in that sympatric populations of I. elegans now share the same niche space as I. graellsii, but not anymore with I. elegans from allopatry. However, we found no evidence for niche displacement in the sympatric I. graellsii populations, which still share a niche with the allopatric I. graellsii populations.

| D ISCUSS I ON
It is recognized that hybrid zones vary greatly in their structures depending on the degree of the molecular and ecological differentiation between the introgressing species and affect associated levels of dispersal (Buggs, 2007;Glotzbecker, Walters, & Blum, 2016). Our molecular results suggest that introgressive hybridization between I. elegans and I. graellsii is chronic, and our niche modeling indicates TA B L E 5 Summary of niche identity and background similarity tests. "Divergence" indicates that the compared populations show significant divergence (overlap is less than expected), while "Conservatism" indicates niche conservatism (overlap values are more similar than expected). "**" indicates a significant (p < .01?), and "n.s." no significant, difference between expected and observed overlap   F I G U R E 3 Niche analyses. Results of background similarity test, showing modified Hellinger distance I. The observed similarity between niches is indicated with the red lines with Hellinger values on top, while histograms indicate the null distribution of ecological niche distance generated randomly. Hellinger distance I ranges from 0 (complete different) to 1 (identical). (a) Map of compared region niches (I. elegans from allopatry, I. elegans from sympatry, I. graellsii from allopatry, and I. graellsii from sympatry). Black unequal symbol indicates that niches are more different that by change; gray equal symbols indicate no differences between niches, and white equal symbol indicates that niches are more equal than by change. (b-d) The observed similarities were lower than their respective null distribution for random niche models.
(e and f) The observed similarities were similar that their respective null distribution for random niche models. (g) The observed similarities were higher than their respective null distribution for random niche models  (Todesco et al., 2016).
Adaptation to local environments can be achieved by selection on standing genetic variation (i.e., variation present in the species at the time of the environmental change) or de novo mutations.
Another less explored path for species to acquire beneficial alleles is through introgression via interspecific hybridization (Arnold, 1997;Glotzbecker et al., 2016;Rieseberg & Wendel, 1993;Whitney, Randell, & Rieseberg, 2006). The outcomes of such introgressive hybridization are manifold and include the fusion of species, genetic swamping of one species by another, elicit reinforcement of reproductive isolation between incompletely isolated species, transfer of genetic material between species, potentially facilitating their adaptive evolution, and ultimately, the origin of new species (Seehausen, 2004;Zemanova, Knop, & Heckel, 2017). Our data suggest that the introgressive hybridization between I. elegans and I. graellsii is potentially providing a new source of alleles for the populations in the hybrid zone and that this may enable their contemporary and possibly also future spread. Introgression can also happen without an obvious impact on the genetic diversity if the recombinant genotypes are retained in a narrow contact zone between two species via a balance between dispersal into the contact zone and selection against hybrids. However, in the current case, it seems that the hybrid zone is dynamic and moving, and this means that the genetic impact will not be confined to a small geographic area. This indicates that the genetic exchange between I. elegans and I. graellsii may have widespread and prolonged consequences for these Ischnura species in southern Europe, and that the final impact of the introgressive hybridization on the species diversity and maintenance is difficult to predict with certainty.
The view of prolonged hybridization in the evolution of biological diversity has changed over time. Some thought of hybrids as the raw materials of evolution and a creative source of functional novelty (Buggs, 2007;Rieseberg & Wendel, 1993). Others, however, thought that hybridization presented an evolutionary dead end (Mayr, 1963).
These contrasting views were based on the fact that at least 25% of plant species, but only 10% of animal species, mostly phylogenetically recent species, are involved in hybridization and potential introgression with other species (Mallet, 2005). Consequently, many plant studies viewed hybridization as an important source of new genetic variation and a frequent component of the evolutionary history of many plant species (Harrison, 1993), whereas hybridization within the animal kingdom was viewed as a rare occurrence (Arnold & Hodges, 1995). Early work (Anderson, 1949;Anderson & Stebbins, 1954;Lewontin & Birch, 1966) emphasized the importance of reshuffling segregating genetic variation via hybridization. Since then, some empirical studies have emerged in support of this view, and these indicate that hybridization provides an important source of genetic variation on which selection can act and that its adaptive role is even more widespread, among both plants and animals, than previously believed (Dowling & Secor, 1997;Hedrick, 2013;Rieseberg & Wendel, 1993;Rius & Darling, 2014). Awareness is also increasing that this mixing of gene pools can in some cases lead to the potential loss or displacement of genotypically distinct species, something which can be especially problematic for rare organisms coming into contact with more abundant ones (Zemanova et al., 2017). Our study adds a new empirical example from a nonmodel animal species to this idea and shows that genetic introgression may be facilitating the range expansion of some hybridizing species via niche expansion.
In conclusion, the molecular data and the niche modeling together suggest adaptive trait introgression from the native I. graellsii into the expanding I. elegans, and that this is enabling the niche shift of introgressed I. elegans (sensu Arnold, 2004). ability to respond to a changing climate. The ongoing interspecific gene flow with I. graellsii may increase the limited adaptive potential that results from standing genetic variation and mutation alone, enabling a quicker demographic recovery in response to changing environments. Taken together, our results corroborate the view that hybridization can play an important and creative role in adaptive evolution.